# User:Ageary/Yilmaz 2

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

## 8.0 Introduction

Strong lateral velocity variations associated with complex overburden structures require earth imaging in depth. Examples of complex overburden include diapiric structures formed by salt tectonics, imbricate structures formed by overthrust tectonics and irregular water-bottom topography. All three are characterized as structure-dependent lateral velocity variations. There also exist structure-independent lateral velocity variations, often associated with facies changes; for instance, changes in lithology from shale to sandstone to carbonate induce lateral changes in acoustic impedance.

Figure 8.0-1a shows a field data example of a diapiric structure associated with salt tectonics. Note how the base-salt (event B) and subsalt reflections at about 2 s are pulled up in the middle of the section. Time migration yields an inaccurate image of the base-salt (Figure 8.0-1b).

Figure 8.0-2a is a field data example of an imbricate structure associated with overthrust tectonics. Time migration is adequate for imaging a target within the imbricate structure itself, but would not yield a correct image of a target below it (Figure 8.0-2b).

Refer to the field data example for a case of irregular water-bottom topography shown in Figure 8.0-3a. Note the false structures along the unconformity T because of irregular water-bottom topography. These are especially noticed below midpoint locations A, B, and C. Sagging of the reflection that is associated with the unconformity at these locations is attributed to the low-velocity overburden, here, the water layer. Once again, time migration does not provide the solution for such a complex overburden as shown in Figure 8.0-3b. Specifically, distortions along the unconformity still exist.

Figure 8.0-1  (a) A CMP stacked section over a salt-dome structure, (b) time migration. Note the similarity of the overmigration at the base-salt event in this section with that seen in the time-migrated section of the synthetic data in Figure 8.2-7.

Earth imaging in depth is achieved by depth migration. The following is a cause-and-effect relation between various factors with regard to depth migration:

1. Complex overburden structures often give rise to strong lateral velocity variations. In the presence of strong lateral velocity variations, an earth image in time derived from time migration is not accurate; and thus, it is imperative to obtain an earth image in depth by depth migration.
2. Strong lateral velocity variations cause significant ray bending at layer boundaries.
3. This then gives rise to nonhyperbolic behavior of reflection times on CMP gathers that correspond to layer boundaries below a complex overburden structure.
4. As a result, amplitudes and traveltimes associated with the reflection events with nonhyperbolic moveout are distorted during conventional CMP stacking which is based on the hypberbolic move-out assumption.
5. This causes CMP stack to depart from an ideal zero-offset wavefield. Therefore, when depth migration is needed, in principle, it must be done before stack and not after.
6. Finally, complex overburden structures often exhibit three-dimensional (3-D) behavior. Therefore, when depth migration is needed, again in principle, it must be done not only before stack but also in three dimensions.

Compare the results of depth imaging in 2-D and 3-D using post- and prestack data shown in Figures 8.0-4 through 8.0-7. Note the improved imaging of the fault planes and the reflectors within the faulted zone attained by combining the advantages of prestack migration and 3-D migration. Figure 8.0-4 shows 2-D poststack depth migration of the 3-D DMO stack in Figure 7.0-16a derived from common-cell gathers along an inline from a 3-D marine survey. Velocities vary in the lateral direction across the steep faults in the middle portion of the image. The imaging within the faulted zone can be improved by prestack depth migration as shown in Figure 8.0-5. Nevertheless, the 3-D behavior of the fault planes requires 3-D imaging (Figure 8.0-6). The complete solution to the imaging problem in the presence of lateral velocity variations and the 3-D behavior of reflector geometries, of course, is 3-D prestack depth migration (Figure 8.0-7). Since lateral velocity variations can change the polarity of fault-plane reflections, the image sections in Figures 8.0-4 through 8.0-7 are displayed both with normal and reverse polarities.

Although the migration methods discussed in migration are based on a layered media assumption, simple modifications of the basic algorithms make them accurate for situations with mild lateral velocity variations. For example, rms velocities can be varied laterally in Kirchhoff migration. In the finite-difference method, as long as lateral velocity variations are mild, the thin-lens term can be dropped (Section D.3), and the velocity function used in the diffraction term can be varied laterally. In the frequency-wavenumber methods such as Stolt migration, lateral velocity variations are accommodated by varying the stretch factor between 0 and 1. Even when velocity varies, the output from these three methods is still a time section, thus, the term time migration.

The situation is different when strong lateral velocity variations are encountered. In that case, simple algorithmic modifications no longer provide adequate accuracy, and depth migration [1] [2], rather than time migration, must be done. While both migration types use a diffraction term for collapsing energy along a diffraction hyperbola to its apex, only depth migration algorithms implement the additional thin-lens term that explicitly accounts for lateral velocity variations (Section D.4). Unlike time migration, the output from the migration algorithms that include the thin-lens term is a depth section, thus the term depth migration.

Sometimes, although not so often, the complex overburden can be defined by a single layer and the boundary between the overburden and the substratum can be determined in the form of an irregular interface with a significant velocity contrast. In that case, the layer replacement method described in layer replacement, in lieu of depth migration, can be used to remove the deleterious effect of the overburden on the geometry of the underlying reflections.

Time migration requires an rms velocity field, whereas depth migration requires an interval velocity-depth model. A velocity-depth model usually is defined by two sets of parameters — layer velocities and reflector geometries. While an rms velocity field does not contain discontinuities, an interval velocity-depth model can include discontinuities associated with layer boundaries.

A velocity-depth model is the seismic representation of an earth model in depth. An earth model and the earth image created from it are an inseparable pair of products of seismic inversion. To obtain an earth image in depth, one has to first estimate an accurate earth model in depth. We shall defer estimation of earth models in depth to the next chapter where we shall demonstrate the use of various inversion methods and depth migration itself to estimate the earth model parameters. Aside from earth modeling and imaging in depth, depth migration also is used to verify and update velocity-depth models. In this chapter, we shall be concerned only with depth migration as a tool for earth imaging in depth.

Historically, depth migration was used in an iterative fashion. Starting with an initial velocity-depth model, depth migration is performed and results are interpreted for updating layer velocities and reflector geometries. Using the updated velocity-depth model, depth migration is repeated until such time as what is inferred from depth migration as the velocity-depth model matches with the velocity-depth model input to depth migration. To achieve rapid convergence to a final velocity-depth model, only one set of parameters — usually, reflector geometries, are altered from one iteration to the next. We shall review the iterative depth migration procedure for zero-offset, CMP-stacked and prestack data, and draw conclusions as to when such a procedure yields an acceptable image in depth (2-D poststack depth migration).

In this chapter, we shall review two methods of two-dimensional (2-D) prestack depth migration (2-D prestack depth migration). In 2-D migration, we assume that the seismic line is along the dip direction and that the recorded wavefield is two dimensional. Shot-geophone migration, which is based on downward extrapolation of common-shot and common-receiver gathers using the double-square root equation (Section D.1), focuses primary reflection energy to zero-offset. The migrated section is obtained by retaining the zero-offset traces and abandoning the nonzero-offset traces. Shot-profile migration is based on migrating each common-shot gather, individually. In this method, the migrated section is obtained by sorting the migrated common-shot gathers into common-receiver gathers and summing the traces in each receiver gather.

We shall review 3-D post- and prestack depth migrations in 3-D poststack depth migration and 3-D prestack depth migration, respectively. While 3-D poststack depth migration can be performed using a wide variety of algorithms based on finite-difference (Sections G.1 and G.2), frequency-wavenumber (Sections G.3 and G.4), and Kirchhoff integral solutions (Section H.1) to the 3-D scalar wave equation, 3-D prestack depth migration often is done using the Kirchhoff integral (Section H.1) or the eikonal equation (Section H.2). This is because Kirchhoff migration is conveniently adaptable to source-receiver geometry irregularities in 3-D seismic surveys. Topographic irregularities also are conveniently accommodated in Kirchhoff migration.

The output of prestack depth migration consists of image gathers, which may be likened to moveout-corrected CMP gathers with the vertical axis in depth. Image gathers, however, consist of traces in their migrated positions. A stack of image gathers represents the earth image in depth obtained from prestack depth migration. If the velocity-depth model used in prestack depth migration is correct, then, events on an image gather would exhibit a flat character with no moveout. An erroneously too low or too high velocity would cause a residual moveout on the image gather. The initial velocity-depth model can be updated by analyzing this residual moveout and correcting for it (model updating).

### Lateral velocity variations

Lateral velocity variations often are associated with steep dips. Hence, a depth migration algorithm must not only handle lateral velocity variations but also must image steeply dipping events accurately. The steep-dip implicit and explicit frequency-space migration algorithms described in frequency-space migration in practice are particularly suitable to accommodate lateral velocity variations. In the case of the implicit schemes, the action of the thin-lens term, which accounts for lateral velocity variations, is achieved by a complex multiplication of the wavefield in the frequency-space domain with a velocity-dependent exponential term (Sections D.3 and D.4). In the case of the explicit schemes, lateral velocity variations are accounted for by designing a velocity-dependent, laterally varying extrapolation filter and convolving it with the wavefield in the frequency-space domain (Section D.5).

The problem of lateral velocity variations will be studied using a point diffractor buried in a medium with five different types of velocity-depth models. The first velocity-depth model is shown in Figure 8.0-8. The corresponding zero-offset section consists of an exact diffraction hyperbola. Therefore, only the diffraction term is needed in imaging the scatterer (Section D.3). The surface projection of the point diffractor (the source) coincides with the location of CMP 240, indicated by the vertical arrow, and is aligned with the apex of the diffraction hyperbola. After time migration, the diffraction hyperbola is collapsed to its apex which, in this case, is coincident with the CMP location of the point diffractor.

Consider what happens when the point diffractor is lowered into the second layer as shown in Figure 8.0-9. Raypaths from the scatterer to the surface are bent at the interface between the first and second layers according to Snell’s law of refraction. The zero-offset section is now approximately hyperbolic. From normal moveout, we recall that in a horizontally layered earth model, travel-times are governed by the hyperbolic moveout equation. However, moveout is hyperbolic only within the small-spread approximation. The velocity associated with this approximate hyperbola is the vertical rms velocity down to the diffractor. Suppose the velocity-depth model of Figure 8.0-9 is replaced with that of Figure 8.0-10, where the velocity of the first layer now corresponds to the rms velocity of the diffractor of Figure 8.0-9 (1790 m/s). The zero-offset section associated with this new model is an exact hyperbola (Figure 8.0-10).

The traveltime trajectory in the zero-offset section derived from the original model (Figure 8.0-9) differs negligibly from the hyperbolic traveltime trajectory in the section associated with the replacement model (Figure 8.0-10) only at the far flanks. Therefore, it is reasonable to assume that the zero-offset traveltime trajectory for a point diffractor in a horizontally layered earth model is a hyperbola. The apex of this approximate hyperbola in Figure 8.0-9 coincides with the surface projection of the diffractor, indicated by the arrow. Therefore, only time migration is needed to image a diffractor that is buried in a horizontally layered earth model. This time migration can be performed either by Kirchhoff summation, which uses the rms velocity, or by frequency-space or frequency-wavenumber methods, which honor raypath bendings at the interfaces associated with horizontal layers.

Now consider the point diffractor situated in the third layer as shown in Figure 8.0-11. Now we no longer have even an approximate hyperbolic diffraction response. The traveltime trajectory is skewed so that the apex A does not coincide with the lateral position B of the diffractor. As expected, time migration partially focuses the energy toward its apex A, which is to the left of the actual lateral position of diffractor B.

To properly focus the energy and place it laterally to its true subsurface position B, depth migration must be performed as shown in Figure 8.0-11. The depth-migrated image is aligned with the true subsurface position B. This lateral positioning is accomplished by the thin-lens term. The amount of lateral shift is the lateral distance AB between the apex A of the skewed diffraction traveltime trajectory and the actual location B of the diffractor. This shift depends on the amount of ray bending that occurs at the interfaces above the point diffractor.

From Figure 8.0-11, note that the apex of the skewed traveltime trajectory A coincides with the surface position of the ray that emerges vertically. This special ray is called the image ray, and it was first recognized by Hubral [3]. The image ray associated with the point diffractor in Figure 8.0-11 is roughly at midpoint 200. The diffractor itself is located beneath midpoint 240. Therefore, the lateral shift is equivalent to 40 midpoints.

There is no lateral shift for the horizontally layered earth model (Figure 8.0-9), since there is no lateral velocity variation. The image ray in this case emerges at the surface location of CMP 240 coincident with that of the diffractor. For a mild to moderate lateral velocity variation, as in Figure 8.0-12, the lateral shift is less than 10 midpoints. For some objectives, this small lateral shift and the complete focusing may not be critical; therefore, time migration may be acceptable in lieu of depth migration. In those cases, the coefficient of the thin-lens term is negligibly small, which is why we often can get away with time migration in areas of mild to moderate lateral velocity variations.

The imaging problem is complicated when the overburden is as complex as that in Figure 8.0-13. Here, because of the bowties, a distorted diffraction traveltime trajectory indicates a false structure. The resulting complexity can give rise to more than one image ray. In this case, three image rays emerge at around midpoints 160, 250, and 370. Time migration fails here and imaging of this scatterer can only be achieved by depth migration.

We studied some examples of lateral velocity variations in Figures 8.0-8 through 8.0-13. The image ray behavior and the quality of focusing determine whether time or depth migration should be performed. If the starting and end points of the image ray have the same CMP location (Figure 8.0-9), only time migration is needed. A small amount of lateral deviation of the image ray (Figure 8.0-12) usually implies a well-focused time migration result and hence, a good representation of the geometric form of the subsurface. Large image-ray deviations imply grossly incorrect focusing, thus requiring depth migration rather than time migration (Figure 8.0-11). Finally, if more than one image ray is associated with a subsurface point (Figure 8.0-13), depth migration is imperative.

These observations on the point diffractor models are extended to a velocity-depth model that involves the reflecting interfaces in Figure 8.0-13. The image rays associated with this model are shown in Figure 8.0-14. There is no deviation from the vertical along the image rays down to horizon 2. Therefore, depth migration is not needed to image this horizon. On the other hand, the image rays significantly deviate from the vertical as they travel down to horizons 3 and 4. For example, the image ray starting at CMP location 140 reaches horizon 4 approximately beneath CMP location 180; a lateral shift of 40 midpoints. Proper imaging of these two horizons is achievable only by depth migration.

From Figure 8.0-11, note that time migration collapses the energy to apex A of the diffraction curve that coincides with the image ray location at the surface. In principle, the time migration output then can be converted to depth along the image rays, rather than along the vertical rays [3]. Mapping along the image rays performs some of the action that is associated with the thin-lens term. Remember that at each downward-continuation step, the action of the thin-lens term is equivalent to a time shift that depends on spatial velocity variation. Since the thin-lens and diffraction terms are applied in an alternate manner as the wavefield is downward continued in depth, the effects of these two terms are strongly coupled when the lateral velocity variation is as severe as shown in Figure 8.0-13.

When lateral velocity variation is moderate to strong, these two terms can, in principle, be separated and applied consecutively without significant error [4]. Full separation implies that a correction for the effects of the thin-lens term can be done either before or after time migration. If the correction is done after time migration, image-ray mapping should be used. If the correction is done before time migration, mapping using vertical time shifts usually is applied. In practice, a correction before time migration often performs better, since it tends to provide a better focused migration result. Nevertheless, all of the above conjectures are related only to now outdated migration algorithms since contemporary implementations of depth migration algorithms are based on splitting, and not separation, of the diffraction and thin-lens terms (Sections D.3 and D.4). It is not a common practice to correct for the effect of the thin-lens term by time-to-depth conversion of time-migrated data using image rays. However, it is a common practice to do time-to-depth conversion of time horizons interpreted from time-migrated data by using image rays (model building).

## 8.1 Layer replacement

In this section, we demonstrate that the problem of a complex overburden that involves only one layer boundary, such as an irregular water bottom, at which significant ray bending takes place can sometimes be addressed by prestack layer replacement followed by NMO correction, CMP stacking, and poststack time migration. In both cases, the approaches are based on the philosophy of revising velocity estimates and obtaining an improved unmigrated stacked section.

Consider the velocity-depth model in Figure 8.1-1a. Note that complex geometry of the boundary (horizon 2) between the overburden and the substratum and the significant velocity contrast across this boundary cause the severe ray bending. This in turn causes distortions and disruptions of the underlying target reflections. Without the velocity contrast (Figure 8.1-1b), the rays would not bend and there would be no need for depth migration. Figure 8.1-1 suggests that replacing the overburden velocity with the substratum velocity can be a viable alternative to using depth migration to remove the deleterious effects of a complex overburden, such as an irregular water-bottom topography, on the substrata. This technique is known as layer replacement.

A technique for layer replacement [5] based on wave-equation datuming [6], [7] is presented here. Berryhill’s datuming technique involves extrapolating a known wavefield at a specified datum of arbitrary shape to another datum, also of arbitrary shape. Wave extrapolation is performed using the Kirchhoff integral solution to the scalar wave equation. It incorporates both the near-field and far-field terms (Section H.1). The velocity used in extrapolation is that of the medium confined between the input datum and the output datum.

Figure 8.1-2  Layer replacement by Kirchhoff summation. (a) Input zero-offset section (at datum level z = 0) based on a velocity-depth model consisting of three point scatterers buried at depths of 800, 1300, and 1900 m. Interval velocities are indicated on the right side of the section. (b) Step 1: Downward continuation of the wavefield from z = 0 to datum level z = 800 m using a velocity of 2000 m/s. (c) Step 2: Upward continuation of (b) back to z = 0 using a velocity of 2500 m/s. (d) The zero-offset section derived independently using interval velocities that are indicated on the left side of the section. This section should be compared to (c).

### Wave-equation datuming

Figure 8.1-2 shows a simple case of datuming. A zero-offset section is computed over three point scatterers buried beneath midpoint location A in a medium with the layered velocity structure as denoted between Figures 8.1-2a and b. The point diffractors are situated at the layer interfaces at 800, 1300, and 1900 m depths. The traveltime trajectory associated with the shallow point diffractor is a hyperbola. The traveltimes associated with the deeper diffractors are nearly hyperbolic (normal moveout).

Extrapolate the zero-offset wavefield at z = 0 (Figure 8.1-2a) using the velocity of the first layer (2000 m/s) and compute the wavefield at the first interface, z = 800 m. The result is shown in Figure 8.1-2b. As expected, the hyperbola associated with the shallow point scatterer largely collapses to its apex since this scatterer is located at the first interface. Because the receivers now are closer to the other two deeper scatterers, events associated with them also are compressed. Figure 8.1-2b shows the zero-offset section that would have been recorded if the receivers were placed along the first interface. The energy from the shallow point scatterer on this section now arrives at t = 0 because the datum for this section is the interface at which the scatterer is located.

Now, extrapolate the wavefield at the first interface (z = 800 m) shown in Figure 8.1-2b back up to the surface (z = 0) using the velocity of the second layer (2500 m/s). Figure 8.1-2c shows the result, which is compared to the zero-offset section in Figure 8.1-2d. The latter was derived independently from the velocity-depth model denoted between Figures 8.1-2c and d. With this two-step wave extrapolation, the first layer with the 2000-m/s velocity was replaced with the second layer with the 2500-m/s velocity.

Figure 8.1-3  Poststack layer replacement involves two steps. Step 1: The zero-offset section (a) is extrapolated down to the water bottom (horizon 2 in Figure 8.1-1a) using the water velocity. Section (b) is obtained when the receivers are placed along the water bottom. Step 2: This intermediate wavefield is extrapolated back up to the surface using the velocity of the stratum below the water bottom (2000 m/s). The resulting zero-offset section (c) can be compared against the zero-offset section derived independently (d) using the same velocity-depth model as in Figure 8.1-1a, except the overburden velocity is the same as that of the substratum (2000 m/s) as shown in Figure 8.1-1b.

In principle, wave-equation datuming can be performed by any extrapolation method based on phase-shift, finite-difference, or Kirchhoff summation techniques. Nevertheless, the Kirchhoff summation is more convenient in handling datum surfaces with arbitrary shapes.

It is important to distinguish between datuming and migration. Datuming produces an unmigrated time section at a specified datum z(x), which can be arbitrary in shape. Migration involves computing the wave-field at all depths from the wavefield at the surface. In this respect, datuming is an ingredient of migration, when migration is done as a downward-continuation process. In addition to downward continuation, migration, of course, requires invoking the imaging principle (t = 0) (migration principles).

Wave-equation datuming has several practical applications — horizon flattening, forward modeling of seismic wavefields, and layer replacement. These are performed in either the prestack or poststack mode. The main difference between the two implementations is that the velocity must be halved when doing poststack datuming to conform to the exploding reflectors model (introduction to migration). Horizon flattening involves downward continuing from one marker horizon to another. When this is done successively for all marker horizons in a section, the technique is useful in reconstructing the past structural history in a given survey area [8]. Seismic modeling involves a series of upward continuations through a specified set of velocity interfaces starting at the bottom of the model and ending at the top. An example is provided in Section K.1.

### Poststack layer replacement

Now consider another practical application of wave-equation datuming to 2-D surface seismic data — removing the degrading effect of an irregular water-bottom topography on the continuity and geometry of reflections below. This problem is particularly severe in areas with a strong velocity contrast between the water layer and the substratum. Despite the usual 3-D nature of the problem, the 2-D interpretation of the target reflections often can be improved by replacing the velocity of the water layer with the velocity of the substratum.

Figure 8.1-4  (a) A modeled common-shot gather from the complicated part of the velocity-depth model in Figure 8.1-1a before (a) and after (b) layer replacement. Limitations in modeling with ray tracing cause abrupt terminations along the moveout curves and amplitude glitches in (a). In turn, these caused the spurious diffractions in (b) during the upward continuation steps.

Consider the zero-offset section shown in Figure 8.1-3a based on the velocity-depth model shown in Figure 8.1-1a. Poststack layer replacement involves two extrapolation steps:

1. The first step in poststack layer replacement involves downward continuing the wavefield at the surface (Figure 8.1-3a) to the water bottom (horizon 2 in Figure 8.1-1a) using the water velocity in extrapolation. The intermediate result is shown in Figure 8.1-3b. Note that in this horizon-flattened section, the water-bottom reflection is at t = 0, which means that all receivers are situated on the irregular water bottom. If we specified the overburden velocity or the water-bottom topography incorrectly, then the water-bottom reflection would not be at t = 0. In this respect, the intermediate section becomes a useful diagnostic tool before moving to the next step. In fact, wave-equation datuming actually can be applied layer by layer for structural model restoration. At each layer boundary, by examining the flatness of the event at t = 0 and observing any arrival-time departures of the event from t = 0, the validity of an estimated velocity-depth model can be verified.
2. The second step in poststack layer replacement involves upward continuation of the intermediate wavefield (Figure 8.1-3b) back to the surface z = 0 using the velocity of the substratum (2000 m/s). Figure 8.1-3c is the zero-offset section at z = 0 after layer replacement.

A zero-offset section was created from the same velocity-depth model (Figure 8.1-1a) as for Figure 8.1-3a, except that the first layer velocity was set to 2000 m/s (Figure 8.1-1b). Compare this zero-offset section as shown in Figure 8.1-3d with the output of layer replacement as shown in Figure 8.1-3c, and note that the two sections are largely equivalent. Both layer replacement and depth migration are processes aimed at removing the effects of the complex overburden. However, note that layer replacement only requires accurate representation of the overburden (horizon 2 in Figure 8.1-1a), while depth migration requires accurate representation of the entire velocity-depth model (Figure 8.1-1a). Also note that the output from depth migration is a migrated depth section, while the output from layer replacement is an unmigrated time section (Figure 8.1-3c). After eliminating the complex overburden effect, this section only requires time migration.

### Prestack layer replacement

As with depth migration, layer replacement after stack can remove the effect of a complex overburden, provided the input section accurately represents a zero-offset section. However, the complex overburden causes raypath distortions that generate anomalous, nonhy-perbolic moveout patterns in prestack data. Poststack depth migration does not produce a completely accurate image of the subsurface, even when the velocity-depth model is known accurately. Similarly, poststack layer replacement does not remove the effect of complex overburden entirely, even if its geometry is known accurately, because the input stacked section differs from the zero-offset section. Nevertheless, based on simple ray tracing we can determine whether poststack layer replacement can delineate the underlying structure. If the effort yields unsatisfactory results, layer replacement should be performed before stack.

Complex moveout is evident in the modeled common-shot and CMP gathers in Figures 8.1-4a and 8.1-5a, respectively. These gathers were modeled from the velocity-depth model in Figure 8.1-1a by using ray tracing that neither includes diffractions nor properly modeled amplitudes. The offset range is 50 to 2387.5 m with 12.5-m receiver spacing. A total of 437 shot gathers was generated, each with 192 traces. The coverage is uniform along the line and is 96-fold.

Figure 8.1-5  (a) A modeled CMP gather and (b) velocity spectrum from the complicated part of the velocity-depth model in Figure 8.1-1a before layer replacement, (c) the same gather as in (a) and (d) velocity spectrum after layer replacement.

Starting with the common-shot gathers, prestack layer replacement involves the following steps:

1. Downward continue all receivers to the output datum using the overburden velocity.
2. Sort the data to common receiver gathers.
3. Downward continue all shots to the same output datum using the overburden velocity.
4. Upward continue all shots back to the surface using the velocity of the substratum.
5. Sort the data back to common-shot gathers.
6. Upward continue all receivers back to the surface using the substratum velocity.

This series of operations eliminates the traveltime distortions associated with the water bottom, as seen in the common-shot and CMP gathers (Figures 8.1-4b and 8.1-5c, respectively). Despite the undesirable effects caused by the limitations of modeling with ray tracing, the complexity of the reflections associated with the two events (horizons 3 and 4 in Figure 8.1-1a) was reduced by layer replacement. Once the complex overburden effect is removed, these events have the hyperbolic moveout as seen in Figure 8.1-4b. Events A, B, and C are associated with horizons 2, 3, and 4, respectively, as shown in Figure 8.1-1a. Note that horizon 3 is flat. Therefore, the apex of the hyberbola is at the near-offset trace (event B). On the other hand, horizon 4 dips down to the right (Figure 8.1-1a). Therefore, the apex of the hyberbola shifts up-dip (as indicated by the arrow in Figure 8.1-4b).

Compare the velocity spectra in Figures 8.1-5b and d, and note the improvement in the velocity estimates after layer replacement for reflections beneath the complex overburden. The velocity analysis before layer replacement yields a good pick for the water-bottom reflection A. However, picks B and C, which are associated with the deeper layers (horizons 3 and 4 in Figure 8.1-1a), are not distinct. Similarly, events B and C also are indistinct on the maximum correlation curve plotted on the right side of the velocity spectrum (Figure 8.1-5b). After layer replacement (Figure 8.1-5d), note that the picks (denoted by ×) associated with the three events are distinct in the velocity spectrum and the correlation curve. The dipping water-bottom reflection A now has considerably higher moveout velocity, 2300 m/s (Figure 8.1-5d), as compared to the original velocity 1600 m/s (Figure 8.1-5b). The velocity pick for the flat event B from Figure 8.1-5d is 2000 m/s, which is the velocity of the medium above this reflector after layer replacement.

Just as it is true for any wave extrapolation-based process, the wave-theoretical layer replacement technique described here suffers from spatial aliasing. With modern data acquisition, we can record with many channels at small shot and group intervals. Therefore, spatial aliasing should not be an issue with more recent data. With this provision, prestack layer replacement followed by poststack time migration is a practical alternative to depth migration before stack in areas where simple geology is overlain by a single-layer overburden with a strong lateral velocity variation, such as an irregular water bottom.

The mathematical details of wavefield extrapolation based on the Kirchhoff integral are given in Section H.1. For poststack layer replacement, we assume that the stack is a zero-offset section. A zero-offset section is equivalent to the exploding reflectors model, provided the medium velocity is halved (introduction to migration). We assume that sources already are situated along the reflectors in the subsurface. Thus, we only need to move the receivers from one datum to another during downward and upward continuation steps. For prestack layer replacement, each common-shot and common-receiver gather is extrapolated, independently. In particular, the wavefield at a point on the output datum is computed using all the traces in the input gather. The output gather should be computed beyond the lateral extent of the input gather to prevent possible loss of steeply dipping events. For prestack layer replacement, the velocity used in the extrapolation is that of the medium between the input datum and the output datum.

### Field data example

Now consider the field data example in Figure 8.0-3a. Layer replacement will be done on these data before and after stack to remove the effect of the irregular water bottom. First, we must define the geometry of the overburden, in this case, the water-bottom topography. Since the overburden velocity is constant (1475 m/s), we migrate the CMP-stacked section (Figure 8.0-3a) using the constant-velocity Stolt method, digitize the water bottom, and convert it to depth as shown in Figure 8.1-6a.

First, consider poststack layer replacement. We downward continue the CMP stack (Figure 8.0-3a) by using the water velocity from the surface to the water bottom defined in Figure 8.1-6a to get the horizon-flattened section in Figure 8.1-6b. As usual, we assume that the CMP stack is a zero-offset wavefield. The water-bottom reflection is approximately flat and is situated at t = 0, which indicates that the velocity-depth model for the water layer in Figure 8.1-6a is fairly accurate. Although incorrect in shape, substratum reflections are quite continuous on this section, indicating that we achieve a focusing effect by downward extrapolation. The bottom of the section reflects the mirror image of the water-bottom topography.

The next step is to take this wavefield back to the surface, this time using the velocity of the substratum (2500 m/s). Luckily, the velocity derived from the seismic data is fairly constant across the section for the substratum. The section after upward continuation is the result of layer replacement and is shown in Figure 8.1-6c. After eliminating the effect of the complex overburden, only time migration is needed to image this section (Figure 8.1-6d).

The results of prestack layer replacement now are examined. Starting with the common-shot gathers and using the velocity-depth model in Figure 8.1-6a, we perform downward and upward continuations of common-shot and common-receiver gathers. The intermediate steps for selected common-shot gathers are shown sequentially in Figure 8.1-7. Note the arrival time of the water-bottom reflection at t = 0 on the gathers in step 2 after downward continuing both the shots and receivers to the water bottom.

A selected set of CMP gathers before and after layer replacement, sorted from the shot gathers in steps 0 and 4, respectively, is shown in Figure 8.1-8. Because the data contain strong diffracted multiples, it is difficult to evaluate the layer replacement results. As in any other one-way extrapolation technique, wave-equation datuming treats multiples as primaries. From the velocity analysis in Figure 8.1-9, we are encouraged by the picking we can do on the velocity spectrum after layer replacement. However, the stack has the final say and is shown in Figure 8.1-10a. Again, once time-migrated (Figure 8.1-10b), the result of layer replacement should provide an accurate image of the substratum, free from the effects of complex overburden (compare Figure 8.1-10b with Figures 8.0-3d and 8.1-6d).

The interpretation of the unconformity T from the results of prestack layer replacement (Figure 8.1-10b) closely agrees with the proposed velocity-depth model in Figure 8.0-3c. The unconformity is continuous beneath midpoint C, where there is a structural high. Below and to the left of midpoint A, the unconformity extends down to the right with some tensional faults into the ancient continental shelf below midpoint B, then finally extends to the continental slope of that age to the right of the structural high below midpoint C. The present-day situation is just the opposite, with the continental shelf on the right. Note that shallow reflector R now is conformable with the rest of the shallow section above the unconformity.

For the irregular water-bottom case, it was shown that prestack layer replacement followed by NMO correction, stack, and time migration after stack with depth conversion basically is equivalent to depth migration before stack. The situation may be more complicated when dealing with a complex overburden problem that involves more than one interface. In principle, datuming can be applied layer by layer and the effect of the overburden can be removed. Nevertheless, there is no practical alternative to depth migration under those circumstances.

## 8.2 2-D poststack depth migration

Figure 8.2-1 shows a velocity-depth model for a salt pillow. The aspect ratio of the horizontal and vertical axes is 1; hence, the diagram exhibits the true shape of the diapiric structure. The model can be treated in three parts — the constant-velocity overburden above the salt, the salt diapir itself, and the substratum that includes the flat reflector below. So far as the flat reflector is concerned, the salt diapir constitutes a complex overburden structure with strong lateral velocity variations. Note the significant velocity contrast across the top-salt boundary and the undulating reflector geometry of the base-salt boundary — both give rise to ray bending that can only be handled by imaging in depth.

A total of 154 shot records were modeled along the lateral extent of the velocity-depth model in Figure 8.2-1 using the two-way acoustic wave equation. Shot and receiver group intervals both are 50 m, and the trace spacings of the CMP-stacked and zero-offset sections are 25 and 50 m, respectively. In Figure 8.2-1, the velocity-depth model, the CMP-stacked and zero-offset sections have been displayed with the correct lateral position with respect to one another. Note the focusing and de-focusing of the reflection amplitudes associated with the base-salt boundary on the zero-offset and CMP-stacked sections.

Selected shot records are shown in Figure 8.2-2. Each shot record consists of 97 channels associated with a split-spread geometry. The shot interval is 50 m and the receiver group interval is 50 m. The maximum offset is 2350 m. Note the complexity of reflection times on shot records located over the salt diapir and the variations in reflection amplitudes caused by the focusing and defocusing of rays.

The complexity of reflection traveltimes surprisingly is not as apparent on the CMP gathers (Figure 8.2-3) as it is on the shot records (Figure 8.2-2), although variations in reflection amplitudes can be observed with ease. CMP stacking, however, clearly shows the effect of nonhyperbolic moveout (Figure 8.2-1). Compare with the zero-offset section and note the traveltime and amplitude distortions on the CMP-stacked section. Departure of the stacked section from the true zero-offset section imposes a limitation on the accuracy of the image we get from poststack depth migration.

We now closely examine the degree of ray bending at layer boundaries in the salt diapir model. Shown in Figure 8.2-4 are the normal-incidence rays from each of the three layer boundaries — top-salt, base-salt and the flat reflector below, and the computed zero-offset traveltimes. Superposition of the modeled traveltime trajectories yields the zero-offset traveltime section associated with the velocity-depth model as shown in Figure 8.2-4. The zero-offset section in Figure 8.2-1 represents a modeled zero-offset wavefield, whereas the zero-offset section in Figure 8.2-4 represents a modeled zero-offset traveltime profile.

### Image rays and lateral velocity variations

Normal-incidence rays are associated with zero-offset traveltimes and therefore can be used to examine the degree of complexity in velocity-depth models as demonstrated in Figure 8.2-4. For a quantitative assessment of lateral velocity variations, however, image rays need to be examined as shown in Figure 8.2-5. By definition, image rays emerge at the right angle to the surface. As shown in Figure 8.0-11, the lateral shift between the point of departure of the image ray at the reflector position and the point of emergence of the image ray at the surface provides a measure of lateral velocity variation.

Consider the image rays departing from the top-salt layer boundary in Figure 8.2-5. These rays show no lateral shift, and therefore, imaging the top-salt boundary does not require depth migration; instead, it can be achieved by time migration. The image rays from the base-salt boundary, however, show significant lateral shifts, especially beneath the flanks of the diapir. The stronger the lateral velocity variations, the more the lateral shifts in image rays. This behavior of the image rays indicate that the lateral velocity variations caused by the salt diapir require depth migration to image the base-salt boundary, accurately.

The image rays associated with the flat reflector below the salt diapir also show significant lateral shifts (Figure 8.2-5). Again, this reflector can only be imaged accurately by depth migration, rather than time migration. Note that image rays do not sample the reflector boundaries uniformly — there are regions that contain densely and sparsely populated image rays.

In principal, an earth image in depth can be obtained by first migrating a stacked section in time, then converting the time-migrated section to depth along image rays using the appropriate velocity-depth model [3] [4]. This ray-theoretical two-step depth migration to obtain an earth image in depth is rarely used in practice. However, it is common practice to perform time-to-depth conversion of time horizons using image rays. Specifically, 3-D volume of stacked data first is migrated in time and selected time horizons are interpreted. These time horizons are then converted to depth horizons along image rays, again, using an appropriate velocity-depth model. Creating depth structure maps using this procedure is called map migration.

In conclusion, by examining the behavior of image rays through the salt diapir model, we can judge as to which layer boundary requires imaging in depth (Figure 8.2-5). The image rays down to the top-salt boundary are not distorted laterally; therefore, time migration is adequate for imaging the overburden above the salt diapir. Significant ray bending, however, takes place at the top-salt boundary; this results in lateral distortions of the image rays down to the base-salt boundary and the deeper reflector. Depth migration, therefore, is needed for accurate imaging of the base-salt boundary and the subsalt region.

### Time versus depth migration

We continue our discussion with the salt diapir model data. Figure 8.2-6 shows time and depth migrations of the zero-offset section associated with the salt diapir model shown in Figure 8.2-1. Ignore the dispersive noise in the vicinity of the top-salt event on the migrated sections; this is caused by the differencing approximations made to the differential operators in the implicit scheme used here. While the top-salt boundary is imaged accurately by both time and depth migration, note the overmigration exhibited by the base-salt and the deeper event in the time-migrated section. Depth migration, on the other hand, images these two reflectors, accurately, if the velocity-depth model input to depth migration is the true model as in the case of Figure 8.2-6. Although we used the true velocity-depth model, time migration produced the incorrect image of the base-salt boundary and the subsalt region. Time migration algorithms do not include the term that accounts for strong lateral velocity variations as manifested by the image rays shown in Figure 8.2-5. On the other hand, depth migration algorithms include this term and thus are able to correct for the lateral shift in image rays.

Figure 8.2-7 shows time and depth migrations of the CMP-stacked section associated with the salt diapir model shown in Figure 8.2-1. Again, the top-salt boundary is imaged accurately by both time and depth migration. As expected, however, time migration fails to produce a correct image of the base-salt boundary and the deeper reflector. Note that even depth migration fails to image these two reflectors with sufficient accuracy, although the velocity-depth model input to depth migration was the true model. Note, for instance, the subtle distortions on the base-salt event and the not-so-flat reflector in the depth-migrated section. This is a direct consequence of the fact that the CMP-stacked section is only a close representation of the zero-offset wavefield in the presence of strong lateral velocity variations associated with complex overburden structures. Since poststack migration algorithms are based on the zero-offset wavefield theory (migration principles), application of zero-offset migration to a CMP-stacked section would produce less-than-ideal results. To circumvent this deficiency in CMP stacking, and to correctly image the substratum that includes the base-salt boundary and the flat reflector, strictly, one needs to do prestack depth migration. Prestack depth migration is reviewed in the next section; however, for comparison, results of the zero-offset, poststack and prestack depth migrations are shown in Figure 8.2-8. Note that prestack depth migration produces an image that is free of the distortions observed on the image produced by poststack depth migration. Imaging accuracy is similar to that of the zero-offset section.

A way to minimize departure of a stacked section from a zero-offset section — that is, to minimize traveltime and amplitude distortions caused by nonhyper-bolic moveout during CMP stacking, is to use partial stacking. Figure 8.2-9 shows a portion of a full-fold CMP-stacked section that contains a salt diapir structure. (The CMP fold is 30.) Note the spurious horst-like structure at the base-salt boundary B. This data set is from the Red Sea where such structures are common. Hence, at first, this horst structure appears to be geologically plausable. Note, however, the horst block at the base-salt is replaced with a continuous event on the single-fold near-offset section. What may seem to be a horst block actually is no more than a manifestation of amplitude and traveltime distortions caused by the full-fold stacking that spans the entire offset range of the recorded data. The conclusion we can draw from this observation is that it is not necessarily the full-fold stack, but the near-offset section that better resembles a zero-offset section. Hence, it is the near-offset section, and not the full-fold section, that is an appropriate input to poststack depth migration. An immediate objection to this conclusion is that the near-offset section contains a significant amount of multiple energy that was naturally attenuated during CMP stacking by way of velocity discrimination between primaries and multiples. Also, note the loss of continuity of the base-salt event on the near-offset section; in contrast, CMP stacking improves the signal-to-noise ratio. A way to benefit from the power of CMP stacking to attenuate multiples and improve the signal-to-noise ratio while circumventing the nonhyperbolic moveout effect is to do partial stacking. By a simple series of tests, one can judge as to what portion of the cable — near offsets, mid-range offsets or far offsets, provides this optimum stack as input to poststack depth migration. It is important to note that, while poststack depth migration may require a stack based on a subset of offsets, prestack depth migration requires all offsets.

We must remind ourselves that an accurate image from depth migration is attainable only when the velocity-depth model is correctly defined, independent of the input data type — zero-offset, stack, or prestack. Figure 8.2-10 shows results of the zero-offset, post-and prestack depth migrations of the salt diapir model data using an incorrect velocity-depth model. An incorrect velocity-depth model causes a poor image produced by not just poststack depth migration, but also by zero-offset and prestack depth migration.

### Iterative depth migration

Historically, depth migration has been used in an iterative manner to obtain an earth image in depth from CMP-stacked data. When performed iteratively, depth migration is done using an initial velocity-depth model and the result is interpreted for the layer boundaries included in the model. The velocity-depth model then is modified accordingly and depth migration is performed once more. The process is continued until convergence is achieved.

Convergence means that what is input to depth migration as the velocity-depth model matches with the velocity-depth model inferred from the output from depth migration. With iterative depth migration, we know that we have achieved convergence when differences between the velocity-depth models from two consecutive iterations are minimal. We shall demonstrate that, by way of convergence, the final velocity-depth model from iterative depth migration, albeit not guaranteed to be accurate, can be made at least consistent with the input data. Consistency means that the modeled zero-offset traveltimes match with the observed reflection traveltimes on the stacked data associated with the layer boundaries included in the velocity-depth model. Convergence and consistency are the two necessary, but not sufficient, conditions for an earth model to be certified as a valid, geologically plausable solution from seismic inversion. For a velocity-depth model to be valid, a further requirement is that it also needs to be consistent with prestack data; thus, the strategic requirement for doing prestack depth migration.

In this section, we shall examine iterative depth migration for three data types — zero-offset, CMP stack, and prestack. We shall consider six different initial velocity-depth models, each having errors in layer velocities and/or reflector geometries. Convergence rates and the end results for each starting model will be evaluated and some practical conclusions about iterative depth migration will be drawn from these experiments.

Figure 8.2-11  Iterative depth migration using the zero-offset section as in Figure 8.2-1 and Model A (same as the true velocity-depth model as in Figure 8.2-1).

### Iteration with zero-offset data

Depth migration of the zero-offset section in Figure 8.2-1 using Model A shown in Figure 8.2-11 produces an accurate image of the salt diapir and yields the same output model as the input model in a single iteration. This happened because the input data set is the true zero-offset section and the input velocity-depth model (Model A in Figure 8.2-11) is the same as the true velocity-depth model (Figure 8.2-1).

Now consider Model B in Figure 8.2-12 in which the salt velocity and the base-salt boundary are specified incorrectly. There are practical reasons behind these deliberate errors introduced into the model. In practice, often the top-salt boundary is determined with reasonable accuracy by way of time migration. However, accurate delineation of the base-salt boundary is almost impossible with time migration. Additionally, the salt velocity may vary because of dolomite or shale intrusions, and these variations are often difficult to determine.

Start with Model B as the initial velocity-depth model and perform depth migration (Figure 8.2-12). The result indicates that the geometry of the base-salt boundary and the flat reflector below has changed (top of the right column). Interpret all the layer boundaries — top-salt, base-salt, and the deeper reflector, from the result of depth migration and create an updated velocity-depth model (center of the left column). Use this new model and perform depth migration once more (center of the right column). Interpret the new image from depth migration, again, for all the three layer boundaries and obtain an updated velocity-depth model (bottom of the left column). Finally, perform depth migration for the third time (bottom of the right column) and interpret for the three layer boundaries, once more. After this third iteration, we find that the velocity-depth model does not change from the previous iteration; hence, convergence is achieved. Nevertheless, the final solution has converged to a velocity-depth model that is different from the true model. We conclude from this experiment that convergence and consistency, albeit required for a model to be certified as a valid solution from depth migration, do not guarantee that the solution is the true model. The result of iterative depth migration obviously is dictated by the parameters of the initial velocity-depth model.

Consider now an initial model in which the geometry of only one reflector — that associated with the base-salt boundary, is in error while layer velocities are specified correctly (Figure 8.2-13). This may be possible in practice if there is an abundance of well data in the area of investigation that enables a reliable specification of layer velocities.

It often is preferable to specify a simple, in this case a flat reflector, geometry for the base-salt boundary. We shall let iterative depth migration produce a final geometry for the base-salt boundary. Now, start with Model C as the initial velocity-depth model and perform iterative depth migration (Figure 8.2-13). Convergence is achieved after three iterations. The final solution has converged to a velocity-depth model that is fairly close to the true model. We conclude that true geometry of reflectors can be recovered by iterative depth migration provided layer velocities are correctly specified in the initial velocity-depth model.

The case of Model D shown in Figure 8.2-14 is similar to that of Model B (Figure 8.2-12). While the layer velocity for the salt diapir in Model B has been specified erroneously too low, it has been specified in Model D erroneously too high. As a result, the final velocity-depth model, once again, departs from the true model, significantly.

Sometimes time migration may yield an inaccurate image of the top-salt boundary because of the incorrect overburden velocity field. Such is the case in Model E (Figure 8.2-15) where not only the overburden velocity, and therefore, the top-salt boundary are incorrect, but also the salt velocity and the base-salt boundary are in error. This model has too many errors in layer velocities and reflector geometries. It has taken four iterations to achieve convergence. The final solution radically departs from the true model.

If the initial model contains significant errors but only in reflector geometries, and layer velocities are specified correctly and thus are not altered from one iteration to the next, then convergence to the true velocity-depth model can be achievable (Figure 8.2-16).

To recap the results of the analysis of the six models (Figures 8.2-11 through 8.2-16), refer to Figure 8.2-17. Starting with an initial velocity-depth model (top left), which contains some errors in model parameters, perform iterative depth migration. After n iterations, obtain an earth image in depth (top right) that corresponds to an earth model in depth (center left) that is different from the true model (center right). Nevertheless, forward modeling of zero-offset traveltimes (bottom left) using the estimated model (center left) from iterative depth migration yields results that are consistent with the traveltimes on the zero-offset wavefield section (bottom right) used as input to iterative depth migration. We have met the convergence and consistency criteria for the final solution from iterative depth migration, but we have only obtained a solution, and not the solution — the true model (Figure 8.2-11).

Figure 8.2-18 shows the initial velocity-depth models (Figures 8.2-11 through 8.2-16) that contain various types of errors in model parameters. Figure 8.2-19 shows the final solutions derived from iterative depth migration using the initial models in Figure 8.2-18. Note that those initial models with correctly specified layer velocities but incorrectly defined reflector geometries (Models A, C, and F) converge to the true model (Model A). Those initial models with incorrect layer velocities (Models B, D, and E), however, do not converge to the true model.

Zero-offset traveltime sections computed from the final velocity-depth models are shown in Figure 8.2-20. Compare these traveltime sections with the input zero-offset wavefield section (Figure 8.2-11), and note that they all are consistent with the latter. Consistency, however, does not guarantee that the final solution from iterative depth migration yields the true model. In practice, we will never know which one of the final solutions shown in Figure 8.2-19 corresponds to the true model.

### Iteration with CMP-stacked data

Depth migration of the CMP-stacked section in Figure 8.2-1 using Model A required two iterations to achieve convergence (Figure 8.2-21). The resulting image and the model inferred from it have some inaccuracies. It appears that starting with the true model (Model A), iterative poststack depth migration does not exactly reproduce the true model. For the zero-offset section, starting with the true model, convergence was achieved after one iteration and the resulting model was the same as the input model (Figure 8.2-11). Again, the reason for the less-than-ideal performance of poststack depth migration is that the CMP stack is only an approximation to the zero-offset wavefield. The more the stacked section departs from the zero-offset wavefield, the more the final model will depart from the true model.

Now consider Model D in Figure 8.2-22 in which the salt velocity and the base-salt boundary are incorrectly specified. After the third iteration with the stacked section, we find that convergence is achieved. Nevertheless, the final solution has converged to a velocity-depth model that is different from the true model.

If the initial model contains errors in reflector geometries only, and layer velocities are specified correctly, then convergence using the stacked section to the true velocity-depth model is nearly achievable (Figure 8.2-23).

### Iteration with prestack data

To compare with iterative depth migration applied to zero-offset and CMP-stacked data, we shall examine the performance of iterative depth migration applied to prestack data. Techniques for prestack depth migration itself, however, will be reviewed in the next section. Depth migration of the prestack data using the true model (Model A) shown in Figure 8.2-11 required only one iteration to achieve convergence (Figure 8.2-24).

Next, consider Model D in Figure 8.2-25 in which the salt velocity and the base-salt boundary are incorrectly specified. After the third iteration with the prestack data, we find that convergence is achieved. Nevertheless, the final solution has converged to a velocity-depth model that is different from the true model.

If the initial model contains errors in reflector geometries only, and layer velocities are specified correctly, then convergence using the prestack data to true velocity-depth model is achievable (Figure 8.2-26).

### Iteration in practice

In this section, we tested iterative depth migration of zero-offset, CMP-stack and prestack data using six different initial velocity-depth models. Aside from the true model, each model was corrupted by errors in layer velocities and/or reflector geometries. We summarize the results of the model experiments as follows:

1. When will the image from iterative depth migration converge to the true velocity-depth model? The answer to this question is many-fold, and it depends on the type of input data and errors in the initial velocity-depth model (Figure 8.2-27). If the input data set is zero-offset section and if the input velocity-depth model is the true model (an impossible case in practice), then the answer is yes. For prestack data, the answer also is yes, whereas for stacked data, the answer is nearly. If the initial velocity-depth model is in error of reflector geometries only, the answer is most likely, very likely, and likely for zero-offset, prestack and stacked data, respectively. Finally, if the initial velocity-depth model is in error of layer velocities, no matter what the input data set is, the answer is invariably no. Needless to say, iterative depth migration never guarantees that the solution converges to a true velocity-depth model.
2. It is wrong to discontinue an iterative application of depth migration without achieving convergence — the intermediate output does not represent a valid earth image in depth and neither does it infer a valid earth model in depth. Zero-offset traveltimes associated with this earth model are not consistent with the traveltimes in the input section.
3. In iterative depth migration, only one set of parameters — either layer velocities or reflector geometries, should be modified from one iteration to the next. Changing the parameter type at some intermediate iteration will only divert the solution to a different end result and will cause the convergence to that end result to take longer. It is advisable to keep layer velocities unaltered from one iteration to the next, and only modify reflector geometries by interpreting the output image from depth migration.
4. Number of iterations depends on how much the initial model departs from the true model. If the initial model contains errors in layer velocities, only, fewer iterations are required than if the model contains errors in reflector geometries. If the model contains errors both in layer velocities and reflector geometries, a large number of iterations is required to achieve convergence.
5. Iterative depth migration converges to an answer, which almost never corresponds to the true velocity-depth model. Therefore, the final velocity-depth model estimated from iterative depth migration needs to be calibrated to well data. Specifically, layer velocities are adjusted so as to match at well locations the depth values of the layer boundaries included in the final velocity-depth model with the well tops corresponding to those layer boundaries.

Figure 8.2-28 shows a portion of a CMP-stacked section and its time migration. The deepest layer with an abundance of diffractions represents salt with anhydrite-dolomite rafts. The objective is to obtain an accurate geometry for the base-salt boundary represented by the strong, deepest reflection.

To illustrate poststack iterative depth migration, start with a simple, horizontally layered earth model with constant layer velocities as the initial velocity-depth model (Figure 8.2-29a) and migrate the stacked section to obtain the depth image shown in Figure 8.2-28b. Superimpose the flat depth horizons associated with the initial velocity-depth model on this image section, and note that the reflector geometries implied by the depth image show discrepancy with the flat depth horizons. Discard the latter and interpret the image section to derive a set of structurally consistent depth horizons (Figure 8.2-28c). This completes the first iteration of poststack depth migration and model updating.

Next, keep the layer velocities the same as for the initial velocity-depth model (Figure 8.2-28a), but use the updated depth horizons (Figure 8.2-28c) to build a new velocity-depth model as shown in Figure 8.2-28d. Now, perform poststack depth migration once more and note that the reflector geometries implied by the depth image (Figure 8.2-28e) are still in discrepancy with the intermediate velocity-depth model (Figure 8.2-28d). Discard the depth horizons and reinterpret them from the image section (Figure 8.2-28f). This completes the second iteration of poststack depth migration and model updating.

Repeat the steps for model updating, poststack depth migration and re-interpreting the depth horizons for the third (Figures 8.2-28g,h,i), and fourth time (Figures 8.2-28j,k,l) until convergence is achieved; that is, the velocity-depth model input to depth migration is consistent with the reflector geometries inferred by the depth image. The final velocity-depth model and the depth image derived from it are shown in Figures 8.2-30a,b. Convergence criterion can be verified by normal-incidence modeling of the zero-offset traveltimes that correspond to the reflectors associated with the layer boundaries included in the velocity-depth model. Superimpose the modeled traveltimes on the unmigrated stacked section and note that the modeled and the actual traveltimes are in good agreement.

While the consistency of the modeled and actual zero-offset traveltimes verifies that the final velocity-depth model (Figure 8.2-30a) from poststack iterative depth migration meets the convergence criterion, the model is not guaranteed to be accurate. In fact, as demonstrated by the synthetic data examples in this section, there exists not just one but a multiple number of velocity-depth models that are consistent with the stacked data. An acceptable model is that which also is consistent with prestack data; thus, one strategic reason for prestack imaging is resolving the uncertainty in the acceptable velocity-depth models and reducing the many possible models to a few that are geologically plausable.

## 8.3 2-D prestack depth migration

In this section, we shall review two methods of prestack depth migration applied to 2-D data. These methods both have historical and conceptual significance. An early method of prestack depth migration is based on downward continuation of sources and receivers (Section D.1). The method commonly is known as shot-geophone migration. A common-shot gather represents a wavefield and thus can be extrapolated in depth at discrete intervals. As a result, receivers are lowered from one depth level to the next. By invoking reciprocity, we may also consider a common-receiver gather representing a wavefield. Thus, by extrapolating a common-receiver gather, sources are lowered from one depth level to the next. By alternating between extrapolation of common-shot and common-receiver gathers at each depth level, all sources and receivers are lowered from the surface to each of the reflectors in the subsurface. While sources and receivers are lowered vertically downward from one depth level to the next, the recorded waves are back-propagated along the raypaths from source to a reflector back to receiver locations at the surface. When sources and receivers are lowered to the reflector, they coincide and the traveltime diminishes to naught. This satisfies the imaging condition.

When the maximum depth of extrapolation is reached, traces at zero-offset from each of the resulting common-shot gathers are extracted and placed side-by-side to produce the image from prestack depth migration. Nonzero-offset traces are abandoned since all primary energy has collapsed to zero-offset provided the velocity-depth model is correct.

Note that in practical implementation of shot-geophone migration, data need to be sorted from one gather type to another (common-shot and common-receiver) at each depth level. This sorting operation increases the cost of the method, formidably. Another practical aspect of this method is that missing near-offset traces in common-shot gathers are filled in with zero traces before downward continuation is started.

A more popular method of prestack depth migration is based on migration of shot records, individually. This is plausable since a shot record is a wavefield generated by a single source. The method commonly is known as shot-profile migration. The migrated shot records are then sorted into common-receiver gathers. Finally, traces in each receiver gather are summed to construct the image below the receiver location. By placing the traces that result from this summation side by side, we obtain the image from shot-profile migration.

A practical advantage of shot-profile migration is its ability to handle irregularities in recording geometry. Since each shot record is migrated independently, missing shots, duplicate shots or just irregularities in shot spacing are irrelevant.

### Shot-geophone migration

Figure 8.3-1 shows selected common-shot gathers associated with a land line recorded over a salt diapiric structure. Note the presence of coherent reflections on shot records to the left and to the right of the salt diapir, and poor reflection quality of the shot records above the salt diapir. Following the sorting to common-receiver gathers, we observe the missing shots along the line (Figure 8.3-2).

Figure 8.3-3a shows the CMP-stacked section with the salt diapir in the middle. The base-salt reflection with the velocity pull-up typical of salt diapirs is the most coherent event in this portion of the section. The poor imaging of the base-salt boundary by poststack depth migration is a direct consequence of the fact that the CMP-stacked section is only an approximation to a zero-offset section (Figure 8.3-3b). Compare with the image in depth derived from shot-geophone migration shown in Figure 8.3-3c. The geometry of the base-salt boundary indicates the presence of fault blocks. Also note some coherent events, albeit weak, in the subsalt region.

A bonus effect of prestack migration — be it time or depth, is its ability to attenuate multiples. Note, for instance, in Figure 8.3-3a the pegleg (m) on the left flank of the diapir. This multiple has been retained during CMP stacking and moved with the primary velocity during poststack migration (Figure 8.3-3b). However, it is not on the section derived from prestack depth migration. To explain this bonus effect of prestack migration, refer to the selected shot records after prestack depth migration (Figure 8.3-4). Since the migration velocity field is associated with the primary events, the primary energy has focused on the zero-offset traces, while the multiple energy has been dispersed to nonzero-offset traces. The focusing of the primary energy also is seen on the common-receiver gathers after migration (Figure 8.3-5). The image from prestack depth migration (Figure 8.3-3) is formed by compiling the zero-offset traces from the shot records after migration (Figure 8.3-4). All the remaining nonzero-offset traces in the shot records are simply abandoned.

The quality of focusing at zero offset by shot-geophone migration understandably depends on the accuracy of the velocity field used in migration. Erroneously too low or too high velocities would cause a partial collapse of the primary energy to zero offset. Note the curvature along some of the events in the vicinity of the zero-offset traces on the shot records (Figure 8.3-4) and receiver gathers (Figure 8.3-5). Some exhibit a curvature upward (erroneously low velocity above the event) and some exhibit a curvature downward (erroneously high velocity above the event). Theoretically, by measuring this curvature as if it is associated with residual moveout, velocities may be updated at each shot or receiver location. Such residual moveout analysis has been developed for shot-profile migration [9] and is discussed later in this section.

We shall extend the discussion on shot-geophone migration by comparing prestack depth migration with prestack time migration results. Shown in Figure 8.3-6 are the images obtained from post- and prestack time migration. First, note that both yield an incorrect image of the base-salt boundary. Since the problem of imaging the base-salt boundary is associated with a complex overburden structure, it only can be handled by depth migration. Strictly speaking, it should be handled by prestack depth migration (Figure 8.3-3).

Second, prestack time migration actually can produce, as in this example, a poorer image of the base-salt boundary compared to poststack time migration (Figure 8.3-6). Doing the migration before stack does not necessarily solve the problem of imaging beneath the salt diapir; doing it as depth migration is the key to solving this problem. The poor performance of prestack time migration can be verified by examining the focusing of energy at zero offset on common-shot (Figure 8.3-7) and common-receiver gathers (Figure 8.3-8). Specifically, note that the focusing of energy to zero offset by prestack time migration is not as good as it is by prestack depth migration (Figures 8.3-4 and 8.3-5).

We now complete our discussion on shot-geophone migration by examining the effect of missing data on imaging quality. Figure 8.3-9 shows a model of a salt diapir. A total of 193 shot records was created by using a nonzero-offset ray-theoretical modeling procedure. The receiver cable is split-spread with an offset range of 501200 m, and 48 receivers at 50-m interval. Putting aside the inadequacies of ray theory for modeling traveltimes, we should be able to make use of the modeled shot records for investigating the missing-data problem.

Selected shot records before and after shot-geophone migration using the true velocity-depth model (Figure 8.3-9a) are shown in Figure 8.3-10. Note that the events have focused onto and at the vicinity of the zero-offset traces. By extracting the zero-offset traces and placing them side by side, we obtain the depth image in Figure 8.3-9b. (The amplitude weakening along the top-salt event is caused by the limitations in the ray-theoretical modeling.)

A total of 15% of the traces from the modeled shot records was discarded arbitrarily and replaced with zero traces. As a result, some shot records contained few zero traces while some contained all zero traces. Shown in Figure 8.3-11 are the same shot records as in Figure 8.3-10 with missing traces. Note that, after prestack depth migration, focusing of the energy to zero-offset and its vicinity (Figure 8.3-11) has been achieved in a manner comparable to the case of the complete data set (Figure 8.3-10). Likewise, the prestack depth-migrated section derived from the missing data (Figure 8.3-9c) is very similar to the section derived from the complete data set (Figure 8.3-9b). We may conclude that shot-geophone depth migration can accommodate missing data resulting from recording geometry irregularities.

### Shot-profile migration

We shall review a method of shot-profile migration [10] using the salt diapir model shown in Figure 8.2-1. A wave theoretical modeling scheme based on the two-way acoustic wave equation was used to generate 154 shot records. Selected shot records are shown in Figure 8.2-2. Shot spacing is 50 m and receiver spacing is 50 m. Each shot record contains 97 traces corresponding to a split-spread recording geometry with a maximum offset of 2350 m.

Figure 8.3-12a shows three shot records — one located away from the main diapiric body, one on the left flank, and one other on the right flank of the diapir. The velocity-depth model (Figure 8.3-12b), although vertically exaggerated, is the same as in Figure 8.2-1. The same shot records after shot-profile migration using the true velocity-depth model (Figure 8.3-12b) are shown in Figure 8.3-13a. Note that each shot record after migration represents partial image of the subsurface within a limited lateral extent.

Now, imagine all 154 shot records after migration placed at their corresponding shot locations along the line. Then, consider one specific receiver location. There will be traces from a number of migrated shot records that will coincide with this receiver location that are common to all. These traces consitute a common-receiver gather after shot-profile migration. Selected common-receiver gathers sorted from the migrated common-shot gathers are shown in Figure 8.3-13b. The number of traces in a common-receiver gather obviously is determined by the recording geometry.

If the velocity-depth model used in shot-profile migration coresponds to the true model, then each shot record should yield a correct image of the subsurface, albeit limited in lateral extent (Figure 8.3-13a). Then, traces from various shot records after migration at the same receiver location should represent the identical image below that receiver location. Stated differently, a common-receiver gather should contain flat events if the velocity-depth model used in shot-profile migration is correct. The final step of shot-profile migration involves the summation of the traces in each receiver gather to create the image in depth (Figure 8.3-13c). Note that the partial images with limited lateral extent from each individual shot record (Figure 8.3-13a) coincide with the complete image within the lateral extent of the line in its entirety.

Examine the common-receiver gathers in Figure 8.3-13b, closely. The receiver locations are labeled as 179, 209, and 249. Note that the events associated with the top-salt and base-salt boundaries, and the flat reflector below are positioned differently in relation to the receiver location of each gather. In case of a flat layer boundary, the event is positioned symmetrically with respect to the receiver location (179). In the case of a dipping layer boundary, such as the top-salt, the event is positioned to the right of the receiver location (209) or to the left of the receiver location (249), but always in the updip direction.

### Sensitivity of image accuracy to velocity errors

Accuracy of the velocity-depth model used in prestack depth migration can be checked by examining event curvature on common-receiver gathers. Figure 8.3-14a shows three common-receiver gathers sorted from the shot records that were migrated using the constant overburden velocity above the salt diapir (Figure 8.3-14c). Note that the top-salt event exhibits flat character at all three receiver locations since the velocity used for migration is the same as the layer velocity above the top-salt boundary (3000 m/s). Whereas the base-salt event and the event corresponding to the flat reflector below do not exhibit flat character since the migration velocity, in this case, is erroneously lower than the true layer velocities.

Use three different trial constant velocities for the overburden assigned to all the layers in the subsurface and migrate the shot records. Consider a common-receiver gather at one specific location (179 in Figure 8.3-14b). The top-salt event exhibits a curvature upward with erroneously low velocity (2500 m/s), flat character with the correct velocity (3000 m/s), and downward curvature with erroneously high velocity (3500 m/s). Moreover, the event depth changes from one trial velocity to another — the low velocity causes the event to appear at shallow depth.

Event curvature on common-receiver gathers can be likened to residual moveout on moveout-corrected common-midpoint gathers caused by incorrect move-out velocities. By measuring the residual moveout, layer velocities can be updated at each receiver location. Residual moveout analysis can be formulated within the context of common-receiver gathers derived from shot-geophone migration [9], common-receiver gathers derived from shot-profile migration [10] [11], or common-depth-point gathers (image gathers) derived from common-offset migration [12] [13].

We shall examine the sensitivity of image accuracy to velocity errors based on event curvature on common-receiver gathers derived from the shot-profile migration of the salt diapir model data (Figure 8.3-12). This analysis also is a prelude to using the event curvature on image gathers as a criterion for layer-by-layer velocity determination (next chapter). Start with the velocity-depth model defined as a half space with three different trial constant velocities that may be considered appropriate for the overburden above the salt diapir (Figure 8.3-15). Perform shot-profile migration using each of the velocity-depth models and sort the output to common-receiver gathers. Then, stack the receiver gathers to obtain the images in depth.

Figure 8.3-16 shows the depth image and a receiver gather (location 249) from each of the three prestack depth migrations. Note that the top-salt event exhibits a strong curvature upward with the trial velocity 2500 m/s for the overburden layer under consideration (Figure 8.3-16a). We conclude that this velocity is too low for the overburden. The top-salt event exhibits a strong curvature downward with the trial velocity 3500 m/s for the overburden layer (Figure 8.3-16e). We conclude that this velocity is too high for the overburden. Finally, the top-salt event is flat with the trial velocity 3000 m/s for the overburden layer (Figure 8.3-16c). We conclude that this velocity is appropriate for the overburden.

Examine the stack power in Figure 8.3-16 for the top-salt event on the depth images and note that the highest stack power is attained with the trial velocity 3000 m/s for the overburden layer (Figure 8.3-16d). In contrast, the stack power with the too-low (2500 m/s) and too-high (3500 m/s) velocities (Figures 8.3-16b and f) is weaker because of the curvature of the top-salt event on common-receiver gathers.

Based on the results of prestack depth migration using the three trial velocities (Figure 8.3-16), we assign the optimum velocity 3000 m/s to the overburden. Next, we interpret the top-salt event on the corresponding depth image (Figure 8.3-16d). Then, we modify the velocity-depth model in Figure 8.3-15b as shown in Figure 8.3-17, accordingly. Here, we have specified a two-layer model — the overburden above the top-salt boundary interpreted from Figure 8.3-16d, and the new half-space below with three different trial velocities appropriate for the diapir layer.

Figure 8.3-18 shows the depth image and a receiver gather (location 249) from each of the three prestack depth migrations. Note that the top-salt event is flat in all three cases since the overburden layer is the same for all the three models. The base-salt event exhibits a moderate curvature upward with the trial velocity 4500 m/s for the diapir layer under consideration (Figure 8.3-18a). We conclude that this velocity is too low for the diapir layer. The base-salt event exhibits a moderate curvature downward with the trial velocity 5500 m/s for the diapir layer (Figure 8.3-18e). We conclude that this velocity is too high for the diapir layer. Finally, the base-salt event is flat with the trial velocity 5000 m/s for the diapir layer (Figure 8.3-18c). We conclude that this velocity is appropriate for the diapir layer.

Examine the stack power in Figure 8.3-18 for the top-salt event on the depth images and note that it is identical in all three cases since the overburden layer is the same for all the three models. The highest stack power for the base-salt event is attained with the trial velocity 5000 m/s for the diapir layer under consideration (Figure 8.3-18d). In contrast, the stack power with the too-low (4500 m/s) and too-high (5500 m/s) velocities (Figures 8.3-18b and f) is slightly weaker because of the curvature of the base-salt event on common-receiver gathers.

Based on the results of prestack depth migration using the three trial velocities (Figure 8.3-18), we assign the optimum velocity 5000 m/s to the diapir layer. Next, we interpret the base-salt event on the corresponding depth image (Figure 8.3-18d). Then, we modify the velocity-depth model in Figure 8.3-17b as shown in Figure 8.3-19, accordingly. Here, we have specified a three-layer model — the overburden above the top-salt boundary interpreted from Figure 8.3-16d, the diapir with the base-salt boundary interpreted from Figure 8.3-18d, and the new half-space below with three different trial velocities appropriate for the subsalt region.

Figure 8.3-20 shows the depth image and a receiver gather (location 249) from each of the three prestack depth migrations. Note that the top-salt and base-salt events are flat in all three cases since the overburden and the diapir layers are the same for all the three models. The deepest event associated with the flat reflector within the half-space exhibits a very mild curvature upward with the trial velocity 3500 m/s for the half-space (Figure 8.3-20a). We conclude that this velocity is low for the half-space. The same event exhibits a very mild curvature downward with the trial velocity 4500 m/s for the half-space (Figure 8.3-20e). We conclude that this velocity is high for the half-space. Finally, the same event is flat with the trial velocity 4000 m/s for the half-space (Figure 8.3-20c). We conclude that this velocity is appropriate for the half-space that represents the subsalt region with the flat reflector.

Examine the stack power in Figure 8.3-20 for the top-salt and base-salt events on the depth images, and note that they are identical in all three cases since the overburden and diapir layers are the same for all three models. It is difficult to distinguish the images of the deepest event associated with the flat reflector within the half-space based on stack power (Figures 8.3-20b, d, and f). Nevertheless, based on the subtle differences in stack power in favor of the trial velocity 4000 m/s and the implausable undulations in the geometry of the flat reflector that resulted from the trial velocities 3500 m/s and 4500 m/s, we may conclude that the optimum velocity for the subsalt region is 4000 m/s.

This concludes the estimation of the velocity-depth model using the event curvature on common-receiver gathers. The conceptual appeal of this analysis can be very attractive for construction of earth models in depth. In practice, however, the computational cost can be prohibitive. Instead, residual moveout analysis of image gathers derived from common-offset migration is used for updating layer velocities (earth modeling in depth).

Refer to the common-receiver gathers in Figures 8.3-16, 8.3-18, and 8.3-20 and note that the deeper the event, the higher the velocity and the shorter the cable length, the poorer the resolving power of curvature analysis for velocity determination. This observation is comparable to the case of conventional stacking velocity analysis, and it also applies to residual moveout analysis of image gathers derived from common-offset migration.

Compare the amplitudes on the images from depth migration of the zero-offset and prestack data (Figure 8.3-21). Note the differences in amplitude distribution along the flat reflector below the salt diapir. Imaging beneath complex structures has certain implications as to the acquisition geometry — specifically, on the choice of the cable length.

### Field data examples

We now review results of shot-profile migration of a 2-D land and a 2-D marine line. Figure 8.3-22 shows the CMP-stacked section of the land line over an imbricate structure associated with overthrust tectonics. The topographic high is the surface expression of the imbricate structure. Refraction and residual statics corrections were applied to the data before shot-profile migration.

Figure 8.3-22 also shows the poststack time-migrated section. The main thrust plane starts from the surface at CMP 340 and dips down to the right. To account for the severe elevation differences, migration was performed from the irregular topography (further aspects of migration in practice).

The recording geometry for this land line has severe irregularities with duplicate shots and missing shots. Here, we only want to discuss results of prestack depth migration and leave out the velocity-depth model estimation. Selected common-receiver gathers sorted from the migrated shot records are shown in Figure 8.3-23. Most events exhibit flat character. Note also the symmetry of some events with respect to the receiver location as seen on gathers in the top row — they are associated with flat or gently dipping reflectors. The asymmetry of some events with respect to the receiver location is seen on gathers in the bottom row in Figure 8.3-23.

Stacking of the common-receiver gathers yields the depth image obtained from shot-profile migration (Figure 8.3-24). For comparison, the depth image from poststack depth migration using the same velocity-depth model is shown also in Figure 8.3-24. Trace spacing in the prestack depth-migrated section corresponds to the average shot spacing, whereas that in the poststack depth-migrated section corresponds to the CMP spacing. Note from the prestack depth-migrated section (Figure 8.3-24) the improved image below 2 km of the imbricate structure to the right of the main thrust plane which starts at the surface in the vicinity of CMP 350 and dips down to the right.

Figure 8.3-25 shows the CMP-stacked section from the marine line over two salt diapirs. The target zone is base-salt and the subsalt region. For comparison, Figure 8.3-25 also shows the poststack time-migrated section. Note the velocity pull-up along the base-salt reflection at 1.9 s below CMP 400 and CMP 800 — typical of salt diapirs.

Selected common-receiver gathers sorted from the migrated shot records are shown in Figure 8.3-26. Note that most events exhibit flat character; this suggests that the velocity-depth model used in prestack depth migration is fairly accurate. Flat or gently dipping reflectors manifest themselves as events on the common-receiver gathers that are symmetric with respect to the receiver locations. Dipping reflectors, on the other hand, manifest themselves as events that are asymmetric with respect to the receiver locations.

Stacking of the common-receiver gathers yields the depth image obtained from shot-profile migration (Figure 8.3-27). For comparison, the depth image from poststack depth migration using the same velocity-depth model also is shown in Figure 8.3-27. Trace spacing in the prestack depth-migrated section corresponds to the shot spacing, whereas that in the poststack depth-migrated section corresponds to the CMP spacing. Note from the prestack depth-migrated section the improved image of the base-salt reflector and the subsalt region. The base-salt reflector appears to be faulted and shows a structural closure at approximately 3.5 km below CMP 630.

## 8.4 3-D poststack depth migration

The fundamentals of 3-D migration are discussed in 3-D poststack migration and its mathematical aspects are provided in Appendix G. By now, we should be familiar with one-pass and two-pass, implicit and explicit 3-D time migration algorithms. In this section, we shall examine aspects of 3-D poststack migration within the context of imaging beneath complex structures: (a) 3-D poststack time versus depth migration, (b) one-pass versus two-pass 3-D poststack depth migration, and (c) implicit versus explicit 3-D poststack depth migration. In the next section, we shall extend the analysis to 3-D depth migration of prestack data.

### 3-D poststack time versus depth migration

Consider a 3-D survey over a hypothetical salt-dome structure. The base map is shown in Figure 8.4-1. The synthetic 3-D survey data consist of 481 inlines and 481 crosslines with 25-m trace spacing in both directions. Selected cross-sections of the 3-D velocity-depth model and the associated zero-offset inline sections are shown in Figure 8.4-2. The salt dome has a circular symmetry with a flat base. The base map shown in Figure 8.4-1 actually is a time slice at 1200 ms to show the circular symmetry.

Keep in mind that the sections in Figure 8.4-2 are the cross-sections of the 3-D zero-offset wavefield along the traverses that coincide with the cross-sections of the 3-D velocity-depth model. In other words, we are not starting with the cross-sections of the velocity-depth model in Figure 8.4-2 and generating a set of 2-D zero-offset wavefields from them. Strictly speaking, the only section that does not contain any sideswipe energy is the center line I-241. Note the energy reflecting off the flank of the salt dome and being recorded on the lines away from the center line.

First, we perform one-pass implicit 3-D time migration (Section G.1) on the entire 3-D zero-offset synthetic data and display the same lines after 3-D time migration (Figure 8.4-3). The migration velocity field is based on the true subsurface velocity-depth model used in computing the 3-D zero-offset wavefield (Figure 8.4-2). Also, for comparison, pretend that two lines — center line I-241 and line I-181 away from the center, are part of a 2-D survey and perform 2-D time migration. Note the following observations:

1. The top of the salt, albeit the vertical scale is in time, has been imaged properly with 3-D time migration, while the base of the salt has not. This is because the salt diapir acts as a complex overburden.
2. The 2-D time migration produced the correct result for the top of the salt only along the center inline (I-241), because there are no sideswipes on this line; therefore, there is no need for 3-D migration. However, on a line away from the center line, say I-181, even the top of the salt has not been imaged properly by 2-D time migration, let alone the base of the salt. This is because the line contains sideswipes off the flank of the salt dome.
3. Now you may have second thoughts about 2-D seismic exploration. Fortunately, most structures have dominant strike and dip directions. Properly oriented 2-D migrated lines often yield an acceptably accurate structural picture. However, one lesson to be learned from comparing the 3-D and 2-D time migrations of line I-181 shown in Figure 8.4-3 is that the migration velocities that yield an acceptable 2-D migrated section may be quite different from the true subsurface velocity model required by 3-D migration.

We now perform one-pass implicit 3-D depth migration (Section G.1) of the 3-D zero-offset synthetic data shown in Figure 8.4-2. After 3-D depth migration of the entire volume of the 3-D zero-offset wavefield using the true 3-D velocity-depth model, we display the same lines (Figure 8.4-4). Also, for comparison, pretend that the two lines — center line I-241 and line I-181 away from the center, are part of a 2-D survey and perform 2-D depth migration, again using the correct velocity-depth model. Note the following observations:

1. The top and base of the salt now have been imaged properly with 3-D depth migration. The complex overburden problem has been solved. (Compare the left column in Figure 8.4-3 with that in Figure 8.4-4.)
2. The 2-D depth migration produced the correct subsurface model only along the center inline (I-241) because there are no sideswipes on this line. However, on a line away from the center line, say I-181, neither the top nor the base of the salt have been imaged properly. This is because the line contains sideswipes off the flank of the salt dome.
3. By performing 2-D depth migration iteratively to converge to assumingly correct depth model, we may be forcing the model to converge to something very different from the truth. This stems from the treatment of the sideswipes as events that are in the plane of profile.

### Two-pass versus one-pass 3-D poststack depth migration

In 3-D poststack migration, we examined the characteristics of two-pass and one-pass 3-D poststack time migration. We concluded that given the choice between one-pass and two-pass schemes, 3-D poststack time migration may be done in two passes provided the vertical velocity gradient is not excessively large and dips are gentle. Is the two-pass strategy also applicable to 3-D poststack depth migration? We shall experiment with the salt-dome synthetic data set (Figure 8.4-2) and conclude that 3-D poststack depth migration has to be done in one pass.

First, we consider two-pass implicit finite-difference 3-D time migration of the salt-dome synthetic data set. Start with the 3-D zero-offset wavefield (Figure 8.4-2), and apply time migration in the inline direction and display four selected inline sections and four selected crossline sections (Figure 8.4-5). After this first-pass migration, the imaging of the center inline (I-241) is complete, since there is no sideswipe energy that needs to be moved out of the plane of this inline. On the other hand, the center crossline (X-241) has not been imaged at all, because no movement of energy took place in this line as yet after the first-pass migration. On other lines, there is still some more imaging to be done — note, for instance, the sideswipe energy on inline I-151.

Now sort the data into crosslines, perform the second-pass migration on the already migrated data, and display the same inline and crossline sections (Figure 8.4-6). These are now the inlines and crosslines after the two-pass 3-D time migration. Since the overburden velocity above the salt layer is constant, the two-pass 3-D time migration correctly images the top-salt boundary. Since we have done time migration, whether it is two-pass or one-pass, we have not imaged the base-salt boundary correctly.

Figure 8.4-7 gives a summary of the results of one-pass (as in Figure 8.4-3) and two-pass (as in Figure 8.4-6) implicit 3-D time migration of the salt-dome model data (Figure 8.4-2). Although not evident from the inline sections in Figure 8.4-7, the top-salt boundary actually has been imaged more accurately by the two-pass migration as compared with the one-pass migration. The constant-velocity of the layer above the salt enables the two-pass scheme to produce an accurate result (Section G.1). On the other hand, even though the overburden velocity above the salt layer is constant, the approximation made in the splitting has caused the one-pass scheme to produce an incorrect image of the top-salt boundary, especially along the directions diagonal to inline and crossline directions (3-D poststack migration). Again, the base-salt boundary has not been imaged correctly by either one-pass or two-pass migrations, since we have done time migration rather than depth migration.

We now examine the plausability of a hybrid two-pass 3-D migration of the salt-dome synthetic data set, wherein the first-pass migration is in time and the second-pass migration is in depth. In areas with a complex overburden structure, such as the overthrust belts, usually velocity varies laterally more in the dip direction perpendicular to the thrust fronts than the strike direction. When that is the case, it might be plausable to do time migration in the strike direction with mild lateral velocity variations, followed by depth migrations of selected lines in the dip direction with strong lateral velocity variations. Such two-step hybrid strategy may be useful in building the 3-D velocity-depth model for a subsequent, proper 3-D depth migration of the 3-D data.

For the synthetic data set in Figure 8.4-2, start with the results of time migration in the inline direction (Figure 8.4-5) and perform depth migration in the crossline direction. Then, display the same inlines and crosslines as in Figure 8.4-6. After time migration as the first-pass and depth migration as the second-pass (Figure 8.4-8), note that we certainly have restored the true geometry of the top-salt boundary correctly. This is because the overburden is constant velocity; therefore, it did not matter whether we did time or depth migration. However, we have not been able to restore the geometry of the base-salt boundary, except for the center crossline (X-241), because this salt structure is truly 3-D in character with no dominant strike or dip direction.

Figure 8.4-9 gives a summary of the results of one-pass (as in Figure 8.4-4) and two-pass (as in Figure 8.4-8) 3-D depth migration of the salt-dome model data (Figure 8.4-2). Again, note that the top-salt boundary actually has been imaged more accurately by the two-pass migration as compared with the one-pass migration. The constant-velocity of the layer above the salt enables the two-pass scheme to produce an accurate result. On the other hand, even though the overburden velocity above the salt layer is constant, the approximation made in the splitting has caused the one-pass scheme to produce an incorrect image of the top-salt boundary. The base-salt boundary, however, has not been imaged correctly by the two-pass scheme, while it has been imaged correctly by the one-pass scheme.

The experiments using the salt-dome synthetic data set (Figure 8.4-2) combined with the model experiments presented in 3-D poststack migration lead us to the following conclusions:

1. If the velocity field is judged to be suitable for time migration, the two-pass strategy for 3-D time migration may be acceptable provided the vertical velocity gradient is not excessively large and dips are not very steep.
2. If the velocity field requires depth migration, the one-pass strategy for 3-D depth migration is imperative.

Figures 8.4-10 and 8.4-11 show a field data example of 3-D depth migration using the one-pass scheme. The inline (top left) and crossline (top right) stacked sections in both figures show a structural high, which is associated with culminations in an overthrust belt. Figure 8.4-12a is an inline cross-section of the 3-D velocity-depth model. The deepest horizon on the velocity model (horizon 8) corresponds to the event below 2 s on the inline stacked section in Figure 8.4-10. The true geometry of this horizon has been distorted by the structure above acting as a complex overburden. Although in two dimensions, the image-ray plot (Figure 8.4-12b) through the velocity-depth model (Figure 8.4-12a) verifies the presence of a complex overburden above horizon 8. (See 2-D poststack depth migration for a discussion on image rays.)

Note the similarity between the inline velocity-depth model (Figure 8.4-12a) and line B after 3-D depth migration (bottom left, Figure 8.4-10). Despite this similarity, there are parts of the survey in which the data are overmigrated (on line A, bottom left, Figure 8.4-11) and parts in which the data are undermigrated (not shown). Differences between the output from depth migration and the velocity-depth model used requires an iterative modification of the velocity-depth model where it departs from the output of depth migration (2-D poststack depth migration). For comparison, Figures 8.4-10 and 8.4-11 show 2-D depth migrations of the selected lines. Note the obvious difference between 2-D and 3-D depth migrations of line D (Figure 8.4-11). Although the 2-D depth-migrated section contains an abundance of reflections as compared to the 3-D depth migrated section, that reflection energy does not belong on line D. From line A or B, note that energy should be migrated in the updip direction from line D to C. Hence, after 3-D migration, line D is depleted of reflection energy (Figure 8.4-11), while line C is enriched (Figure 8.4-10).

The question of how to supply a 3-D velocity field to do 3-D depth migration is beyond the scope of this chapter. Earth modeling in depth is discussed in earth modeling in depth.

### Implicit versus explicit 3-D poststack depth migration

Finally, we compare the performance of the implicit and explicit schemes using the circularly symmetric salt-dome model data set of Figure 8.4-2. The implicit scheme uses the 45-degree extrapolator in a split mode, and the explicit scheme uses a one-dimensional (1-D) explicit filter combined with the 5 × 5 McClellan filter template (Section G.2). Selected inline sections are shown in Figure 8.4-13, and depth slices are shown in Figure 8.4-14. The top-salt boundary is imaged more accurately by the explicit scheme because of the near-circular symmetry of its impulse response (Figure 7.3-14). Positioning errors by the implicit scheme implied by its impulse response (Figure 7.3-3) are better observed on the depth slices. Note in Figure 8.4-14 that the top-salt boundary image by the implicit scheme is not circularly symmetric; instead, significant undermigration especially along the two diagonals has resulted. In contrast, the explicit scheme has preserved the circular character of the salt dome.

### 3-D poststack datuming

As a byproduct of 3-D migration algorithms based on wavefield extrapolation, 3-D stacked data can be datumed from the surface to a specified depth or time level in a 3-D sense. This capability can be particularly useful in reservoir studies. Basically, the surface wavefield is downward continued to a desired depth without invoking the imaging principle along the way. Figure 8.4-15 shows the 3-D zero-offset data (Figure 8.4-2) datumed from the surface z = 0 to a 1000-m depth. The following conclusions can be made:

1. Since there are no sideswipes on the center line, I-241, datuming in a 2-D or 3-D sense is identical.
2. However, there is a significant difference between 2-D datuming and 3-D datuming for lines with sideswipe energy — those that are increasingly farther from the center line, say, Line I-181.

The wave-equation datuming described here takes the input wavefield from one flat constant horizontal datum to another. Prestack and poststack wave-equation datuming using arbitrary 2-D datum interfaces is discussed in layer replacement. The constant datum level should not be a limitation, particularly in reservoir studies. For example, the 3-D stacked data can be datumed to the top of the reservoir level followed by detailed imaging of the target zone only.

## 8.5 3-D prestack depth migration

We remind ourselves of the two prominent circumstances that require migration of seismic data before stack and in three dimensions:

1. The compelling reason for doing time migration is dipping events. The compelling reason for doing time migration before stack is to account for conflicting dips with different stacking velocities. And the compelling reason for doing time migration in three dimensions is to account for the 3-D behavior of fault-plane reflections and reflections within fault blocks that give rise to the problem of conflicting dips. In a strict theoretical sense, in the presence of conflicting dips with different stacking velocities, you need to image the subsurface by migration of seismic data in time, before stack and in three dimensions.
2. The compelling reason for doing depth migration is lateral velocity variations. The compelling reason for doing depth migration before stack is to account for nonhyperbolic moveout caused by lateral velocity variations. And the compelling reason for doing depth migration in three dimensions is to account for the 3-D behavior of complex overburden structures that give rise to lateral velocity variations. In a strict theoretical sense, in the presence of lateral velocity variations, you need to image the subsurface by migration of seismic data in depth, before stack and in three dimensions.

Figures 8.5-1 and 8.5-2 show subsurface images along the same inline traverse obtained from 2-D and 3-D, post- and prestack, and time and depth migrations. Examine the quality of imaging of the fault planes and the reflectors within the faulted zone in each of the sections. The fault-plane reflections and the gently dipping reflections within the fault blocks give rise to the problem of conflicting dips, and the lateral velocity variations across the fault blocks themselves give rise to the problem of nonhyperbolic moveout. These two problems also manifest themselves in three dimensions. Therefore, to get the best possible image, as predicated by the data set in Figures 8.5-1 and 8.5-2, one should do the migration before stack, in depth, and in three dimensions.

What then is an appropriate algorithm for 3-D prestack depth migration — Kirchhoff-summation, finite-difference, or frequency-wavenumber? The algorithm of choice must meet the following requirements:

1. For it to qualify as a depth migration algorithm, first and foremost, the algorithm of our choice must be able to image steeply dipping reflectors in the presence of lateral velocity variations.
2. 3-D prestack seismic data invariably suffer from irregular spatial sampling. Therefore, the algorithm of choice for prestack migration must cope with irregularly sampled data.
3. Just as we use prestack time migration to estimate rms velocities (migration velocity analysis and 3-D prestack time migration), we often wish to use prestack depth migration to estimate interval velocities. When used as a velocity estimation tool, we do not have to generate image gathers from prestack depth migration at all CMP locations along a 2-D line or at all bin locations over a 3-D survey area. Instead, it often is sufficient to generate image gathers at sparse intervals along the line, or along selected lines or even on a sparse grid over the 3-D survey area.

Whatever the type of algorithm, requirement (a) cannot be waived for depth migration. Requirement (b) to handle irregular spatial sampling may be fulfilled by azimuth-moveout (AMO) correction [14] or inversion to common offset (ICO) [15]. Once data are spatially regularized so that the resulting prestack data have uniform fold of coverage and have been corrected for source-receiver azimuth, then, in principle, any of the three categories of migration algorithms — Kirchhoff-summation, finite-difference, or frequency-wavenumber, can be used to perform 3-D prestack depth migration. While it has the advantage of producing an image in depth along a set of line traverses without having to produce an image volume in depth, Kirchhoff summation technique lacks the rigor to handle amplitudes that the frequency-wavenumber techniques can provide. Finally, the finite-difference and frequency-wavenumber migration methods are global methods; as such, they are not suitable to meet requirement (c). The better treatment of amplitudes by the frquency-wavenumber algorithms compared to the Kirchhoff summation technique, however, has greatly increased their use in practice [16].

### Kirchhoff summation

Kirchhoff’s integral solution to the scalar wave equation

 ${\displaystyle {\frac {\partial ^{2}P}{\partial x^{2}}}+{\frac {\partial ^{2}P}{\partial y^{2}}}+{\frac {\partial ^{2}P}{\partial z^{2}}}={\frac {1}{v^{2}(x,y,z)}}{\frac {\partial ^{2}P}{\partial t^{2}}}}$ (1)

gives the pressure wavefield P(x, y, z; t) propagating in a medium with velocity v(x, y, z) at a location (x, y, z) and at an instant of time t. As such, the Kirchhoff solution is a mathematical statement of Huygen’s principle which states that the pressure disturbance at time tt is the superposition of the spherical waves generated by point sources at time t [17].

The discrete form of the integral solution to equation (8-1) as used in practical implementation of Kirchhoff migration is given by (Section H.1)

 ${\displaystyle P_{out}={\frac {\Delta x\Delta y}{4\pi }}\sum _{A}{\frac {\cos \theta }{vr}}{\frac {\partial }{\partial t}}P_{in},}$ (2)

where Δx and Δy are inline and crossline trace intervals, Pout = p(xout, yout, z; τ = 2z/v) is the output of migration using the input wavefield Pin = P(xin, yin, z = 0; τ = t − r/v) within an areal aperture A. The geometry associated with equation (8-2) is given in Figure H-1.

The Kirchhoff summation method requires

1. computing nonzero-offset traveltimes through a 3-D, spatially varying velocity medium, and
2. scaling and summation of the amplitudes along the computed traveltime trajectory based on the Kirchhoff integral solution to the scalar wave equation.

The scaling of amplitudes before summation includes application of the obliquity factor cos θ, the spherical divergence factor 1/vr, and the amplitude and phase corrections, |ω| exp(/2), implied by the derivative operator ∂P/∂t (migration principles). Additionally, as for any migration, undersampling of the input data in x and y directions needs to be compensated for by a suitable antialiasing filter. Nevertheless, it is the traveltime computation that poses numerical accuracy and efficiency challenges when implementing the Kirchhoff summation method for 3-D prestack depth migration.

### Calculation of traveltimes

A direct method to compute traveltimes is ray tracing through the specified velocity-depth model. A bundle of rays emerging from a source location at the surface can be sprayed down into the earth and traced through the subsurface while accounting for ray bending caused by changes in velocity gradient and refraction at layer boundaries with velocity contrast. Reflection points along each of the raypaths are identified as the intersection points of the rays with the layer boundaries. The traveltime from the source location at the surface and a reflection point at the subsurface is then calculated by integrating the elements of distance along the raypath divided by the velocity associated with that element. By applying reciprocity, the traveltime from a receiver location at the surface to a reflection point in the subsurface can be computed in the same manner. Finally, for a given source-receiver pair at the surface and a reflection point in the subsurface, the total traveltime is computed by adding the traveltime from the source to the reflection point to the traveltime from the reflection point to the receiver.

The two-point ray tracing described above is conceptually simple, but is computationally intensive. Efficient ray tracing through complex velocity-depth models is not a trivial task. Alternatives to two-point ray tracing, however, have been developed and implemented with sufficient accuracy. Examples include paraxial ray tracing [18] [19] and Gaussian beam ray tracing [20].

Given the circumstances described above, one can identify potential problems with the ray tracing. There will not always be a raypath combination associated with a source-receiver pair and a reflection point. Depending on the recording geometry and the complexity of the velocity-depth model, there may be zones through which some rays may be missed. This complication is compounded by the computational load involved in ray tracing itself. As such, direct ray tracing is rarely used for traveltime computations required for prestack depth migration.

Figure 8.5-3 shows the traveltime contours through a velocity-depth model that includes a salt sill with velocity higher than the surrounding sediments. The discontinuities along the traveltime contours correspond to locations where ray bending is implied by ray tracing. Note, however, that the locations where ray bending takes place do not coincide with the salt sill boundary where the largest velocity contrast exists. Problems with ray tracing through a complex model as outlined above also cause physically implausable rapid variations in the traveltime contours.

An alternative to two-paint ray tracing is wavefront construction [21], which involves tracing not just one ray but a fan of rays together. As such, the medium represented by the velocity-depth model used in ray tracing is covered adequately by controlling the ray density along wavefronts [22]. In areas with low ray density, additional ray bundles may be created by paraxial rays.

### The eikonal equation

Consider a plane wave function P(x, y, z; t) with a spatially varying amplitude P0(x, y, z) and spatially varying traveltime T(x, y, z) (Section H.2)

 ${\displaystyle P(x,y,z;t)=P_{0}(x,y,z)\exp {\Big \{}-i\omega {\Big [}t-T(x,y,z){\Big ]}{\Big \}}.}$ (3)

Assume for the moment that the wave amplitude P0(x, y, z) does not vary spatially, but is a constant. Then, the plane-wave solution of the form defined by equation (8-3) satisfies the scalar wave equation (8-1). Substitution of equation (8-3) into equation (8-1) yields

 ${\displaystyle {\left({\frac {\partial T}{\partial x}}\right)}^{2}+{\left({\frac {\partial T}{\partial y}}\right)}^{2}+{\left({\frac {\partial T}{\partial z}}\right)}^{2}={\frac {1}{v^{2}(x,y,z)}},}$ (4)

which is called the eikonal equation. It is a ray-theoretical approximation to the scalar wave equation (8-1). Derivation of equation (8-4) is provided in Section H.2.

While the solution of the scalar wave equation (8-1) represents the wavefield P(x, y, z; t) at a point in space (x, y, z) and at an instant of time t, the solution of the eikonal equation (8-4) represents the traveltime T(x, y, z) for a ray passing through a point (x, y, z) in a medium with velocity v(x, y, z). Specifically, T(x, y, z) = constant represents the wavefront of constant phase at an instant of time. The wave is propagated from one wavefront to the next by way of raypaths, which are perpendicular to the wavefronts.

When the medium velocity is not constant but is an arbitrary function of space variables v(x, y, z), and the wave amplitude P0(x, y, z) is not constant but also varies spatially, then the traveltime function T(x, y, z) of equation (8-3) is not a solution to the eikonal equation (8-4). For a wave function, such as that given by equation (8-3), with spatially varying amplitudes, the eikonal equation is a good approximation to the wave equation only at a high-frequency limit (Section H.2). The high-frequency limit is equivalent to small wavelengths. How small should the wavelength be for a solution of the eikonal equation to be a good approximation to the wave equation? The approximation is valid if the fractional change in the velocity gradient is much less than the wave frequency [17]. In practice, this means that the eikonal equation can be used to compute traveltimes if the velocity-depth model does not contain large velocity gradients.

Just as a solution to the scalar wave equation can be formulated using a finite-difference scheme (migration principles), the eikonal equation (8-4) also can be solved using the finite-difference technique (Section H.3). Consider the 3-D computation mesh sketched in Figure H-2b. We want to compute the traveltime T of the eikonal equation (H-27) at grid point (x + Δx, y, z + Δz) using the known traveltimes at grid points (x, y, z), (x + Δx, y, z), and (x + Δx, y + Δy, z).

Computing the traveltimes at depth z + Δz from those at depth z means extrapolating T in the z-direction. Rewrite equation (8-4) in the form of an extrapolation equation as [10]

 ${\displaystyle v{\frac {\partial T}{\partial z}}={\sqrt {1-\left(v{\frac {\partial T}{\partial x}}\right)^{2}-\left(v{\frac {\partial T}{\partial y}}\right)^{2}}}.}$ (5)

As shown in Section H.3, the 3-D eikonal equation (8-5) can be solved for the traveltime values T(x, y, z) for the propagating wavefront through a velocity field v(x, y, z) in the subsurface using a finite-difference scheme.

Figure 8.5-4 shows the traveltime contours through the same velocity-depth model as in Figure 8.5-3 computed by using the eikonal equation. These may be considered as wavefronts expanding from a source at the surface located at 3000 m from the origin. Note that the shape of the traveltime contours honors the boundaries with velocity contrast. Nevertheless, the traveltime contours do not show the abrupt changes as would be the case at a boundary such as the salt-sediment interface with sharp velocity contrast.

To explain the smooth behavior of the traveltime contours at layer boundaries in Figure 8.5-4, refer to the sketch in Figure 8.5-5 of an expanding wavefront through a flat interface with a large velocity contrast. At the critical angle of refraction, waves travel along the layer boundary with the faster velocity of the underlying layer, Eventually, these waves are refracted back into the overlying layer and are recorded in the form of first arrivals; thus, they often are called head waves. Note that the wavefront associated with the head wave tends to smooth out the sharp change in the expanding wavefront as it crosses over the layer boundary with velocity contrast. It is this effect that we see in Figure 8.5-4 where the traveltime contours make a smooth transition from the sediment to the salt layer.

Ideally, the traveltime contours should pronounce the sharp layer boundaries. While some solutions to the eikonal equation include the head waves [24] [25] as shown in Figure 8.5-4, others exclude the head wave [10] [26]. Figure 8.5-6 shows the traveltime contours derived from the solution of the eikonal equation that has excluded the head wave. Note the sharp changes in the traveltime contours, which at some locations appear to coincide nicely with the salt-sediment layer boundary. Nevertheless, the match is not consistent at all locations and the traveltime contours show rapid fluctuations that are not physically plausable.

Figure 8.5-7a shows traveltime contours that were computed using a finite-difference solution to the eikonal equation. The velocity-depth model consists of layers with velocities that vary from 1500 m/s to 4000 m/s. Note that the eikonal solution, although it is associated only with the fastest arrival, always yields a raypath from a source at the surface to a point in the subsurface through the gridded velocity-depth model. Some of the raypaths are indicated by trajectories that have a solid circle at the end. The smooth behavior of the traveltime contours is associated with the head wave that often is the first arrival.

In contrast with the eikonal solution, the wavefront construction yields multiple arrivals as shown in Figure 8.5-7b. Again, some of the raypaths are indicated by trajectories that have a solid circle at the end. Note the multiple raypaths associated with a single point in the subsurface. Where a head wave or diffraction develops, the wavefront construction leaves a gap in the traveltime contour. With wavefront construction, however, you can be selective in the type of rays that you would like to include in the summation of amplitudes. Figure 8.5-7c shows traveltime contours associated with the reflecting boundary represented by the thick curve. In fact, the solution from wavefront construction in this case includes arrivals associated with reflected waves as well as transmitted waves.

The eikonal solution does not necessarily yield the maximum energy along the single-arrival raypath. In wavefront construction, on the other hand, amplitudes associated with multiple arrivals can be included in the summation. Thus, the chance of capturing the maximum-energy arrival increases, albeit at an increase in computational cost. By including multiple arrivals in the summation, the likelihood of attaining a complete image of a complex structure also is higher. Figure 8.5-8a shows raypaths from a single source at the surface that illuminate portions of the steep flank of a salt dome. By using the arrivals associated with these raypaths, only, in Kirchhoff summation, the resulting partial image shown in Figure 8.5-8b is obtained. By including in the summation the arrivals associated with the reflections off the salt flank as well as the arrivals associated with the reflections from a deeper interface that bounce off the salt flank (Figure 8.5-8c), a more complete image can be obtained (Figure 8.5-8d).

### Fermat’s principle

The traveltime along a raypath from one point to another has an extremum value which, for most physical problems, is a minimum. This is the formal statement of Fermat’s principle [17]. The raypath also defines the direction of energy flow. Among a bundle of rays from one point to another, Fermat’s principle can be applied to discard all but one raypath that corresponds to a minimum time of travel from one point to the other. This practical concept can be used to perform traveltime computations for prestack depth migration [27] [28].

Consider the raypath geometry depicted in Figure 8.5-9. The velocity-depth model may be split into a set of horizontal slabs with a specified thickness, say 50 to 300 m. Assume that traveltimes from a source at the surface z = 0 to a depth z through a velocity-depth model already have been computed. Now compute the traveltime from the grid point 6 at depth z to each of the grid points at depth z + Δz within a specified aperture, say grid points 1 to 11. An average velocity between the grid points at depth z and z + Δz along each of the raypaths may be used in the computation. Among the 11 raypaths from the grid point 6 at depth z to grid points 1 to 11 at depth z + Δz, choose that which corresponds to the minimum traveltime.

The process may be continued for all grid points at depth z, and then from one depth to the next, and the traveltimes associated with the minimum-traveltime raypaths through each of the horizontal slabs may be added to compute the total traveltime from a source or receiver point at the surface z = 0 to a reflection point at some depth in the subsurface.

### Summation strategies

Figure 8.5-10 shows the zero-offset wavefield responses of a point diffractor buried in media with varying degrees of complexity. The traveltime trajectories associated with the point diffractors buried in a constant-velocity medium, beneath an overburden with mild to moderate lateral velocity variations, and beneath an overburden with strong lateral velocity variations, all are single-valued. Therefore, ray tracing through such models would produce unambiguous traveltimes for Kirchhoff summation.

The zero-offset traveltime trajectory associated with the point diffractor buried beneath an overburden with severe lateral velocity variations, however, is multivalued (Figures 8.5-10d). One would have to decide as to what summation path to choose:

1. The travelpath that corresponds to the first arrivals; that is, minimum-time summation trajectory [24] [25],
2. The travelpath that corresponds to the bowties that contain the most significant portion of the energy associated with the zero-offset wavefield response; that is, maximum-energy summation trajectory [23],
3. The travelpath that corresponds to the shortest distance between the source or receiver point at the surface and the reflection point at the subsurface; that is, minimum-distance summation trajectory [29], or
4. The entire multivalued travelpath.

From the zero-offset wavefields associated with a diffractor buried in a medium with varying complexity shown in Figure 8.5-10, it can be inferred that the minimum-time strategy may be suitable for cases of moderate to strong lateral velocity variations (Figure 8.5-10b,c), whereas the maximum-energy strategy may be imperative for a case of a complex overburden with severe lateral velocity variations (Figure 8.5-10d). Ideally, it would be desirable not to exclude any portion of the traveltime trajectory and use a multivalued summation path. This, however, can be formidably costly and often is not needed in practice. Efficient traveltime calculation and choice of a summation path are important considerations for the 3-D prestack depth migration of large volumes of seismic data [27] [30].

### Migration aperture

Because of cost considerations in 3-D prestack depth migration, it is compelling to make a careful choice of aperture width in the inline and crossline directions. While an excessively large aperture unnecessarily increases run time, a small aperture can produce a poor image from 3-D prestack depth migration. Migration aperture in the Kirchhoff summation was discussed already in kirchhoff migration in practice. Recall that a small aperture causes destruction of steeply dipping events. Excessively small aperture width also organizes random noise, especially in the deeper part of the section, as horizontally dominant spurious events.

Figure 8.5-11 shows aperture tests for 3-D prestack depth migration. A small aperture in both inline and crossline directions causes poor imaging of the steeply dipping fault planes and the steep reflector that defines the base of the sedimentary basin. Note also the organizing effect of small aperture on random noise.

### Operator antialiasing

In the 2-D Fourier transform we reviewed the phenomenon of spatial aliasing caused by undersampling of recorded data, and in further aspects of migration in practice we examined its adverse effects on migration. Prestack data may be made spatially uniform, for instance, by azimuth-moveout correction (AMO) [14], thus enabling use of finite-difference or frequency-wavenumber algorithms for prestack depth migration. When it becomes necessary to use Kirchhoff summation algorithm, operator aliasing is an issue that needs to be dealt with.

Shown in Figure 8.5-12 is a low-velocity hyperbolic summation path. Note that it intercepts more than one sample in a given input trace, each sample indicated by a shaded grid block. Thus, the summation along this path should include more than one sample per input trace. There are three different approaches to handle operator aliasing caused by mulitple samples per input trace included in the summation:

1. The summation trajectory that represents the kinematics of the migration operator may be truncated to exclude the steep flanks that suffer from aliasing.
2. Trace interpolation may be used to create additional traces so as to avoid multiple samples per input trace included in the summation (processing of 3-D seismic data).
3. The frequency components that are aliased at a given dip along the summation trajectory are filtered out [32].

The first approach is undesirable since truncation is equivalent to limiting the migration aperture which can destroy steeply dipping events (kirchhoff migration in practice). The second approach can be taken but at a high cost in the case of prestack data. Trace interpolation also suffers in accuracy when applied to data with dipping events that conflict with one another. The third approach requires multiple copies of the input trace with different band-widths; this can be troublesome for very large input prestack data. Lumley [33] applied a triangular filter [34] to input traces in the time domain without creating multiple copies. An improved form of the triangular filtering method for 3-D Kirchhoff summation is given by Abma [35].

### 3-D common-offset depth migration

In 3-D prestack time migration we defined a strategy for 3-D prestack time migration based on the application of 3-D zero-offset migration theory to 3-D common-offset data. Specifically, following NMO and 3-D DMO correction, prestack data will have been corrected for dip and source-receiver azimuth effects. As a result, the common-offset volumes of data are decoupled representations of a 3-D zero-offset wavefield; thus, each common-offset volume can be migrated independently using a 3-D zero-offset time migration algorithm. For the subsequent processing sequence, refer to 3-D prestack time migration.

Decoupling of common-offset data volumes, in a strict theoretical sense, is only possible for events with hyperbolic moveout. For events with nonhyperbolic moveout that calls for depth migration, each of the resulting common-offset volumes following the applications of NMO and 3-D DMO corrections would not be a representation of a 3-D zero-offset wavefield. As such, and again in theory, a 3-D common-offset volume of data could not be depth-migrated.

Despite this theoretical mandate, it is surprisingly pleasant to see that indeed, in practice, the common-offset strategy for 3-D prestack time migration sometimes can also be applicable to 3-D prestack depth migration. While, following NMO and 3-D DMO correction, each of the common-offset volumes is time-migrated using a 3-D rms velocity field, depth migration of the common-offset volumes is done using a 3-D velocity-depth model.

There are two reasons why you may choose to apply the common-offset strategy to 3-D prestack depth migration:

1. You may wish to use a 3-D zero-offset depth migration algorithm other than the Kirchhoff summation, and thus avoid the troublesome task of computing 3-D nonzero-offset traveltimes needed for the latter.
2. For line output or for selected image gathers Kirch-hoff summation is the appropriate algorithm. For volume output from 3-D prestack depth migration other algorithms can be more efficiently applied to common-offset data.

Figure 8.5-13 shows selected inline sections derived from 3-D prestack depth migration using the Kirchhoff summation and the common-offset strategy. For the latter, a 3-D zero-offset phase-shift-plus-correction (PSPC) [36] was used to migrate all 32 volumes of 3-D DMO-corrected common-offset data; thus a full volume of image in depth was created. Comparison of the sections from the two approaches indicate that the imaging quality is largely comparable.

The common-offset strategy for depth migration will not work strictly for a case of complex overburden structure that gives rise to complex, nonhyperbolic moveout. Shown in Figure 8.5-14a is the image obtained from prestack depth migration based on Kirchhoff summation with a finite-difference solution to the eikonal equation to compute the traveltimes. The synthetic data [37] are associated with a salt sill model (Section K.2). The objective is to image the subsalt region accurately. Note that the common-offset strategy for prestack depth migration with or without DMO correction (Figure 8.5-14b,c) included in the workflow fails to produce an accurate image of the subsalt region with a level of accuracy that is comparable to the result from the Kirchhoff summation (Figure 8.5-14a).

As a final note, efficient and accurate implementation of 3-D prestack depth migration is a topic of active and intensive research in the industry. A choice of strategies for 3-D prestack depth migration is summarized below:

1. You may wish to first apply azimuth-moveout (AMO) correction [14] to regularize the prestack data, then use a frequency-wavenumber algorithm [16] to perform the 3-D prestack depth migration. The method involves the design and application of a wavefield extrapolation operator on the whole of the 3-D prestack data without decoupling the common-offset volumes.
2. Alternatively, you may wish to perform 3-D prestack depth migration using the Kirchhoff summation technique discussed in this section. The Kirchhoff summation implicitly handles geometry irregularities associated with the prestack data.
3. A third strategy, which also is described in the latter part of this section, is based on migration of the decoupled 3-D common-offset volumes of data that have been corrected for geometry irregularities and reduced to zero offset by way of 3-D DMO correction. Again, the preferred choice for the algorithm to perform the 3-D common-offset depth migration is one based on wave extrapolation.

A general guideline for the choice of strategy for 3-D prestack depth migration is given below:

1. If you are dealing with a low-relief structure and moderate lateral velocity variations, use 3-D prestack depth migration for conducting image-gather analysis for model verification and update, but not necessarily for imaging in depth. For the purpose of model verification and update, you need to work with only a sparse grid of image gathers; and the Kirchhoff summation is the suitable algorithm to produce such a set of image gathers without creating a whole image volume. As for imaging in depth, you may be content with 3-D poststack depth migration; and a finite-difference (Sections G.1 and G.2) or frequency-wavenumber algorithm (3-D poststack migration) that uses wave extrapolation is the suitable algorithm.
2. If you are dealing with a complex structure and moderate-to-strong lateral velocity variations, use a 3-D prestack depth migration algorithm based on decoupling of common-offset volumes of data as demonstrated in this section.
3. Finally, if you are dealing with a complex overburden structure, use of a 3-D prestack depth migration algorithm without decoupling the common-offset data is imperative.

## Exercises

Exercise 8-1. Consider the zero-offset x − t section depicted in Figure 8.E-1. Using the velocity-depth model in the same figure, sketch the depth-migrated x − z section.

Exercise 8-2. Identify the prominent sideswipe energy on Lines I-181 and I-151 in Figure 8.4-2 (right column).

Exercise 8-3. In areas with complex overburden, which one is more important — 3-D poststack depth migration or 2-D prestack depth migration?

Exercise 8-4. How would the geometry of the subsalt target reflectors labeled in yellow in Figure 8.E-2 behave on a section from an image volume created by 3-D poststack depth migration?

## Appendix H: Diffraction and ray theory for wave propagation

### H.1 The kirchhoff integral

Kirchhoff’s integral solution to the scalar wave equation

 ${\displaystyle \left[{\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial x^{2}}}-{\frac {1}{v^{2}(x,y,z)}}{\frac {\partial ^{2}}{\partial t^{2}}}\right]P(x,y,z;t)=0}$ (H-1)

is a mathematical statement of Huygen’s principle. In equation (H-1), P(x, y, z; t) is the pressure wavefield propagating in a medium with velocity v(x, y, z). Huygen’s principle states that the pressure disturbance at time t + Δt is the superposition of the spherical waves generated by point sources at time t.

Consider the geometry in Figure H-1 of a point diffractor at a location S(x, y, z) and an observation surface A for the diffraction wavefield generated by the source at the diffractor. Actually, the surface area A is only a portion of a closed surface and is the aperture of the observation made over that closed surface. For convenience, we shall choose the receiver location R(0, 0, 0) on the observation surface A to be at the origin of our coordinate system.

Also for convenience, we shall apply Fourier transform to the wavefield in the time direction

 ${\displaystyle P(x,y,z,;\omega )=\int P(x,y,z;t)\ \ \exp(-i\omega t)\ dt,}$ (H-2a)

where ω is the angular frequency. The inverse transform is given by

 ${\displaystyle P(x,y,z,;t)=\int P(x,y,z;\omega )\ \ \exp(i\omega t)\ d\omega .}$ (H-2b)

Now apply Fourier transform to equation (H-1) in the time direction to obtain

 ${\displaystyle \left(\nabla ^{2}+{\frac {\omega ^{2}}{v^{2}}}\right)P(x,y,z;\omega )=0,}$ (H-3)

where ${\displaystyle \nabla ^{2}}$ is the Laplacian operator

${\displaystyle \nabla ^{2}P=\left({\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial x^{2}}}\right)P.}$

We may intuitively state that what we observe at the surface A is what is generated at the source S. This statement is mathematically expressed by Gauss’s divergence theorem as [17] [38] [39] [40]

 ${\displaystyle \int _{V}\nabla ^{2}PdV=\int _{A}{\frac {\partial P}{\partial n}}\ dA,}$ (H-4)

where V is the volume of the region enclosed by the surface A and the derivative ∂P/∂n is taken normal to the surface A in the outward direction. Note that Gauss’s divergence theorem expressed by equation (H-4) transforms a surface integral into a volume integral. For our specific application, this means imaging the earth’s interior from seismic waves observed at the earth’s surface.

We shall solve equation (H-3) for each frequency component ω and sum the resulting solutions over all frequency components to compute the wavefield at the source P(x, y, z; t = 0). The solution obtained by Kirchhoff in 1882 requires Green’s function which describes the propagation outward from a point source with spherical symmetry as [39] [38]

 ${\displaystyle G(r,\omega )={\frac {1}{r}}\ \ \exp \left(-i{\frac {\omega }{v}}r\right),}$ (H-5a)

where

 ${\displaystyle r={\sqrt {x^{2}+y^{2}+z^{2}}}}$ (H-5b)

is the distance between the observation point and the source location. Green’s function given by equation (H-5a) is also a valid solution to equation (H-3):

 ${\displaystyle \left(\nabla ^{2}+{\frac {\omega ^{2}}{v^{2}}}\right)G(x,y,z;\omega )=0.}$ (H-6)
Figure H-1  Geometry of a point diffractor to derive the Kirchhoff integral solution to the scalar wave equation.

We now rewrite equation (H-4) by multiplying both sides with Green’s function G of equation (H-5a) as

 ${\displaystyle \int _{V}G\nabla ^{2}PdV=\int _{A}G{\frac {\partial P}{\partial n}}\ dA,}$ (H-7a)

and rewrite, in turn, equation (H-7a) by interchanging our wave function P with Green’s function G to obtain

 ${\displaystyle \int _{V}P\nabla ^{2}GdV=\int _{A}P{\frac {\partial G}{\partial n}}dA.}$ (H-7b)

Now subtract equation (H-7b) from (H-7a) to get

 ${\displaystyle \int _{V}(G\nabla ^{2}P-P\nabla ^{2}G)\ dV=\int _{A}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dA.}$ (H-8)

Equation (H-8) is known as Green’s theorem. Albeit the algebraic steps that we took from Gauss’s theorem given by equation (H-4) to Green’s theorem given by equation (H-8) may appear informal, the intention here is to provide the basic concept underlying the derivation.

Substitute equations (H-3) and (H-6) into the left-hand side of equation (H-8)

${\displaystyle \int _{V}(G\nabla ^{2}P-P\nabla ^{2}G)\ dV=\int _{V}\left(-G{\frac {\omega ^{2}}{v^{2}}}P+P{\frac {\omega ^{2}}{v^{2}}}G\right)\ dV,}$

and hence note that

 ${\displaystyle \int _{V}(G\nabla ^{2}P-P\nabla ^{2}G)\ dV=0.}$ (H-9)

Now we turn our attention to the right-hand side of equation (H-8). Because Green’s function defined by equation (H-5a) becomes infinite at the source location S, we need to place it inside an infinitesimally small enclosed surface E. This would then require computing the right-hand side of equation (H-8) in two parts — once for the surface E and once for the surface A.

Substitute equation (H-5a) into the right-hand side of equation (H-8) and note, from Figure H-1, that for the surface E, (/∂n) = −(/∂r):

 ${\displaystyle \int _{E}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dE=\int _{E}\left[-{\frac {1}{r}}\ \exp {\Bigg (}-i{\frac {\omega }{v}}r{\Bigg )}{\frac {\partial P}{\partial r}}+P{\frac {\partial }{\partial r}}\left({\frac {1}{r}}\ \exp {\Bigg (}-i{\frac {\omega }{v}}r{\Bigg )}\right)\right]\ r^{2}d\Omega ,}$ (H-10a)

where dE = r2dΩ and Ω is the solid angle around the point source S in Figure H-1. Carry out the differentiation with respect to r and simplify the right-hand side of equation (H-10a) to obtain

 ${\displaystyle \int _{E}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dE=-\int _{E}\exp {\bigg (}-i{\frac {\omega }{v}}r{\bigg )}\left(r{\frac {\partial P}{\partial r}}+P+i{\frac {\omega }{v}}rP\right)\ d\Omega .}$ (H-10b)

Finally, take the limit r → 0 to obtain the contribution of the surface E:

 ${\displaystyle \int _{E}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dE=-4\pi P.}$ (H-11)

Now substitute, again, equation (H-5a) into the right-hand side of equation (H-8) and note from Figure H-1 that for the surface A, (/∂n) = −(/∂z):

 ${\displaystyle \int _{A}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dA=\int _{A}\left\{-{\frac {1}{r}}\ \exp \left(-i{\frac {\omega }{v}}r\right){\frac {\partial P}{\partial z}}+P{\frac {\partial }{\partial z}}\left[{\frac {1}{r}}\ \exp \left(-i{\frac {\omega }{v}}r\right)\right]\right\}\ dA.}$ (H-12)

Carry out the differentiation with respect to z while noting from Figure H-1 that ∂r/∂z = cos θ, and simplify the right-hand side of equation (H-12) to get

 ${\displaystyle \int _{A}\left(G{\frac {\partial P}{\partial n}}-P{\frac {\partial G}{\partial n}}\right)\ dA=-\int _{A}\exp \left(-i{\frac {\omega }{v}}r\right)\left({\frac {1}{r}}{\frac {\partial P}{\partial z}}+{\frac {\cos \theta }{r^{2}}}P+i{\frac {\omega }{v}}{\frac {\cos \theta }{r}}P\right)\ dA.}$ (H-13)

The total contribution to the right-hand side of equation (H-8) actually is the sum of equations (H-11) and (H-13). The left-hand side of equation (H-8) vanishes by way of equation (H-9). Hence, the resulting expression from equation (H-9) is

 ${\displaystyle 4\pi P=\int _{A}\exp \left(-i{\frac {\omega }{v}}r\right)\left({\frac {1}{r}}{\frac {\partial P}{\partial z}}+{\frac {\cos \theta }{r^{2}}}P+i{\frac {\omega }{v}}{\frac {\cos \theta }{r}}P\right)\ dA.}$ (H-14)

Recall that P = P(x, y, z; ω) in equation (H-14) and multiply both sides by exp(iωt) and integrate over frequency ω. The left-hand side of the equation then becomes P(x, y, z; t) by way of the inverse Fourier transform as in equation (H-2b). Thus the resulting expression is

 ${\displaystyle P(x,y,z;t)={\frac {1}{4\pi }}\int _{\omega }\int _{A}\left({\frac {1}{r}}{\frac {\partial P}{\partial z}}+{\frac {\cos \theta }{r^{2}}}P+i{\frac {\omega }{v}}{\frac {\cos \theta }{r}}P\right)\ \exp \left[-i\omega \left(t-{\frac {r}{v}}\right)\right]\ dA\ d\omega .}$ (H-15)

Define the variable τ = t − r/v as retarded time. Recall from Appendix A that if P(ω) is the Fourier transform of P(t), then exp(−iωr/v)P(ω) is the Fourier transform of P(τ = t − r/v). Also, if P(ω) is the Fourier transform of P(t), then iωP(ω) is the Fourier transform of ∂P/∂t. Incorporate these relations into equation (H-15) and apply inverse Fourier transform to the right-hand side to obtain

 ${\displaystyle P(x,y,z;\tau )={\frac {1}{4\pi }}\int _{A}\left\{{\frac {1}{r}}\left[{\frac {\partial P}{\partial z}}\right]+{\frac {\cos \theta }{r^{2}}}[P]+{\frac {\cos \theta }{vr}}\left[{\frac {\partial P}{\partial t}}\right]\right\}\ dA,}$ (H-16)

where [P] means that the integration over the area A is done using the wavefield P at retarded time τ = t − r/v.

The first term depends on the vertical gradient of the wavefield ∂P/∂z. The second term is called the near-field term since it decays with 1/r2. Both terms are neglected in seismic migration. The remaining third term is called the far-field term and it is the foundation of Kirchhoff migration. Writing it in discrete form, we have

 ${\displaystyle P_{out}={\frac {\Delta x\Delta y}{4\pi }}\sum _{A}{\frac {\cos \theta }{vr}}{\frac {\partial }{\partial t}}P_{in},}$ (H-17)

where Δx and Δy are inline and crossline trace spacings, Pout = P(xout, yout, z; τ = 2z/v) is the output of migration using the input wavefield Pin = P(xin, yin, z = 0; τ = t − r/v) within an areal aperture A.

### H.2 The eikonal equation

The eikonal equation is a ray-theoretical approximation to the scalar wave equation. Its solution represents wavefronts of constant phase. The wave is propagated from one wavefront to the next by way of raypaths which are perpendicular to the wavefronts.

Consider a compressional plane wave P(x, y, z; t) in 3-D Cartesian coordinates (x, y, z)

 ${\displaystyle P(x,y,z;t)=P_{0}\ \exp(-i\omega t+ik_{x}x+ik_{y}y+ik_{z}z),}$ (H-18)

where P0 is the wave amplitude, t is the traveltime, and kx, ky, kz and ω are the Fourier duals of the variables x, y, z and t, respectively. That this equation is a solution to the 3-D scalar wave equation

 ${\displaystyle {\frac {\partial ^{2}P}{\partial x^{2}}}+{\frac {\partial ^{2}P}{\partial y^{2}}}+{\frac {\partial ^{2}P}{\partial z^{2}}}={\frac {1}{v^{2}(x,y,z)}}{\frac {\partial ^{2}P}{\partial t^{2}}}}$ (H-19)

can be verified by computing the partial derivatives of the wavefield P(x, y, z; t) given by equation (H-18) and direct substitution into equation (H-19). The result is the dispersion relation

 ${\displaystyle k_{x}^{2}+k_{y}^{2}+k_{z}^{2}={\frac {\omega ^{2}}{v^{2}}}}$ (H-20)

of the scalar wave equation (H-19). Since equation (H-18) satisfies this relation, it is a valid solution to the scalar wave equation. In equation (H-19), v(x, y, z) is the propagation velocity of the compressional plane wave P(x, y, z; t).

We shall rewrite the phase term in the plane-wave solution given by equation (H-18) as

 ${\displaystyle P(x,y,z;t)=P_{0}\ \exp \left\{-i\omega \left[t-\left({\frac {k_{x}}{\omega }}x+{\frac {k_{y}}{\omega }}y+{\frac {k_{z}}{\omega }}z\right)\right]\right\},}$ (H-21)

so that we can define a 3-D traveltime surface T(x, y, z) as

 ${\displaystyle T(x,y,z)={\frac {k_{x}}{\omega }}x+{\frac {k_{y}}{\omega }}y+{\frac {k_{z}}{\omega }}z.}$ (H-22)

Substitute this definition into equation (H-21) to obtain the expression for the plane wave in terms of the traveltime surface T(x, y, z)

 ${\displaystyle P(x,y,z;t)=P_{0}\ \exp {\Big \{}-i\omega {\Big [}t-T(x,y,z){\Big ]}{\Big \}}.}$ (H-23)

Now, we need to verify that this form of the plane wave solution satisfies the scalar wave equation (H-19). Compute the partial derivative of the wave function P(x, y, z; t) of equation (H-23) with respect to the space variable x

 ${\displaystyle {\frac {\partial P}{\partial x}}=P_{0}(i\omega )\left({\frac {\partial T}{\partial x}}\right)\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}}$ (H-24)

and differentiate the result once more to get

 ${\displaystyle {\frac {\partial ^{2}P}{\partial x^{2}}}=-P_{0}\left[\omega ^{2}\left({\frac {\partial T}{\partial x}}\right)^{2}-i\omega {\frac {\partial ^{2}T}{\partial x^{2}}}\right]\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.}$ (H-25a)

Repeat the differentiation for the space variable y

 ${\displaystyle {\frac {\partial ^{2}P}{\partial y^{2}}}=-P_{0}\left[\omega ^{2}\left({\frac {\partial T}{\partial y}}\right)^{2}-i\omega {\frac {\partial ^{2}T}{\partial y^{2}}}\right]\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}}$ (H-25b)

and for the space variable z

 ${\displaystyle {\frac {\partial ^{2}P}{\partial z^{2}}}=-P_{0}\left[\omega ^{2}\left({\frac {\partial T}{\partial z}}\right)^{2}-i\omega {\frac {\partial ^{2}T}{\partial z^{2}}}\right]\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.}$ (H-25c)

Next compute the partial derivative of the wave function P(x, y, z; t) of equation (H-23) with respect to the time variable t

 ${\displaystyle {\frac {\partial ^{2}P}{\partial t^{2}}}=-P_{0}\omega ^{2}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.}$ (H-25d)

Substitute equations (H-25a,H-25b,H-25c,H-25d) into equation (H-19) and combine the terms into real and imaginary parts

 ${\displaystyle \omega ^{2}\left[\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}\right]-i\omega \left({\frac {\partial ^{2}T}{\partial x^{2}}}+{\frac {\partial ^{2}T}{\partial y^{2}}}+{\frac {\partial ^{2}T}{\partial z^{2}}}\right)={\frac {\omega ^{2}}{v^{2}(x,y,z)}}.}$ (H-26)

Since the term on the right-hand side is real, the imaginary part of the left-hand side has to vanish, leading to the final expression

 ${\displaystyle \left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}={\frac {1}{v^{2}(x,y,z)}},}$ (H-27)

which is called the eikonal equation. It gives the traveltime value T(x, y, z) for a ray passing through a point (x, y, z) in a medium with velocity v(x, y, z).

A solution to the eikonal equation, T(x, y, z) = constant represents the wavefront at an instant of time. Kinematically, a solution to the eikonal equation (H-27) should also be a solution to the wave equation (H-19). This raises the important question as to under what circumstances the eikonal equation may be considered a good approximation to the wave equation. To investigate this matter, we shall consider a plane wave function as in equation (H-23) but with a spatially varying amplitude P0(x, y, z) [17]:

 ${\displaystyle P(x,y,z;t)=P_{0}(x,y,z)\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}}$ (H-28)

and test its validity as a solution to the wave equation (H-19).

Compute the partial derivative of the wave function P(x, y, z; t) of equation (H-28) with respect to the space variable x

 ${\displaystyle {\frac {\partial P}{\partial x}}=\left({\frac {\partial P_{0}}{\partial x}}+i\omega P_{0}{\frac {\partial T}{\partial x}}\right)\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}}$ (H-29)

and differentiate the result once more, then simplify to get

 ${\displaystyle {\frac {\partial ^{2}P}{\partial x^{2}}}=\left\{\left[{\frac {\partial ^{2}P_{0}}{\partial x^{2}}}-\omega ^{2}P_{0}\left({\frac {\partial T}{\partial x}}\right)^{2}\right]+i\omega \left[2{\frac {\partial P_{0}}{\partial x}}{\frac {\partial T}{\partial x}}+P_{0}{\frac {\partial ^{2}T}{\partial x^{2}}}\right]\right\}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.}$ (H-30a)

Repeat the differentiation for the space variables y

 ${\displaystyle {\frac {\partial ^{2}P}{\partial y^{2}}}=\left\{\left[{\frac {\partial ^{2}P_{0}}{\partial y^{2}}}-\omega ^{2}P_{0}\left({\frac {\partial T}{\partial y}}\right)^{2}\right]+i\omega \left[2{\frac {\partial P_{0}}{\partial y}}{\frac {\partial T}{\partial y}}+P_{0}{\frac {\partial ^{2}T}{\partial y^{2}}}\right]\right\}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}}$ (H-30b)

and z

 ${\displaystyle {\frac {\partial ^{2}P}{\partial z^{2}}}=\left\{\left[{\frac {\partial ^{2}P_{0}}{\partial z^{2}}}-\omega ^{2}P_{0}\left({\frac {\partial T}{\partial z}}\right)^{2}\right]+i\omega \left[2{\frac {\partial P_{0}}{\partial z}}{\frac {\partial T}{\partial z}}+P_{0}{\frac {\partial ^{2}T}{\partial z^{2}}}\right]\right\}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.}$ (H-30c)

Next compute the partial derivative of the wave function P(x, y, z; t) of equation (H-28) with respect to the time variable t

 ${\displaystyle {\frac {\partial ^{2}P}{\partial t^{2}}}=-P_{0}\omega ^{2}\ \ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}.}$ (H-30d)

Substitute equations (H-30a,H-30b,H-30c,H-30d) into equation (H-19), and combine the terms into real and imaginary parts to obtain

 ${\displaystyle {\begin{array}{l}\left\{\omega ^{2}P_{0}\left[\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}\right]-\left({\frac {\partial ^{2}P_{0}}{\partial x^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial y^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial z^{2}}}\right)\right\}\\\quad -i\omega \left[2\left({\frac {\partial P_{0}}{\partial x}}{\frac {\partial T}{\partial x}}+{\frac {\partial P_{0}}{\partial y}}{\frac {\partial T}{\partial y}}+{\frac {\partial P_{0}}{\partial z}}{\frac {\partial T}{\partial z}}\right)+P_{0}\left({\frac {\partial ^{2}T}{\partial x^{2}}}+{\frac {\partial ^{2}T}{\partial y^{2}}}+{\frac {\partial ^{2}T}{\partial z^{2}}}\right)\right]={\frac {\omega ^{2}}{v^{2}}}P_{0}.\end{array}}}$ (H-31)

Equate the real parts of both sides of equation (H-31) to get

 ${\displaystyle \left[\left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial y}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}\right]-{\frac {1}{\omega ^{2}P_{0}}}\left({\frac {\partial ^{2}P_{0}}{\partial x^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial y^{2}}}+{\frac {\partial ^{2}P_{0}}{\partial z^{2}}}\right)={\frac {1}{v^{2}}},}$ (H-32a)

and equate the imaginary parts to get

 ${\displaystyle {\frac {2}{P_{0}}}\left({\frac {\partial P_{0}}{\partial x}}{\frac {\partial T}{\partial x}}+{\frac {\partial P_{0}}{\partial y}}{\frac {\partial T}{\partial y}}+{\frac {\partial P_{0}}{\partial z}}{\frac {\partial T}{\partial z}}\right)+\left({\frac {\partial ^{2}T}{\partial x^{2}}}+{\frac {\partial ^{2}T}{\partial y^{2}}}+{\frac {\partial ^{2}T}{\partial z^{2}}}\right)=0.}$ (H-32b)

Now we analyze the implications of equations (H-32a,H-32b). To reduce equation (H-32a) to the eikonal equation (H-27), it follows that the second term on the left-hand side of equation (H-32a) must vanish. This is only possible if we let 1/ω go to zero, or in other words, if we make the high-frequency assumption. Thus, for a wave function, such as that given by equation (H-23), with spatially varying amplitudes, the eikonal equation is a good approximation to the wave equation at the high-frequency limit. Additionally, for T(x, y, z) to be a good approximation to the solution of the wave equation, the first term on the left-hand side of equation (H-32b) must be negligible compared to the second term. In fact, the wave function given by equation (H-23) is a special form of the generalized asymptotic ray series solution given by [41]

 ${\displaystyle P(x,y,z;t)=\sum _{n=0}^{\infty }P_{n}\ \exp \left\{-i\omega \left[t-T(x,y,z)\right]\right\}(i\omega )^{-n},}$ (H-33)

where n = 0, 1, 2, …. Note that the terms for n ≥ 1 in equation (H-33) become insignificant for a high-frequency source. Specifically, in the high-frequency limit, the summation given by equation (H-33) reduces to the wave function of equation (H-23).

Since λ = 2πv/ω, where λ is the wavelength, the high-frequency limit is equivalent to small wavelengths. How small the wavelength should be for a solution of the eikonal equation to be a good approximation to the wave equation is an important consideration in the practical validity of the eikonal equation itself. The approximation is valid if the fractional change in the velocity gradient Δv is much less than the frequency v/λ [17]. In practice, this approximation will not be met across a layer boundary with a sharp velocity contrast or in a layer with velocity variations that occur within a spatial extent much less than the wavelength. As a direct consequence of this approximation, the eikonal equation is a good approximation to the wave equation only for a medium in which velocities do not vary rapidly, and especially, do not exhibit sharp discontinuities.

### H.3 Finite-difference solution to the eikonal equation

Just as the scalar wave equation itself can be solved using a finite-difference scheme to extrapolate a wavefield from one depth level to another (migration principles), the eikonal equation can be solved using a finite-difference scheme to compute the traveltimes along wavefronts expanding from a source [42] [10] [24] [25] [21]. For simplicity, we shall rewrite the eikonal equation (H-27) for the 2-D case in (x, z) coordinates

 ${\displaystyle \left({\frac {\partial T}{\partial x}}\right)^{2}+\left({\frac {\partial T}{\partial z}}\right)^{2}={\frac {1}{v^{2}(x,z)}}.}$ (H-34)

Consider the 2-D computation mesh sketched in Figure H-2a. We want to compute the traveltime T of the eikonal equation (H-34) at grid point (x + Δx, z + Δz) using the known traveltimes at grid points (x, z) and (x + Δx, z).

Computing the traveltimes at depth z + Δz from those at depth z means extrapolating T in the z-direction. Rewrite equation (H-34) in the form of an extrapolation equation [10]

 ${\displaystyle {\frac {\partial T}{\partial z}}=\pm {\sqrt {{\frac {1}{v^{2}}}-\left({\frac {\partial T}{\partial x}}\right)^{2}}},}$ (H-35a)

which can be rearranged to take the form

 ${\displaystyle v{\frac {\partial T}{\partial z}}={\sqrt {1-\left(v{\frac {\partial T}{\partial x}}\right)^{2}}}.}$ (H-35b)

Note that in equation (H-35b) we also have dropped one of the two solutions and only considered that which yields an increase in traveltime as we march in depth.

Now replace the differential operators with difference operators

 ${\displaystyle {\frac {v(z+\Delta z)T(z+\Delta z)-v(z)T(z)}{\Delta z}}=\left[1-\left({\frac {v(x+\Delta x)T(x+\Delta x)-v(x)T(x)}{\Delta x}}\right)^{2}\right]^{1/2},}$ (H-36a)

and rewrite explicitly in terms of T(z + Δz)

 ${\displaystyle T(z+\Delta z)=T(z)+{\frac {\Delta z}{v(z+\Delta z)}}\left[1-\left({\frac {v(x+\Delta x)T(x+\Delta x)-v(x)T(x)}{\Delta x}}\right)^{2}\right]^{1/2},}$ (H-36b)

where we have assumed that v(z)/v(z + Δz) is approximately unity. Also, we have avoided in our algebra crowding the terms in equation (H-36b) with the dual variables (x, z). Now introduce the dual variables explicitly to get the final expression for a finite-difference solution to the 2-D eikonal equation

 ${\displaystyle T_{i+1}^{k+1}=T_{i+1}^{k}+{\frac {\Delta z}{v_{i+1}^{k+1}}}\left[1-\left({\frac {v_{i+1}^{k}T_{i+1}^{k}-v_{i}^{k}T_{i}^{k}}{\Delta x}}\right)^{2}\right]^{1/2},}$ (H-37)

where ${\displaystyle T_{i}^{k}}$ is the discrete form of T(x, z) as labeled in Figure H-2a.

Figure H-2  Finite-difference mesh used for (a) 2-D and (b) 3-D solution to the eikonal equation. See text for details.

Now consider the 3-D computation mesh sketched in Figure H-2b. We want to compute the traveltime T of the eikonal equation (H-27) at grid point (x + Δx, y, z + Δz) using the known traveltimes at grid points (x, y, z), (x + Δx, y, z) and (x + Δx, y + Δy, z). The 3-D equivalent of equation (H-35b) is

 ${\displaystyle v{\frac {\partial T}{\partial z}}={\sqrt {1-\left(v{\frac {\partial T}{\partial x}}\right)^{2}-\left(v{\frac {\partial T}{\partial y}}\right)^{2}}}.}$ (H-38)

Following the same algebra as for the 2-D case, we obtain the finite-difference solution to the 3-D eikonal equation

 ${\displaystyle {\begin{array}{l}T_{i+1,j}^{k+1}=T_{i+1,j}^{k}+{\frac {\Delta z}{v_{i+1,j}^{k+1}}}\left[1-\left({\frac {v_{i+1,j}^{k}T_{i+1,j}^{k}-v_{i,j}^{k}T_{i,j}^{k}}{\Delta x}}\right)^{2}-\left({\frac {v_{i+1,j+1}^{k}T_{i+1,j+1}^{k}-v_{i+1,j}^{k}T_{i+1,j}^{k}}{\Delta y}}\right)^{2}\right]^{1/2},\end{array}}}$ (H-39)

where ${\displaystyle T_{i,j}^{k}}$ is the discrete form of T(x, y, z) as labeled in Figure H-2b.

Note that in equations (H-36), (H-37), and (H-39), unlike the previous authors [25] [26], the indices for the velocity v and traveltime T are kept the same. This was done with the stability of the finite-difference solution of the eikonal equation in mind. Actually, Kjartannson [43] has applied the same notion to the finite-difference solution to the scalar wave equation.[44]

## References

1. Judson et al., 1980, Judson, D. R., Lin, J., Schultz, P.S., and Sherwood, J.W.C., 1980, Depth migration after stack: Geophysics, 45, 376–393.
2. Schultz and Sherwood, 1980, Schultz, P.S. and Sherwood, J.W.C., 1980, Depth migration before stack: Geophysics, 45, 361–375.
3. Hubral (1977), Hubral, P., 1977, Time migration — some ray-theoretical aspects: Geophys. Prosp., 25, 738–745.
4. Larner et al., 1981, Larner, K. L., Hatton, L., and Gibson, B., 1981, Depth migration of imaged time sections: Geophysics, 46, 734–750.
5. Yilmaz and Lucas, 1986, Yilmaz, O. and Lucas, D., 1986, Prestack layer replacement: Geophysics, 51, 1355–1369.
6. Berryhill, 1979, Berryhill, J.R., 1979, Wave-equation datuming: Geophysics, 44, 1329–1333.
7. Berryhill, 1984, Berryhill, J.R., 1984, Wave-equation datuming before stack: Presented at the 54th Ann. Internat. Mtg., Soc. Expl. Geophys.
8. Taner et al., 1982, Taner, M.T., Cook, E.E., and Neidell, N.S., 1982, Paleo-seismic and color acoustic impedance sections: applications of downward continuation in structural and stratigraphic context: Presented at the 52nd Ann. Internat. Mtg., Soc. Expl. Geophys.
9. Al-Yahya, 1989, Al-Yahya, K., 1989, Velocity analysis by iterative profile migration: Geophysics, 54, 718–729.
10. Reshef and Kosloff, 1986, Reshef, M. and Kosloff, D., 1986, Migration of common-shot gathers: Geophysics, 51, 324–331.
11. Lee and Zhang, 1992, Lee, W. B. and Zhang, L., 1992, Residual shot-profile migration: Geophysics, 57, 815–822.
12. Deregowski, 1990, Deregowski, 1990, Common-offset migrations and velocity analysis: First Break, 8, 225–234.
13. Cox and Wapenaar, 1992, Cox, H. L. H. and Wapenaar, C. P. A., 1992, Macromodel estimation by common-offset migration and by shot-record migration: J. Seis. Expl., 1, 29–37.
14. Biondi et al., 1998, Biondi, B., Fomel, S., and Chemingui, N., 1998, Azimuth moveout for 3-D prestack imaging: Geophysics, 63, 574–588.
15. Chemingui and Biondi, 1999, Chemingui, N. and Biondi, B., 1999, Data regularization to inversion to common offset (ICO): 69th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1398–1401.
16. Biondi and Palacharla, 1996, Biondi, B. and Palacharla, G., 1996, 3-D prestack migration of common-azimuth data: Geophysics, 61, 1822–1832.
17. Officer, 1958, Officer, C. B., 1958, Introduction to the theory of sound transmission: McGraw-Hill Book Co.
18. Keho and Beydoun, 1988, Keho, T. H. and Beydoun, W. B., 1988, Paraxial ray Kirchhoff migration: Geophysics, 53, 1540–1546.
19. Červený et al., 1982, Červený, V., Popov, M. M., and Pšenčík, I., 1982, Computation of wavefields in inhomogeneous media — Gaussian beam approach: Geophys. J. Roy. Astr. Soc., 70, 109–128.
20. Červený et al., 1984, Červený, V., Klimes, L., and Pšenčík, I., 1984, Paraxial ray approximation in the computation of seismic wavefields in inhomogeneous media: Geophys. J. Roy. Astr. Soc., 79, 89–104.
21. Vinje et al., 1993, Vinje, V., Iversen, E., and Gjoystdal, H., 1993, Traveltime and amplitude estimation using wavefront construction: Geophysics, 58, 1157–1166.
22. Lecomte, 1999, Lecomte, I., 1999, Local and controlled prestack depth migration in complex areas: Geophys. Prosp., 47, 799–818.
23. Nichols, 1996, Nichols, D. E., 1996, Maximum-energy traveltimes calculated in the seismic frequency band: Geophysics, 61, 253–263.
24. Vidale, 1988, Vidale, J., 1988, Finite-difference calculation of traveltimes: Bull. Seis. Soc. Am., 78, 2026–2076.
25. Podvin and Lecomte, 1991, Podvin, P. and Lecomte, I., 1991, Finite-difference computation of traveltimes for velocity-depth models with strong velocity contrast across layer boundaries — a massively parallel approach: Geophys. J. Int., 105, 271–284.
26. Reshef, 1991, Reshef, M., 1991, Prestack depth imaging of three-dimensional shot gathers: Geophysics, 56, 1158–1163.
27. Meshbey et al., 1993, Meshbey, V., Kosloff, D., Ragoza, Y., Meshbey, O., Egozy, U., and Cozzens, J., 1993, A method for computing traveltimes for an arbitrary velocity model: Presented at the 45th Ann. EAGE Mtg.
28. Vesnaver, 1996, Vesnaver, A., 1996, Ray tracing based on Fermat’s principle in irregular grids: Geophys. Prosp., 44, 741–760.
29. Moser, 1991, Moser, T. J., 1991, Shortest-path calculation of seismic rays: Geophysics, 56, 59–67.
30. Sethian and Popovici, 1999, Sethian, J. A. and Popovici, A. M., 1999, 3-D traveltime computation using the fast marching method: Geophysics, 64, 516–523.
31. Claerbout, 1985, Claerbout, J.F., 1985, Imaging the earth’s interior: Blackwell Scientific Publications.
32. Gray, 1992, Gray, S., 1992, Frequency-selective design of the Kirchhoff migration operator: Geophys. Prosp., 40, 565–571.
33. Lumley et al. (1994), Lumley, D. E., Claerbout, J. F., and Bevc, D., 1994, Antialiased Kirchhoff 3-D migration: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1282–1285.
34. Claerbout, 1992, Claerbout, J.F., 1992, Earth soundings analysis: Blackwell Scientific Publications.
35. Abma et al. (1999), Abma, R., Sun, J., and Bernitsas, N., 1999; Antialiasing methods in Kirchhoff migration: Geophysics, 64, 1783–1792.
36. Kosloff and Kessler, 1987, Kosloff, D. and Kessler, D., 1987, Accurate depth migration by a generalized phase-shift method: Geophysics, 52, 1074–1084.
37. O’Brien and Gray, 1996, O’Brien, M. and Gray, S., 1996, Can we image beneath salt?: The Leading Edge, 17–22.
38. Lass, 1950, Lass, H., 1950, Vector and tensor analysis: McGraw-Hill Book Co.
39. Coulson, 1965, Coulson, C. H., 1965, Waves: A mathematical account of the common type of wave motion: Oliver and Boyd.
40. Jackson, 1962, Jackson, J. D., 1962, Classical electrodynamics: John Wiley and Sons.
41. Hubral and Krey, 1980, Hubral, P. and Krey, T., 1980, Interval velocities from seismic reflection time measurements: Soc. Expl. Geophys.
42. Qin et al., 1992, Qin, F., Luo, Y., Olsen, K. B., Cai, W., and Schuster, G. T., 1992, Finite-difference solution of the eikonal equation along expanding wavefronts: Geophysics, 57, 478–487.
43. Kjartannson (1979), Kjartansson, E., 1979, Modeling and migration by the monochromatic 45-degree equation: Stanford Expl. Proj. Rep. No. 15.
44. Chemingui, N. and Biondi, B., 1999, data regularization by inversion to common offset (ICO): Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1398–1401.