# User:Ageary/Yilmaz

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

## 7.0 Introduction

Subsurface geological features of interest in hydrocarbon exploration are three dimensional in nature. Examples include salt diapirs, overthrust and folded belts, major unconformities, reefs, and deltaic sands. A two-dimensional (2-D) seismic section is a cross-section of a three-dimensional (3-D) seismic response. Despite the fact that a 2-D section contains signal from all directions, including out-of-plane of the profile, 2-D migration normally assumes that all the signal comes from the plane of the profile itself. Although out-of-plane reflections (sideswipes) usually are recognizable by the experienced seismic interpreter, the out-of-plane signal often causes 2-D migrated sections to mistie. These misties are caused by inadequate imaging of the subsurface resulting from the use of 2-D rather than 3-D migration. On the other hand, 3-D migration of 3-D data provides an adequate and detailed 3-D image of the subsurface, leading to a more reliable interpretation.

A typical marine 3-D survey is carried out by shooting closely spaced parallel lines (line shooting). A typical land or shallow water 3-D survey is done by laying out a number of receiver lines parallel to each other and placing the shotpoints in the perpendicular direction (swath shooting).

Figure 7.0-1  A 2-D CMP-stacked section from an area with a subsurface structure associated with overthrust tectonics.

In marine 3-D surveys, the shooting direction (boat track) is called the inline direction; for land 3-D surveys, the receiver cable is along the inline direction. The direction that is perpendicular to the inline direction in a 3-D survey is called the crossline direction. In contrast to 2-D surveys in which line spacing can be as much as 1 km, the line spacing in 3-D surveys can be as small as 25 m. This dense coverage requires an accurate knowledge of shot and receiver locations.

The size of the survey area is dictated by the areal extent of the subsurface target zone and the aperture size required for adequate imaging of that target zone. This imaging requirement means that the areal extent of a 3-D survey almost always is larger than the areal extent of the objective.

A few hundred thousand to a few hundred million traces normally are collected during a 3-D survey. In a modern marine 3-D seismic survey, typically, more than 100 000 traces per square km are recorded. Most 3-D surveys are aimed at detailed delineation of already discovered oil and gas fields. Additionally, 3-D surveys are repeated over the same area at appropriate intervals, say every few years, to monitor changes in fluid saturation which may be inferred from changes in seismic amplitudes. By mapping changes in fluid saturation, changes in fluid flow directions may also be inferred and used for planning of production wells. To make use of seismic amplitudes for reservoir monitoring, however, data from all vintages must be processed consistently using a processing sequence aimed at preserving relative amplitudes. Seismic monitoring of oil and gas reservoirs by using time-lapsed 3-D surveys has come to be called the 4-D seismic method (4-D seismic method).

The basic principles of 2-D seismic data processing still apply to 3-D processing. In 2-D seismic data processing, traces are collected into common-midpoint (CMP) gathers to create a CMP stack. In 3-D data processing, traces are collected into common-cell gathers (bins) to create common-cell stacks. A common-cell gather coincides with a CMP gather for swath shooting. Typical cell sizes are 25 × 25 m for land surveys and 12.5 × 25 m for marine surveys.

Conventional 3-D recording geometries often complicate the process of stacking the data in a common-cell gather. Cable feathering in marine 3-D surveys can result in traveltime deviations from a single hyperbolic moveout within a common-cell gather. For land 3-D surveys, azimuth-dependent moveout within a common-cell gather is an issue.

After stacking, the 3-D data volume is migrated. Before migration, the data sometimes need to be trace-interpolated along the crossline direction to avoid spatial aliasing. The migrated 3-D data volume then is available to the interpreter as vertical sections in both the inline and crossline directions and as horizontal sections (time slices). The interactive environment with powerful 3-D visualization tools provides an efficient means for interpretation of the sheer volume of 3-D migrated seismic data. Fault correlation, horizon tracking, horizon flattening, and some image processing techniques can be adapted to the interactive environment to help improve interpretation.

Figure 7.0-2  A 2-D CMP-stacked section (a) before, and (b) after migration.

### The need for imaging in three dimensions

Figure 7.0-1 shows a CMP-stacked section from an area with very complex subsurface geology. As we examine the subsurface structure inferred by this section starting from the left side, we see a section from top to bottom with good-quality reflections. Then, most likely because of surface conditions, there is a zone of virtually no reflection. The remaining part of the section gives the appearance of very tight anticlinal features cascaded on top of one another. Nevertheless, it is most likely that these features have a strong 3-D aspect. The important point is that even if we were able to supply the correct migration velocity field for the section in Figure 7.0-1, 2-D migration of this 2-D cross-section of a 3-D wavefield would not yield the correct subsurface image.

We have always recorded 3-D seismic wavefields since the beginnings of the seismic method. A 2-D seismic profile is merely a cross-section of a truly 3-D propagating wavefield. Until the mid-1970s, we have been recording cross-sections of 3-D wavefields and processing as if the wavefield is 2-D, rather than 3-D. Although we still continue conducting 2-D seismic surveys for regional exploration, since the mid-1980s, 3-D seismic surveys have been used for prospect generation and detailed delineation of oil and gas reservoirs.

3-D effects can be minimized on 2-D lines by choosing the dominant dip direction as the shooting direction. For instance, in areas with structures that have resulted from overthrust tectonics, 2-D lines are traversed perpendicular to thrust fronts. Nevertheless, you will certainly need strike lines to tie the dip lines.

If there are any 3-D effects on a 2-D line, they can often be recognized. Figure 7.0-2 shows a CMP-stacked section before and after migration. Follow the unconformity A on the migrated section — events appear geologically plausable until we encounter an interfering feature B at the steep flexure. This event is not immediately correlatable with other events on the section, and is almost certainly out-of-plane. Although coherent on its own, it simply does not belong to the plane of this section. Only after 3-D migration of a 3-D volume of stacked data, will this event be imaged at its correct subsurface location. As a result of 3-D migration, this sideswipe event will have moved out of the plane of the section in Figure 7.0-2.

The section shown in Figure 7.0-3 unambiguously demonstrates the need for 3-D migration of 3-D data to achieve a complete and accurate imaging of the subsurface. The stacked section (Figure 7.0-3a) indicates irregular topography associated with a sea-bottom canyon. At first impression, there should be no problem in imaging the water-bottom reflector correctly. All you need to do is Stolt migration using the constant velocity associated with the water layer. The result shown in Figure 7.0-3b exhibits significant undermigration of the event associated with the water bottom. In other words, we have not achieved a complete and accurate imaging even though we used the correct migration velocity.

If it is not an erroneously too-low velocity, then what is the cause of this undermigration? To seek the answer to this question, consider the earth model in Figure 7.0-4, which consists of a dipping plane interface in a homogeneous medium. Examine line A along the dip direction. If the survey consisted of lines parallel to the dip direction, then a 2-D assumption about the subsurface would be valid. No out-of-plane signal would be recorded and 2-D migration of these dip lines would be correct, as shown in Figure 7.0-5a. After 2-D migration, point D beneath surface point X is migrated updip to its true subsurface position D′. Now consider line B, which is in the strike direction in Figure 7.0-4. Line B ties line A at surface position X. The reflection from the subsurface point D′ in Figure 7.0-5a is recorded on both lines A and B at their intersection point X. The event on line B is an out-of-plane reflection, while the event on line A is not. However, the reflection on the strike line from the dipping interface shows no dip (Figure 7.0-5b). Since migration does not alter the position of flat events, the migrated section for the strike line is identical to the corresponding unmigrated section. While the unmigrated sections associated with the dip line A and the strike line B would tie, after migrating both sections, a mistie would occur (Figure 7.0-5b).

The subsurface generally dips in many directions in areas where structural traps of interest exist. Thus, it is not possible to identify the inline direction as either a dip or strike direction. This is the case for line C in Figure 7.0-4. Examine the positioning after migration of point D beneath intersection point X for all three lines in the plan view in Figure 7.0-6. Point D is moved to the true subsurface position D′ along dip line A. The same point does not move after migration along strike line B. It moves to D″ along line C. Note that the lateral displacement DD″ along the oblique line C is less than the lateral displacement DD′ along the dip line. This is because the apparent dip perceived along the oblique line is less than the true dip of the plane interface that is perceived along the dip line. For accurate positioning, a second pass of migration is suggested in the direction perpendicular to line C to move the already migrated energy from D″ to its true subsurface position D′.

Principles of 3-D migration are discussed in 3-D poststack migration. Here, we shall assess the interpretational differences between 2-D and 3-D migrations. Figure 7.0-7 shows an inline (left column) and a crossline (right column) stacked section from a land 3-D survey. Pretend that these two sections are from a 2-D survey and migrate them individually to get the 2-D migrated sections. Then, take the entire volume of 3-D stacked data and perform 3-D migration, and subsequently display the data corresponding to these lines. Note that 3-D migration has yielded a better definition of the top T of the salt dome and better delineation of the faults along the base B of the salt dome. There is no doubt that the interpretation based on 2-D imaging is significantly different from interpretation based on 3-D imaging. It is not really a question of which section has the most lateral continuity of events — what matters is the accuracy of imaging we achieve from 2-D versus 3-D migration.

Figure 7.0-8 shows another example of the significant improvement available from the interpretation of a 3-D migrated section. Note that the two salt domes and the syncline between are delineated better after 3-D migration. The zone between the two salt domes exhibits undermigrated events on the 2-D migrated section (Figure 7.0-8b) in the same manner as described in Figure 7.0-6.

Three-dimensional migration often produces surprisingly different sections from 2-D migrated sections. The example in Figure 7.0-9 shows a no-reflection zone on the 2-D migrated section, while the same zone contains a series of continuous reflections on the 3-D migrated section that are easily correlated with reflections outside that zone. The 2-D migrated section exhibits the shadow events associated with the sideswipe energy masking the continuous reflections in the middle portion of the section. When we do 2-D migration, we confine the movement of the energy into the plane of the line itself. So the energy contained in the unmigrated stacked section in Figure 7.0-9 is indeed the same as the energy contained in the 2-D migrated section — only that it has been moved somewhere else on the section. As a result of moving the energy during 3-D migration within the 3-D volume, some moved into the plane of the section as in Figure 7.0-9 from others and some moved out of the section and migrated into the others.

Figure 7.0-9  (a) Another inline stacked section from the same marine 3-D survey as in Figure 7.0-8, (b) 2-D migration, (c) 3-D migration. (Data courtesy Amoco Europe and West Africa, Inc.)

As noted earlier, 2-D migration can introduce misties between 2-D lines in the presence of dipping events. Two-dimensional migration cannot adequately image the subsurface, while 3-D migration eliminates these misties by completing the imaging process. From the field data examples, we see that 3-D migration provides complete imaging of the 3-D subsurface geology. In contrast, 2-D migration can yield inadequate results. The difference between 2-D seismic and 3-D seismic is the way in which migration is performed. Dense coverage on top of a target zone, say a 12.5-m inline trace spacing and a 25-m crossline trace spacing, will not necessarily provide adequate subsurface imaging unless migration is performed in a 3-D sense.

Cable feathering in marine 3-D surveys and swath shooting geometry in land 3-D surveys give rise to source-receiver azimuthal variations. As a direct consequence of this acquisition-related phenomenon, stacking velocities become not only dip dependent but also azimuth dependent. Three-dimensional dip-moveout correction fortunately accounts for both dip and azimuth effects on stacking velocities. Figure 7.0-10 shows an inline section from a 3-D poststack migrated data from a 3-D marine survey with 3-D DMO correction applied. This data set represents a case of events with conflicting dips — specifically, the steeply dipping fault plane reflections and the gently dipping reflections associated with the sedimentary strata. Note that 3-D DMO correction has preserved both of these events. This then enabled 3-D migration to produce a crisp image of the fault planes.

Figure 7.0-11 shows an inline section from a land 3-D survey. Note that the inline stacked section without 3-D DMO correction (Figure 7.0-11a) exhibits poor stacking along the steep flanks of the salt dome, particularly above 1 s. There is no apparent case of conflicting dips in this section along the salt flanks. What causes the poor stacking performance is the large source-receiver azimuthal variations that arise from the swath shooting geometry. After the application of 3-D DMO correction, the reflections along the steep flanks of the salt dome are now preserved on the stacked section (Figure 7.0-11b). Clearly, the quality of 3-D migration hinges upon the quality of the stacked volume of data. Compare the sections in Figures 7.0-11c and 7.0-11d, and note that 3-D migration of the data with 3-D DMO correction yields a better delineation of the salt dome geometry. The improvement resulting from 3-D DMO correction is primarily within the overburden above the salt dome. The subsalt region should not be affected fundamentally by 3-D DMO correction since this region is adversely influenced by the raypath distortions caused by the complex overburden associated with the top-salt boundary. As such, the imaging problem in the subsalt region is not within the realm of 3-D DMO correction and 3-D time migration; instead, it is a problem that needs to be handled by 3-D imaging in depth.

Figure 7.0-10  A cross-section from a 3-D poststack time-migrated volume of 3-D DMO-stacked data.

Recall from introduction to migration that different migration strategies are required for different types of subsurface geological circumstances (Table 4-1). Consider the case of conflicting dips with different stacking velocities associated with fault-plane reflections and the surrounding sedimentary strata. The marine data set shown in Figures 7.0-12 through 7.0-17 exhibits such a case. The section shown in Figure 7.0-12a represents the CMP stack along an inline traverse with no DMO correction. As such, fault-plane reflections have not been preserved and the 2-D poststack time migration of this stacked section yields a poor image of the fault planes (Figure 7.0-12b).

Since this data set is an inline extracted from a 3-D marine survey, traces in the CMP gathers have source-receiver azimuthal variations that result from cable feathering and multi-cable recording geometry. To correct for the source-receiver azimuthal effect on moveout and therefore stacking velocities, one needs to do 3-D, rather than 2-D, DMO correction (processing of 3-D seismic data). In contrast with the conventional stack (Figure 7.0-12a), note the presence of fault-plane reflections on the 3-D DMO stack shown in Figure 7.0-13a. If, however, the fault planes have a 3-D geometry, then 2-D poststack time migration still falls short of providing a good image of the subsurface (Figure 7.0-13b).

To correct for conflicting dips with different stacking velocities, strictly, one needs to do prestack, and not poststack, time migration (prestack time migration). However, prestack time migration, when performed in 2-D, once again falls short of meeting the accurate imaging requirements (Figure 7.0-14a). Because of the acoustic impedance changes along the fault planes caused by the juxtaposition of velocities across the fault blocks, it is useful to examine the migrated section in Figure 7.0-14a also with its polarity reversed (Figure 7.0-14b).

Imaging in 3-D alone would not provide a good image if the input stacked volume of data does not contain the fault-plane reflections (Figure 7.0-15). By combining 3-D DMO correction and 3-D poststack time migration (Figure 7.0-16), the imaging quality for the fault-plane reflections is greatly enhanced. The ultimate imaging solution in the time domain to the problem of conflicting dips with different stacking velocities is 3-D prestack time migration (Figure 7.0-17).

Figure 7.0-11  (a) A cross-section from a 3-D volume of CMP-stacked data associated with a land 3-D survey without 3-D DMO correction; (b) with 3-D DMO correction; (c) same cross-section as in (a) after 3-D poststack time migration; (d) same cross-section as in (b) after 3-D poststack time migration.

## 7.1 3-D survey design and acquisition

The ultimate goal of conventional processing of a 3-D survey data is to obtain a 3-D seismic image of the subsurface. In this chapter, as in migration, we shall expound primarily on time migration strategies. The image quality from migration depends on stack quality and accuracy in velocity estimation. However, two other factors control the fidelity of migration — aperture and spatial sampling. They also can dictate the design of the field survey.

### Migration aperture

Figure 4.1-1a shows a velocity-depth model that contains a dipping reflector segment CD buried in a homogeneous medium. Zero-offset modeling using normal-incidence rays yields the time section in Figure 4.1-1b. Although not shown in the figure, the time section also should include diffractions off the edges of the reflecting segment.

Migration moves event C′D′ on the time section to its true subsurface position CD, which is overlaid on the time section for comparison. The horizontal extent of the zone of interest is OA. If the line length were confined to OA during recording, then the time section would be blank. On the other hand, if the recording were confined to line segment AB, then event C′D′ would be absent from the migrated section. Although the target is confined to line segment OA, the time section must be recorded over a longer segment OB. The line length also must be long enough to include a significant part of the diffractions that may be present in the data. Additionally, the recording time must be long enough to include diffraction tails and all of the dipping events of interest. The spatial (in the horizontal direction) and the temporal (in the vertical direction) displacement of a point on a dipping event resulting from migration depends on medium velocity, depth, and dip of the event (Table 4-2). Thus, line length and line position on the surface must be chosen carefully based on the migration aperture needed to adequately image the target zone of interest.

These considerations apply just as well to 3-D surveys. Figure 7.1-1 is a depth contour map of a fictitious structural high. The subsurface extent of the objective portion of the structure is indicated by the smaller rectangle. Using the principles discussed above and in migration principles (based on Figure 4.1-14), the actual survey size needed to define the objective area is outlined by the larger rectangle.

Note that the survey area does not have to be extended equally in all directions. The northern flank of the structure is the steepest part. Therefore, the survey area must be extended most in that direction. Extensions in other directions are determined accordingly. Another consideration in extending the survey area is the required additional length in profile to achieve full-fold coverage over the already extended survey area. A typical subsurface anomaly with a lateral extent of, say, 4 × 4 km may require a 3-D survey over an area as large as 10 × 10 km.

### Spatial sampling

Spatial aliasing was discussed in detail in the 2-D Fourier transform and in relation to migration in further aspects of migration in practice. The spatial aliasing problem is caused by spatial undersampling of the wavefield to be migrated — for example, the stacked section. The spatial sampling of stacked data (without trace interpolation) is defined by the recording parameters. Therefore, receiver spacing, crossline spacing, and the crossline direction in relation to dominant dip direction used in the field must be chosen carefully.

From Figure 1.2-21, note that a relationship exists between the trace spacing on a stacked section, dip, and the frequency at which spatial aliasing begins to occur. Imagine normal-incidence rays recorded at two receivers, A and B. In the constant-velocity case, the angle between the surface and the wavefront is the true dip of the reflector from which these rays emerged. There is a time delay equivalent to travelpath CB between the receivers at A and B. If this time delay is half the period of a given frequency component of the signal arriving at the receivers, then that frequency is at the threshold of being aliased.

From the relationship given by equation (1-7), note that the maximum frequency that is not aliased gets smaller at increasingly steeper dips, lower velocities, and coarser trace spacing. From this relationship, an optimum trace spacing can be derived for the inline and crossline directions based on the knowledge of a regional velocity field and subsurface dips. Typical trace spacings in the inline and crossline directions in 3-D surveys are 12.5 to 25 m, and 25 to 50 m, respectively. Even if the trace spacing in the crossline direction is as small as possible, for economic reasons, it usually is greater than that in the inline direction. Because of this, trace interpolation may be required along the crossline direction before migrating the data.

Figure 7.1-2  A diagram showing the marine geometry for multicable and dual-source recording. See text for details. (Courtesy PGS Tensor.)

### Other considerations

Almost all of the field operational aspects of 2-D acquisition are applicable to 3-D surveys. For example, selection of navigation and recording equipment depends on field conditions. The operating environment also must be considered. In the marine environment, water depth, tides, currents, sea conditions, fishing and shipping activity, and obstacles such as drilling platforms, wrecks, reefs, and fish traps must be considered. Modern marine 3-D surveys are conducted by deploying up to 12 cables and multiple source arrays. The logistics of a multicable operation requires taking careful account of such circumstances and obstacles. On land, environmental restrictions, accessibility, topography, cultivation, and demographic restrictions are factors that can affect survey design and acquisition. Because of these restrictions, careful planning and some adjustment of nominal shooting geometry often are required to achieve acceptable fold and offset distribution. Accurate surveying is a necessity for 3-D surveys, since data are collected with such dense spatial sampling; statics resulting from line-to-line surveying errors can seriously degrade the image quality obtained from 3-D migration. In fact, some claim that positioning error, rather than economics, is the limiting factor for line spacing in marine 3-D surveys.

### Marine acquisition geometry

Old 3-D surveys were conducted using a single receiver cable and a single source array. This is known as line shooting. A modern marine 3-D survey involves shooting a number of closely spaced parallel 2-D subsurface lines. This is achieved by using multicables and multi-source arrays. While some surveys are known to have been recorded using 12 cables, the more common configuration is 4-8 cables with dual source arrays. The recording geometry for multicable marine surveys is similar to that of the multireceiver line recording geometry — known as swath shooting, used in land, shallow-water, and transition-zone surveys. Figure 7.1-2 shows marine recording geometry that involves 12 cables and dual source array. Since each source-cable combination yields midpoint locations along one subsurface line, this recording geometry yields 24 subsurface lines, simultaneously. As a result, multicable recording increases the productivity in acquisition by greatly reducing the time in the field. Nevertheless, issues with multicable recording arise in relation to large variations of source-receiver azimuths in relation to velocity estimation, migration and amplitude variation with offset analysis.

### Cable feathering

In reality, receiver cables are not straight and parallel to one another as illustrated in Figure 7.1-2. Instead, receiver cables are subject to a certain amount of sideways drift, called feathering, from the ideal cable lines. Cable feathering is caused by cross currents. The cable shapes associated with a selected set of shot points from a marine 3-D survey are shown in Figure 7.1-3. The angle between the actual cable position and the shot-line direction (the boat track) is called the feathering angle. This angle is not always constant, even along the same cable associated with a single shot.

Figure 7.1-4 shows two types of source-cable combination — a single-source and dual-cable, and dual-source and dual-cable. As a result of cable feathering, midpoints associated with each of the source-cable combination depart from a straight subsurface line as illustrated in Figure 7.1-2. Instead, midpoints are scattered within the hatched areas. For a typical 10-degree feathering angle and a 2400-m-long cable, the midpoint associated with the far receiver is offset more than 200 m off the shot line. This is four lines off at a 50-m line spacing. Note that the direction and strength of cross-currents determines the crossline width of the midpoint distribution. In the ideal case of zero feathering, midpoint distribution would have zero crossline width. The single-source, dual-cable configuration simultaneously yields two bin lines, whereas the dual-source, dual-cable configuration yields four bin lines. The latter is conducted by alternate pops produced by the two sources.

Because of significant variations in cable shape during recording, in reality, midpoint distribution in the crossline direction can be quite irregular. We must know exactly where each of the receivers is located along the cable, as well as the shot-point location. Navigation data collected on the survey vessel normally include the boat location, source location, and cable compass readings. There are 8 to 12 digital compasses along a typical marine cable. Readings from these devices allow computation of the (x, y) coordinates of the cable compasses. Cable shape then is computed based on a curve-fitting procedure that rejects any anomalous measurements. Navigation data are analyzed during processing, and quality control is carried out to derive the final shot-receiver locations.

### 3-D binning

As described in basic data processing sequence, processing of 2-D seismic data is done in midpoint-offset coordinates. This requires, first, sorting the recorded data into common-midpoint gathers. Similarly, processing of 3-D seismic data requires binning the recorded data into common-cell gathers. Albeit we shall devote the next section to processing of 3-D data, it is imperative to discuss 3-D binning here, for it is fundamentally related to the acquisition geometry.

To perform 3-D binning, first, a grid is superimposed on the survey area as illustrated in Figure 7.1-5. This grid consists of cells with dimensions of half the receiver group spacing in the inline direction, equivalent to the CMP spacing in 2-D processing, and the line spacing in the crossline direction.

By precisely determining the shot and receiver coordinates, we can determine the midpoint locations and gather them into bins (or cells) for stacking and migration. In reality, midpoint distribution within a cell is not necessarily uniform since cable shape varies from shot to shot and line to line (Figure 7.1-5). Midpoints may be clustered at one part of the cell; that is, the centroid of the midpoints is not necessarily at the center of the cell. Midpoint distribution also can vary from cell to cell. Some cells may contain more traces than others, while some cells may have less uniformly distributed midpoints than others.

Figure 7.1-5 shows the actual midpoint locations over a portion of a 3-D survey area generated by an 8-cable configuration. The grid with dimensions 12.5 × 25 m represents the bins, and the traces at the midpoint locations that fall within each bin constitute a common-cell gather. The number of the midpoints within each bin defines the fold for that bin. Note that the fold and midpoint distribution vary significantly from one bin to another.

Figure 7.1-6a shows the fold of coverage map of a marine 3-D survey. While the average fold is 30, note that the fold of coverage is not uniform over the survey area. Fold of coverage varies from one range of offsets to another. Figures 7.1-7a, 7.1-8a and 7.1-9a show the coverage maps for the near-offset, mid-offset, and far-offset ranges. Irregularities in recording geometry — specifically, nonuniform fold of coverage, irregular offset distribution per bin, and irregular midpoint distribution within each bin, can cause problems in processing of the 3-D data. Lack of uniformity in the fold of coverage causes inconsistency in the accuracy of velocity estimation from one analysis location to another and variations in the stacked amplitudes. Processes such as 3-D dip moveout (DMO) correction, 3-D prestack time and depth migration, and amplitude variation with offset analysis are adversely affected by the acquisition footprint.

Figure 7.1-3  Actual cable shapes for a marine recording using six cables and a single source array. There are eight shot locations in this diagram. (Courtesy PGS.)

Low-fold areas in a 3-D survey must be filled in during acquisition by shooting more lines at appropriate locations. If these deficiencies are discovered in processing, it is costly to send the vessel back for further acquisition. Quality control work therefore must be carried out on board to monitor the fold of coverage.

Consider some modifications to binning of 3-D data. A slight translation and rotation of the grid imposed on the survey area sometimes significantly reduces problems associated with binning. This type of grid optimization may yield a more uniform midpoint distribution within each cell and even improve the uniformity of the fold of coverage over the survey area.

It is common not to restrict the gridding to cells of equal size. More midpoints may be included in the cell from the neighboring cells by expanding the cell size in the crossline direction as needed, typically up to 50%, to achieve a uniform fold of coverage. Although the same midpoint could be used in more than one cell, restrictions often are imposed on the offset range and multiplicity to prevent an excessive use of each midpoint. Adjustment to cell size to achieve a uniform fold of coverage and offset distribution is called flexible binning.

Figure 7.1-6b shows the fold of coverage map as in in Figure 7.1-6a after flexible binning. Except for a few locations, note that the fold of coverage is fairly uniform over the survey area with a nominal value of 30. Following the flexible binning, the fold of coverage maps for the near-offset, mid-offset, and far-offset ranges are shown in Figures 7.1-7b, 7.1-8b, and 7.1-9b. Although, there still exist some variations in the fold over the survey area within the specified ranges of offsets, flexible binning has helped to remove the anomalous fold distribution present in the maps shown in Figures 7.1-7a, 7.1-8a and 7.1-9a.

### Crossline smearing

The common-cell sorting problem is not solved completely by restoring uniformity in the fold of coverage. As mentioned earlier, the centroid of the midpoints may not coincide with the center of the cell. Theoretically, if the centroid departs significantly from the center, then placing the stacked trace at the centroid, rather than at the cell center, may be considered. This destroys equal spacing of the stacked traces, primarily in the crossline direction. In principle, however, 3-D poststack migration based on the Kirchhoff integral method may be used to produce migrated data volume with uniform trace distribution.

Figure 7.1-4  Sketches of (a) a dual-cable (RC1, RC2) and single-source (S) marine recording geometry, (b) a dual-cable (RC1, RC2) and dual-source (S1, S2) marine recording geometry. Crosscurrents cause the marine cable to feather and drift sideways from the shot-line direction. This causes the midpoints to spread in the crossline direction over subsurface strips as denoted by the hatched areas. When data are sorted into common-cell gathers, each cell contains midpoints associated with more than one source line. Circles with arrows on the cable trajectories represent compasses used for cable locationing. (Courtesy Western Geophysical).

Before we embark on searching for a way to alleviate the problems associated with 3-D binning as a result of cable feathering as described above, first, we shall investigate the effect of cable feathering itself on stacking the data. Figure 7.1-10 shows a single cell that is associated with the field geometry sketched in the same figure. The cell is 12.5 m in the inline (labeled as shot line) direction and 50 m in the crossline direction. Different symbols represent the midpoints that are associated with different shot lines. This cell contains midpoints from six different shot lines. In Figure 7.1-10, the midpoint distribution corresponds to the ideal case of a constant feathering angle of 10 degrees over the entire survey. The shooting direction is assumed to be the same for all the shot lines that contribute midpoints to this cell. Assume a simple earth model with a single, 30-degree dipping event. Also assume a constant-velocity medium above the dipping interface.

With common-cell sorting, the traveltimes of the arrivals in a single cell may not follow a single hyperbolic moveout curve [1], [2]. Figure 7.1-11 shows the traveltime curves for the cell under consideration (Figure 7.1-10) for three different shooting directions:

1. Strike-line shooting — no dip perceived along the inline direction (maximum cross-dip case).
2. Shooting direction at a 45-degree azimuth with respect to dip direction.
3. Dip-line shooting — no dip perceived along the crossline direction.

Note that data from different shot lines contribute to different portions of the traveltime curves. For a single dipping event, the 2-D recording geometry normally yields a hyperbolic moveout curve. As the line orientation becomes more parallel to the reflector strike, cross-dip increases and traveltimes deviate from the ideal hyperbolic moveout curve. This ideal hyperbola corresponds to the case of no cable feathering, in which all midpoints in the cell coincide with the cell center. The traveltime deviation is worse when shooting in the strike direction (Figure 7.1-11a). The traveltime deviation increases with increasing feathering angle, cross-dip, and cell dimension in the crossline direction. It is also more significant in the case of large moveout that occurs at lower velocities and shallow depths. If a common-cell stack were done along the best-fit hyperbolic path, then a loss of high frequencies would be expected.

How significant is the high-cut action caused by midpoint scatter in the crossline direction? Measure the time differences between the ideal hyperbolic path and the actual moveout curve in Figure 7.1-11a to derive the stacking operator plotted with its amplitude spectrum in Figure 7.1-12. Note the high-cut filtering action of the stacking operator; the -6 dB amplitude level is at about 70 Hz. Whether or not the effects of midpoint scatter are significant on a particular data set depends on the amount of cross-dip, cable feathering, and desired bandwidth of the stacked data.

A summary of the cause-and-effect relationship for cable feathering is given in Figure 7.1-13. Crosscurrents cause cable feathering, which makes the midpoints scatter within a cell in the crossline direction. If the shooting direction is such that a dipping interface has a cross-dip component, then the traveltimes associated with that interface in a common-cell gather deviate from a single hyperbolic moveout curve. This causes amplitude smearing during stacking, which acts as a high-cut filter as illustrated in Figure 7.1-12b. The cutoff frequency primarily is a function of the amount of cross-dip, reflection time, and velocity. The lower the velocity, that means usually the shallower the event, the larger the cross-dip (maximum along the strike line direction), and the larger the cell size the more likelihood of crossline smearing.

The most effective way to avoid crossline smear is to have the cell dimension in the crossline direction sufficiently small, which means shooting with close line spacing. The multicable recording geometry used in modern marine 3-D surveys provides the means for achieving sufficiently small crossline spacing to minimize crossline smearing on stacked data.

### Strike versus dip shooting

It is interesting to raise the question of which direction is better in recording — dip line or strike line, since we just noted that crossline smearing is most severe when shooting in the strike direction. Both types have advantages and disadvantages.

Aspects of dip-line shooting are summarized below:

1. With dip-line shooting, you have a better spatial sampling in the direction that you need most.
2. Since there is no cross-dip, there is no crossline smearing, and you can afford coarser spacing between lines.
3. With dip-line shooting, you have the disadvantage of complex moveout on common-cell data that can cause problems in velocity analysis [4].
4. You also have the problem of inline smear because of the reflection point dispersal along the dipping reflector associated with a nonzero-offset recording (Section E.1). Nevertheless, DMO correction handles this problem.
Figure 7.1-11  Traveltimes and the least-squares-fit hyperbolic moveout curves associated with (a) strike-line shooting, (b) shooting in the 45-degree direction from the downdip azimuth, and (c) dip-line shooting. These traveltimes were derived from a single planar interface with a 30-degree dip in a constant-velocity medium. The feathering angle is 10 degrees and midpoint distribution in the cell is shown in Figure 7.1-10. Numbers along the moveout curves denote the lines that contribute midpoints to the cell under study [3].

Aspects of strike-line shooting are summarized below:

1. With strike-line shooting, you have better attenuation of coherent noise because of the moveout behavior of side scatterers.
2. Since there is no dip perceived on a strike line, in principle, you do not need DMO correction.
3. With strike-line shooting, however, you have the disadvantage of perceiving the largest cross-dip, hence you would have the problem of crossline smear to deal with.
4. More importantly, strike-line shooting would require closer line spacing to prevent spatial aliasing of steep dips in the orthogonal direction; otherwise trace interpolation becomes important.

Among these various factors, adequate spatial sampling of steep dips is a requirement that takes precedence over the others. This means that, given the choice between the two types of shooting, dip-line shooting must be employed almost always in data acquisition. Nevertheless, modern 3-D surveys are conducted using sufficiently small inline and crossline spacings which more than accommodate the requirements for adequate spatial sampling of steep dips, minimizing crossline smearing and accurate implementation of 3-D DMO correction for removing inline smearing.

### Land acquisition geometry

Land 3-D acquisition commonly is carried out by swath shooting in which receiver cables are laid out in parallel lines (inline direction) and shots are positioned in a perpendicular direction (crossline direction). Figure 7.1-14 shows swath shooting with 6 receiver cables, each with 80 receiver groups, 50-m apart. The receiver cables are spaced apart from left to right at 100, 200, 100, 200, and 100 m distances. Shooting usually is perpendicular to the swath starting from the far left and moving in and away to the right of the swath. As one shot line is completed, the receiver cables are "rolled along" the swath a number of stations, equivalent to the shot-line spacing, and shooting is repeated.

Figure 7.1-12  (a) The stacking operator associated with midpoint scatter along the crossline direction. This operator was derived from traveltime deviations from the best-fit hyperbolic moveout curve in Figure 7.2-3a. (b) The amplitude spectrum suggests that high frequencies are attenuated (causing smearing of stacked amplitudes) because of midpoint scatter [3].

The recording geometry in Figure 7.1-14 provides a 25 × 25-m bin size. Once one swath is completed, another swath parallel to it is recorded; and this procedure is repeated over the entire survey area. Although not exercised in some 3-D surveys, it is important to leave some receiver cables on the ground to ensure proper coupling of statics between swaths. A complete survey plan, including the 12 swaths of receivers and shot locations, is shown in Figure 7.1-15.

The swath shooting method yields a wide range of source-receiver azimuths, which can be a concern during velocity analysis (processing of 3-D seismic data). The source-receiver azimuth is the angle between a reference line, such as a receiver line or a dip line, and the line that passes through the source and receiver stations. So why use a recording geometry that introduces problems in analyzing data? The main advantage of swath shooting is that it is economical.

Figure 7.1-16 shows the spider diagram at each bin center that describes the source-receiver azimuthal variations induced by two different recording geometries — shots placed in the direction orthogonal to the receiver lines, and shots placed in the diagonal direction at 45 degrees to the receiver lines. Note that the first type of recording geometry yields uniform azimuthal distribution from one bin location to another. The more unusual second type of recording geometry yields variations in the azimuthal distribution from once bin location to another. On the other hand, the offset distribution is more uniform in the case of the diagonal shooting as illustrated by the histogram plots in Figure 7.1-17.

Fold of coverage maps for the recording geometries illustrated in Figure 7.1-16 are shown in Figure 7.1-18. Because of operating conditions, uniform coverage usually is not achievable over the entire survey area. In the next section, we shall follow through the entire processing sequence for a land 3-D survey data set. We shall confront the challenges of an irregular recording geometry right at the outset of the analysis.

## 7.2 Processing of 3-D seismic data

Almost all concepts of 2-D seismic data processing apply to 3-D data processing. Additional complications do arise in 3-D geometry quality control, statics, velocity analysis, and migration. Editing traces with high-level noise, geometric spreading correction, deconvolution and trace balancing, field statics applications (for land and shallow-water data) are done as for 2-D surveys. In conventional 2-D processing, traces are collected into common-midoint (CMP) gathers, while in 3-D processing, traces are collected into common-cell gathers. A common-cell gather coincides with a CMP gather for swath shooting when the lines are straight. Sorting into common-cell gathers introduces special problems. For a dipping reflector, there is the problem of azimuthal variations of the normal moveout (NMO) within the cell for most land data and for marine data with significant cable feathering. While in this section, important aspects of marine and land 3-D processing are discussed, 3-D poststack migration is devoted to 3-D migration.

Figure 7.2-1 is the base map for a land 3-D survey that covers a surface area of nearly 270 km2. The bin size is 25 × 25 m, and there are 520 inlines and 820 crosslines. The 3-D survey presented here comprises more than five million recorded traces, and more than 400,000 stacked traces. Table 7-1 provides a list of the survey parameters.

 Shot interval in m 25 Group interval in m 50 Receiver line spacing in m 500 Number of receiver lines per swath 6 Number of groups per receiver line 80 Bin size in m × m 25 × 25 Fold of coverage 12 Number of bins per km2 1600 Number of bins in the survey 426,400 Sampling interval in ms 4 Maximum time in ms 4000 Number of prestack traces 5,116,800 Data volume in gigabytes 20.5
Figure 7.1-15  Shot and receiver locations for an entire land 3-D survey. There are 12 swaths like the one in Figure 7.1-6. Coverage varies over the survey between 12-fold and 24-fold.

Note the irregularities in vibroseis source-array locations — this is mainly a result of cultural restrictions, such as the presence of farm houses and croplands. Figure 7.2-1 also shows a sketch of the recording geometry. There are 6 receiver lines at 500-m intervals, each comprising 80 groups at 50-m intervals. Shots are placed perpendicular to the swath and at the centers of the receiver lines.

Figure 7.2-2 shows selected shot records from the 3-D survey of Figure 7.2-1. Each shot record comprises subrecords each of which corresponds to one of the six receiver lines in the swath. Depending on the shot location relative to the receiver lines in the swath, arrival times of refraction and reflection events vary from one subrecord to another. In some records, ground-roll energy is present; nevertheless, shot records contain generally contain fairly coherent reflection events.

Figure 7.2-3 shows the source and receiver elevations, and the full-fold of coverage over the survey area. Elevations generally increase from southwest to north-east. Note the irregular pattern of the source-array locations and the more-or-less regular layout of the receiver lines in the east-west direction. The average fold of coverage is 12 over the survey area. Note however the variations in fold that follow a regular pattern of striations, and the low-fold area in the center.

Aside from the full-fold coverage map (Figure 7.2-3), it is necessary to examine the fold of coverage for a range of offsets. Figure 7.2-4 shows three coverage maps for offset ranges 0-1000, 1000-2000, and 2000-4000 m. The near-offset fold of coverage appears to be fairly uniform over the survey area. The mid-offset fold of coverage appears to be low and exhibits some variations. Finally, the far-offset fold of coverage map shows that there are bins with missing far-offset traces. Coverage maps are essential to quality control in processing and interpretation. Variations in fold of coverage obviously have undesirable effects on velocity estimation, multiple attenuation, noise attenuation and amplitude variation with offset (AVO) analysis. With missing far offsets, for instance, success in multiple attenuation is compromised.

Fold of coverage should also be examined within the context of source-receiver azimuthal variations. Figure 7.2-5 shows coverage maps for source-receiver azimuths that fall within three ranges — range 1: (340-40) and (160-220); range 2: (40-100) and (220-280); and range 3: (280-340) and (100-160) degrees measured from the north. As a direct consequence of the recording geometry (Figure 7.2-1), the near offsets are largely confined to azimuths that fall in range 1, the middle offsets are largely confined to azimuths that fall in range 2, and the far offsets are largely confined to azimuths that fall in range 3.

### 3-D refraction statics corrections

Figures 7.2-6a and 7.2-7a show the shot and receiver statics, respectively. Field statics were computed in the manner described in refraction statics corrections. Specifically, shots and receivers were lowered from the topographic surface to a shallow refractor using the weathering velocity, then moved up to a flat datum using the refractor velocity. The information required to apply the field statics corrections — the weathering velocity, refractor (bedrock) velocity, and depth to bedrock, were derived from the uphole surveys which were carried out in the field at close spatial intervals.

If field statics corrections are judged to be inadequate, then a near-surface model should be estimated from inversion of refracted arrivals (refraction statics corrections). This obviously is possible only if first-break picking for deriving a set of refracted arrival times can be achieved with confidence.

The surface-consistent refraction statics model discussed in refraction statics corrections does not restrict shot and receiver stations to a survey line — they can be situated anywhere, as in a 3-D survey (Figure 7.2-1). Hence, equation (3-52), which is rewritten below, can be used to model the intercept time anomalies at all shot and receiver stations over the entire 3-D survey area as

 ${\displaystyle t_{ij}^{\prime }=T_{j}+T_{i}+2s_{b}h_{ij}.}$ (1)

Here, ${\displaystyle t_{ij}^{\prime }}$ are the modeled first-break picks associated with a shallow refractor, Tj are the intercept time anomalies at shot stations, Ti are the intercept time anomalies at receiver stations, sb is the bedrock slowness, and 2hij the shot-receiver separation.

As in the case of 2-D data (refraction statics corrections), the intercept time anomalies at shot and receiver stations, and the bedrock slowness (inverse of the bedrock velocity) are estimated from the model equation (1). The procedure is based on a formulation (refraction statics corrections), where the near-surface model parameters are estimated such that the difference between the modeled times t′ and the actual times t picked from first breaks is minimum in the least-squares sense.

### Azimuth dependence of moveout velocities

Following the application of field statics, data are sorted into common-cell gathers. The swath shooting geometry used commonly in land 3-D surveys produces common-cell gathers in which all midpoints coincide with the cell center. However, this technique also can result in large variations in the source-receiver azimuths and, thus, can result in traveltime deviations similar to or greater than those caused by midpoint scattering in marine surveys.

Consider a land 3-D recording geometry that comprises a range of shot-receiver azimuths as in Figure 7.2-5. Much as in the case of cable feathering (Figure 7.1-11), reflections from a dipping interface do not align along a single hyperbolic moveout curve. Instead, source-receiver azimuthal variations give rise to travel time deviations from the ideal hyperbolic moveout. Using a single velocity for the NMO correction would yield a stacked trace with high frequencies attenuated, again as in the case of cable feathering (Figure 7.1-12).

Stack attenuation increases with the increasing shot-receiver azimuth range. If this attenuation is serious, correction for the shot-receiver azimuth before stacking is needed. At first, it appears that this may be achieved by grouping the traces in a CMP gather into various ranges of azimuths and applying different velocities for moveout correction of each group of traces. Nevertheless, as we shall discuss shortly, 3-D DMO process implicitly corrects for the effect of source-receiver azimuth on moveout velocities.

But first, because of its historical significance, we shall briefly review the azimuth-dependent velocity analysis. Levin [5] in his classic paper developed the equations for moveout velocity and traveltime associated with a dipping reflector and for a source-receiver pair along an arbitrary azimuthal direction measured from the dip line direction. The geometry of Levin’s equations is shown in Figure 7.2-8. In normal moveout, we discussed the normal moveout for the case of a 2-D dipping reflector and derived the respective equations in Section C.3.

The 3-D equivalent of equation (3-10) is

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {4h^{2}(1-\sin ^{2}\phi \cos ^{2}\theta )}{v^{2}}},}$ (2)

while the 3-D equivalent of equation (3-11) is

 ${\displaystyle v_{NMO}={\frac {v}{\sqrt {1-\sin ^{2}\phi \cos ^{2}\theta }}},}$ (3)

where θ is the azimuth angle between the structural dip direction and the direction of the profile line, ϕ is the dip angle (Figure 7.2-8), t is the two-way traveltime associated with the nonzero-offset raypath from the source location to the reflection point on the dipping interface back to the receiver location, t0 is the two-way zero-offset time t0 associated with the normal-incidence ray-path at the midpoint location, vNMO is the moveout velocity, and v is the velocity of the medium above the dipping reflector.

Figure 7.2-9 shows a plot of the ratio of vNMO/v based on equation (3) as a function of azimuth angle for a range of dip angles [5]. For a dip line, the azimuth is zero; and for a strike line, the azimuth is 90 degrees. The ratio vNMO/v is significant when the line is shot at or near the structural dip direction.

Equation (3) describes an ellipse in polar coordinates (Figure 7.2-10). The radial coordinate is the NMO velocity vNMO, while the polar angle is the azimuth θ. Orientation of the major axis of the velocity ellipse is in the true dip direction.

The velocity ellipse can be constructed based on moveout velocities measured in three different directions. Lehmann and Houba [6] discuss several practical aspects of this measurement. By subgrouping traces from a CMP gather into three different azimuths (α1, α2, α3), the stacking velocities (vs1, vs2, vs3) can be estimated along those azimuths. If the stacking velocities in three different directions are known, then the semi-major and semi-minor axes a and b, respectively, and the orientation (downdip azimuth angle, ε) of the velocity ellipse can be defined (Figure 7.2-10). Once the velocity ellipse is constructed for a dipping reflector at a CMP location, then the following parameters can, in principle, be determined:

1. The true dip of the reflector, ϕ = cos−1(b/a).
2. The downdip azimuth ε.
3. The stacking velocity at any azimuth, ${\displaystyle v_{s}=v/{\sqrt {1-\sin ^{2}\phi \cos ^{2}\theta }},}$ where v = b (medium velocity above the reflector), and θ is the azimuth angle measured from the major axis a.

This is called the three-parameter velocity analysis. Each trace in the CMP gather is moveout corrected using the velocity along the corresponding shot-receiver azimuth.

A resolution problem arises in practice for the three-parameter velocity analysis. Land 3-D data often are of low fold; so partitioning the data into azimuth ranges with lower folds of coverage may not yield good velocity estimates. Azimuth variations can be minimized on land swath recording by minimizing perpendicular offset of sources from the receivers. Azimuth variation also often is largest on the near offsets where total moveout is small. Thus, the effect of source-receiver azimuth on stacking velocity (equation 2) can be fairly small. Another issue is offset distribution — ideally, it is desirable to have a wide range of offsets for each common-cell gather. Otherwise small errors in moveout correction may cause trace-to-trace chatter on stacked sections. Anyway, we leave the three-parameter velocity analysis and move onto 3-D dip-moveout correction to account for source-receiver azimuth and dip effect on stacking velocities.

### 3-D dip-moveout correction

Equation (2) indicates that each trace in a common-cell gather would have to be moveout-corrected using a different moveout velocity (equation 3) which accounts for the source-receiver azimuth associated with that trace.

In dip-moveout correction and prestack migration, we learned that the dip-moveout (DMO) process corrects the moveout velocity of a dipping reflector for the true dip angle and yields a dip-independent moveout velocity at the unmigrated position of the dipping reflector. Rewrite equations (2) and (3) for the 2-D case by setting the azimuth angle θ = 0 to get

 ${\displaystyle t^{2}=t_{0}^{2}+{\frac {4h^{2}\cos \phi }{v^{2}}}}$ (4)

and

 ${\displaystyle v_{NMO}={\frac {v}{\cos \phi }}.}$ (5)

So, the 2-D DMO correction maps velocity vNMO to velocity v. A single velocity function at a given CMP location can then be used to stack events with conflicting dips that otherwise would have different stacking velocities.

Again, refer to equation (3) and this time define an apparent dip angle ϕ′:

 ${\displaystyle \sin \phi ^{\prime }=\sin \phi \sin \theta .}$ (6)

Using this definition, the NMO velocity given by equation (3) is rewritten as

 ${\displaystyle v_{NMO}={\frac {v}{\cos \phi ^{\prime }}}.}$ (7)

This equation is identical to the 2-D equation (5), except that equation (7) refers to apparent dip, while equation (5) refers to true dip angle.

The implication of equation (7) is that the 3-D DMO process corrects for the effect of apparent dip on stacking velocities. Note from equation (6) that the apparent dip accounts for the combined effect of the true dip as measured along the dip-line direction and the source-receiver azimuth. The apparent dip is equal to true dip when the source-receiver azimuth coincides with the direction of the dip line. We conclude that 3-D DMO not only corrects for the true dip but also for the source-receiver azimuth effects on moveout velocity.

Application of 3-D DMO correction followed by 2-D inverse DMO correction have the effect of correcting for stacking velocities only for the azimuth angle θ in equation (3). This is called azimuth-moveout correction (AMO). The resulting CMP gathers can be treated as if they belong to a 2-D seismic line with shots and receivers placed along a straight traverse.

In dip-moveout correction and prestack migration, we learned that the 2-D DMO process removes the relection-point smear Δ defined as (Section E.1)

 ${\displaystyle \Delta ={\frac {h^{2}}{vt_{0}}}\sin 2\phi .}$ (8)

Equation (8), by way of equation (5-6), is the same as equation (5-10) of dip-moveout correction and prestack migration with offset 2h. The 3-D equivalent of equation (8) is [5]

 ${\displaystyle \Delta ={\frac {2h^{2}}{vt_{0}}}\sin \phi \cos \theta {\sqrt {1-\sin ^{2}\phi \cos ^{2}\theta }},}$ (9)

where, again, θ is the azimuth angle between the structural dip direction and the direction of the profile line (Figure 7.2-8).

Old marine 3-D surveys used to be conducted with single-cable or dual-cable geometry which did not incur large source-receiver azimuthal variations. In such cases, source-receiver azimuths are confined to a feathering angle, generally less than 7 degrees. (For two opposite shooting directions, it is twice the feathering angle.) As a result, the 3-D aspect of the DMO process is not significant for old marine 3-D surveys. In contrast, modern marine 3-D surveys are conducted using multicable recording geometries that incur large source-receiver azimuthal variations. Therefore, the 3-D aspect of the DMO process becomes important. Similarly, land 3-D surveys are usually carried out using swath shooting. This type of recording also gives rise to large source-receiver azimuthal variations; the DMO process, therefore, has to be applied in a 3-D sense.

Clearly, if the true dip is zero, irrespective of the source-receiver azimuth, the DMO process does not affect the data. Additionally, note from equations (2) and (3) that regardless of the true dip, the DMO process does not affect the data along the 90-degree source-receiver azimuth (strike-line shooting). The largest effect of DMO is along the zero-degree azimuth (dip-line shooting). Also, recall from principles of dip-moveout correction that the larger the offset, the slower the velocity, the shallower the reflector, and the steeper the dip, the more the action of the DMO process.

The 2-D DMO operator (principles of dip-moveout correction and E.4) maps the moveout-corrected amplitudes on a common-offset section along truncated elliptical trajectories that narrow at increasingly late times [7]. The 3-D DMO operator maps the amplitudes on the plane of source-receiver azimuths over an ellipsoid of revolution (Figure 7.2-11) [8]. The fidelity of 3-D DMO correction, therefore, depends on source-receiver azimuthal coverage, in particular, how well the ellipsoid is defined.

The 3-D DMO process can be conceptualized as the extension of the 2-D DMO process (Figure 7.2-11a). Take a moveout-corrected trace from a common-cell gather, which for the swath shooting coincides with a common-midpoint gather, and map the amplitude at time sample A to neighboring cells that are coincident with the source-receiver azimuthal direction associated with that input trace, along the elliptical trajectory of the 2-D DMO operator. The DMO ellipse is described by equation (5-21) and the amplitude mapping is described by equation (5-24). The excursion from the center trace depends upon the source-receiver separation, and the NMO-corrected time at the center trace. The amplitude at an output cell location along the mapping trajectory is the amplitude at the center trace scaled by the excursion from the center trace and the source-receiver separation. Notice that mapping is done over the plane of NMO-corrected time versus the source-receiver azimuthal direction.

Continue the process for all the traces from the same common-cell gather and map the amplitudes in the same manner. The output of the DMO process at any one cell location is simply a superposition of the contributions from the mapping of all the input traces.

Consider the hypothetical recording geometry with a common-cell gather coincident with a common-midpoint gather. Also, consider the traces in the gather covering a 360-degree source-receiver azimuthal range, but having the same source-receiver separation. As a result of the DMO correction described above, the amplitude at some specific input time sample A of each input trace is mapped along an elliptical trajectory. When combined, the elliptical trajectories associated with all the traces constitute an ellipsoid of revolution as shown in Figure 7.2-11a. This ellipsoid represents the kinematic aspect of the impulse response of the 3-D DMO operator in a constant-velocity medium. As for the 2-D DMO correction, vertical velocity variations may be accounted for in the design of a 3-D DMO operator [9].

In reality, there never exists a situation that corresponds to the recording geometry of Figure 7.2-11a. Instead, 3-D recording geometries give rise to a nonuniform source-receiver azimuthal coverage with each trace in a common-cell gather associated with a different source-receiver separation as illustrated in Figure 7.2-11b. Again, the amplitude at some specific input time sample A of each input trace is mapped along an elliptical trajectory. When combined, however, the elliptical trajectories associated with all the traces constitute only a sparse representation of the surface of the ideal DMO ellipsoid.

Aside from nonuniform source-receiver azimuthal and offset coverage, in the case of marine 3-D geometry, the fact that a common-cell gather does not coincide with a common-midpoint gather causes nonuniform spatial sampling of the data output from 3-D DMO correction. Figure 7.2-11b shows the surface projections of the output samples represented by the solid circles along each of the DMO ellipses. A similar distribution of the output samples is also shown in a map view in Figure 7.2-12 for the general case of irregular spatial sampling caused by the arbitrary source-receiver azimuthal and offset coverage with midpoint scattering over the survey area. As illustrated by the source-receiver pair S1R1 at the top portion of the cell grid in Figure 7.2-12, if there were no cable feathering, the output samples (again, denoted by the solid circles) from DMO correction applied to the trace at midpoint location M1 could be placed at the centers of the cells along the receiver cable. The output samples from DMO correction applied to a trace with cable feathering, however, do not necessarily coincide with the cell centers along the source-receiver direction. This realistic situation is illustrated by the source-receiver pair S2R2 in the middle portion of the cell grid in Figure 7.2-12.

Now imagine the superposition of the distribution of the output samples from 3-D DMO correction applied to all input traces in a common-cell gather. Just as the midpoint distribution in a common-cell gather input to 3-D DMO correction can be nonuniform because of cable feathering (Figure 7.1-5), the distribution and population of the output samples that fall within a cell from 3-D DMO correction also are nonuniform. The resulting fold of coverage from 3-D DMO correction may exhibit significant variations over the survey area. Moreover, the fold of coverage is time-dependent since the parameters of DMO ellipse are time dependent. As for conventional CMP stacking, the amplitudes after DMO correction need to be compensated for the variations in fold of coverage in a time-dependent manner [10] [11]. Coverage and binning issues are reviewed by Goodway and Ragan [12] for land, and by Brink and Ronen [13] for marine 3-D surveys.

Note that, ideally, the output samples, denoted by the solid circles in Figure 7.2-12, from DMO correction should be placed equidistant away from the midpoint location M2 along the line that connects the source S2 and the receiver R2 locations. In practice, however, the output samples are assigned to the cell centers denoted by the open circles. The unfortunate consequence of this output sampling is the vulnerability of the DMO-corrected data to spatial aliasing.

A way to correct for the irregular spatial sampling and undersampling of the 3-D DMO ellipse is to spread the output amplitudes away from the source-receiver azimuthal direction. [14]. Specifically, as illustrated in Figure 7.2-12, amplitudes not only are mapped onto the cell centers, denoted by the open circles, nearest the line that joins the source S2 and receiver R2 locations, but also to the neighboring cell centers denoted by × using a set of scalars that are inversely proportional to the distance of the cell from the source-receiver line. For instance, the output amplitude at A is not only scaled and placed at the cell center B but is distributed also to cell centers denoted by × using the appropriate scaling factors.

### Inversion to zero offset

To compensate the DMO operator for very severe geometry irregularities, a least-squares inversion scheme for data modeling can be implemented [15]. The underlying theory is identical to that of the discrete Radon transform (the radon transform). Recall that the discrete Radon transform maps amplitudes on a CMP gather that follow a hyperbolic moveout trajectory onto a point in the velocity-stack gather. The CMP gather can be reconstructed by applying inverse moveout correction and summing over the velocity axis. Let this inverse transformation be represented by the operator L. The transpose of this operator LT represents moveout correction and summing over the offset axis. If only the transpose operator LT were used to construct the velocity-stack gather from the CMP gather, hyperbolas would not map onto points. Instead, the amplitudes on the velocity-stack gather would exhibit smearing along the velocity axis. This is caused by discrete sampling along the offset axis and finite cable length. The operator LT alone does not account for these effects. Instead, we must use its generalized linear inverse (LTL)−1LT, which is a representation of the discrete Radon transform (Section F.3). Reconstruction of the CMP gather from the velocity-stack gather is one example of data modeling. Data modeling using the velocity-stack gather computed by the operator LT does not restore the amplitudes of the original CMP gather. In contrast, data modeling using the operator (LTL)−1LT accounts for the effects of the CMP geometry and faithfully restores the amplitudes of the CMP gather. For a formal mathematical treatise of the discrete Radon transform, see Section F.3.

Figure 7.2-12  Spatial aliasing in 3-D DMO correction. See text for details.

Now consider a 2-D operator LT that corresponds to DMO correction and stacking. As a result, a time sample on a moveout-corrected trace of a common-offset section is mapped onto the DMO ellipse, then is summed onto the zero-offset section. The inverse transformation is represented by the operator L. This operator maps the amplitudes on the zero-offset section back to the common-offset section by applying inverse DMO correction to the zero-offset data. Such reconstruction of the common-offset data from the zero-offset data is another example of data modeling. Now, suppose that there are geometry irregularities associated with the common-offset data. The modeling operator L then would not faithfully restore the amplitudes along the DMO ellipse on the zero-offset section back onto a point on the common-offset section. For the modeling operator L to work properly, we need to account for the geometry irregularities during the application of the DMO operator LT. Specifically, instead of applying LT, we must use its generalized linear inverse (LTL)−1LT. The covariance matrix LTL has the footprint of the acquisition geometry which is corrected for by the operator (LTL)−1LT. If LT is the DMO operator, the generalized linear inverse (LTL)−1LT may be termed as the inversion-to-zero-offset (IZO) operator [16] [15] [17] [18] [19]. Albeit its theoretical elegance, cost of the application of the IZO operator can be formidable. Instead, we may be content with the more modest spatial dealiasing scheme that is based on spreading the amplitudes away from the source-receiver azimuthal trajectory of the DMO operator [14] as described above.

### Aspects of 3-D DMO correction — a summary

Because of irregular spatial sampling associated with 3-D recording geometries, the 3-D DMO process is best applied in the time-space domain using the integral method (principles of dip-moveout correction). A velocity analysis is then performed following the DMO correction to estimate dip- and azimuth-corrected velocities.

In principles of dip-moveout correction, we discussed the principles of DMO correction and studied its practical aspects using synthetic and field data. We now extend aspects of 2-D DMO correction outlined in principles of dip-moveout correction to 3-D DMO correction.

1. The process of 3-D dip-moveout corrects for the dip and source-receiver azimuth effects on stacking velocities.
2. Thus, it preserves conflicting dips with different stacking velocities during CMP stacking.
3. The 3-D DMO stack, therefore, is a closer representation of a 3-D zero-offset section as compared to a conventional CMP stack volume of data based on normal-moveout correction, only.
4. The 3-D DMO stack can then be migrated using a 3-D zero-offset migration algorithm with greater accuracy.
5. Conflicting dips with different stacking velocities give rise to multivalued velocity picks from velocity spectra. Velocity analysis of 3-D DMO-corrected data alleviates this problem and increases the accuracy of picking an unambiguous velocity function from a velocity spectrum.
6. Velocities estimated from 3-D DMO-corrected data are dip and azimuth independent; therefore, they are more suitable to derive a migration velocity field as compared to velocities estimated from data without 3-D DMO correction.
7. 3-D DMO correction actually is a process of partial migration before stack. Specifically, it maps normal-moveout-corrected data to normal-incidence reflection points in the subsurface. As a result, the midpoint is a variant under DMO correction.
8. As a direct consequence of aspect (g), 3-D DMO correction removes the reflection point dispersal associated with nonzero-offset recording in the presence of dipping reflectors.
9. By 3-D DMO correction, prestack data can be implicitly regularized into common-midpoint gathers. This then facilitates sorting of prestack data into a set of common-offset volumes, each of which can be considered a replica of a 3-D zero-offset wavefield.
10. Following 3-D DMO correction, prestack data can be migrated so as to create CMP gathers in their migrated position (next section). This then enables us to conduct velocity analysis to derive a migration velocity field with greater confidence.
11. Finally, the CMP gathers from prestack time migration of 3-D DMO-corrected data can be used for amplitude variation with offset analysis.

### Velocity analysis

Once the data are sorted into common-cell gathers, velocity analysis is performed. In this respect, there is no difference between 2-D and 3-D data processing. For 2-D data processing, a number of neighboring CMP gathers are included in the velocity analysis to increase the signal-to-noise ratio. Similarly, a number of common-cell gathers are included in the 3-D velocity analysis. This number can be, say, 5 in the inline and 5 in the crossline direction, for a total of 25 common-cell gathers. For the 3-D example, velocity analyses are performed at certain intervals (e.g., 0.5 km) along selected inlines that may be as much as 0.5 km apart. The results of velocity analyses at selected control points are used to derive the 3-D velocity field for all common-cell gathers over the entire survey. This is achieved by performing a 3-D interpolation of the velocity functions between the control points.

### 3-D residual statics corrections

The surface-consistent residual statics model that we discussed in residual statics corrections does not restrict shot and receiver stations to a survey line; they can be situated anywhere as in a 3-D survey. Hence, equation (3-25), which is re-written below, can be used to model the static shifts at all shot and receiver stations over the entire 3-D survey area as

 ${\displaystyle t_{ij}^{\prime }=s_{j}+r_{i}+G_{k}+4M_{k}h_{ij}^{2}.}$ (10)

Here, ${\displaystyle t_{ij}^{\prime }}$ are the modeled traveltime deviations associated with a moveout-corrected reflection event within a specified time gate, sj is the shot residual statics, ri is the receiver residual statics, Gk is the structure term at the kth midpoint location, where k = (i + j)/2, and ${\displaystyle 4M_{k}h_{ij}^{2}}$ is the residual parabolic moveout term, with 2hij being the shot-receiver separation. The idea is to decompose the reflection time deviations picked from moveout-corrected cell-gathers into surface-consistent shot and receiver residual statics shifts, as well as structure and residual moveout terms. The procedure is based on a formulation (residual statics corrections), where the residual statics terms in equation (10) are estimated such that the difference between the modeled times t′ and the actual times t picked from moveout-corrected cell-gathers is minimum in the least-squares sense.

As it is for the 2-D case, the critical step is in picking the traveltime deviations from the common-cell gathers. For 3-D seismic data, an initial pilot trace is built from the common-cell gather of a selected line where the signal-to-noise ratio is good. Pilot traces for other common-cell gathers are computed from both local input traces and neighboring pilot traces. Note that land recording geometry must accommodate some overlap between adjacent swaths to ensure statics coupling between the swaths.

A fundamental question arises regarding 3-D DMO correction in relation to residual statics corrections. Specifically, which one should follow the other? We note from Figure 7.2-8 and Levin’s equation (2) that reflection times need to be corrected for the source-receiver azimuth effect by way of the 3-D DMO process. If 3-D DMO correction were to follow residual statics estimation, traveltime variations caused by the source-receiver azimuth effect would influence the statics estimates. In contrast, if 3-D DMO were to precede residual statics estimation, the moveout-corrected times used in DMO correction would be contaminated by statics errors. Nevertheless, this problem is of a lesser degree than the problem caused by the effect of traveltimes that were not corrected for source-receiver azimuth on statics estimates. Therefore, it is advisable to apply 3-D DMO correction prior to residual statics estimation.

Figures 7.2-6b and 7.2-7b show the shot and receiver residual statics over the survey area. Figure 7.2-13 shows a crossline from the 3-D survey of Figure 7.2-1 with and without residual statics corrections. Note the improvement in the reflection continuity at times below 2 s.

Following 3-D DMO correction and residual statics corrections, a final velocity analysis is carried out using a grid of the size that may be as small as 500 × 500 m. Figure 7.2-14 shows velocity panels from selected analysis locations over the survey area. Velocity analysis for land data often requires more than just a velocity spectrum (Figure 7.2-14d). A set of guide velocity functions that differ by a certain percentage from one another is specified. These functions are used to apply moveout correction to a set of CMP gathers (along, say, the inline direction) on both sides of the analysis location. The center CMP gather (Figure 7.2-14a) after moveout correction using the guide functions constitute one part of the velocity panel (Figure 7.2-14b). The CMP stack of the selected gathers at the analysis location with the guide functions constitute another part of the panel (Figure 7.2-14c). The velocity spectrum associated with the center CMP gather also is included in the panel (Figure 7.2-14d). Combined analysis of the gather, stack, and spectrum enables reliable picking of a velocity function at the analysis location.

By combining the velocity functions picked at analysis locations over the survey area, a 3-D velocity field is created. Figure 7.2-15 shows selected time slices from the 3-D velocity field. A thorough check on the velocity field is essential for stacking and time migration. For the latter, DMO-stacking velocity field (Figure 7.2-15) usually is smoothed spatially so as to eliminate lateral velocity variations that are judged to be unacceptable for time migration.

Figure 7.2-16 shows selected crosslines from the volume of stacked data without and with 3-D DMO correction. For the location of these crosslines, refer to the base map in Figure 7.2-1. Note that the differences between the sections with and without 3-D DMO correction are seen to the right of inline location 300. Specifically, note the enhanced flanks of diffractions associated with the tight imbricate structure that conflict with the near-flat reflections on the DMO-stacked data.

### 3-D migration

By using a smooth version of the 3-D velocity field (Figure 7.2-15) derived from the 3-D DMO-corrected data, the 3-D DMO-stack volume is then time-migrated. Figure 7.2-17 shows the selected crosslines as in Figure 7.2-16 before and after 3-D poststack time migration (see 3-D poststack migration for 3-D poststack migration). Following 3-D migration, the folded, imbricate structure has been delineated. Note, however, to the left of inline location 200, the left flank appears to terminate abruptly along a spurious fault line that dips down from left to right (see, for instance, crossline X-200). Recall the discussion on migration and line length (further aspects of migration in practice) and examine closely this zone on crossline X-260. Where there is no image of the folded structure to the left of inline location 100, there exist the trails of the migration paths — the wavefronts. So, the incomplete image of the left flank of the folded structure (the no-event zone to the left of inline location 200) is because of missing data — insufficient areal extent of the survey to the left of the imbricate structure.

We now turn our attention to the overthrust fault to the right of the imbricate structure. Refer to the left of inline location 200 on crossline X-260. Again, we observe a no-event zone — but, this one arises for a different reason. The culmination along the overthrust fault has caused significant ray bending to the extent that no data were recorded at the surface. Alternatively, we may be dealing with turning rays along the overthrust fault, which were not accounted for in migrating the data.

Figure 7.2-18 shows the selected inlines before and after 3-D poststack time migration. Note on inline I-80, there is an absence of events below 2 s after 3-D migration. This precisely corresponds to the missing data zone on the crosslines after migration (Figure 7.2-15, see for instance crossline X-260, left of inline location 200). We also observe a contrasting situation — events appear on migrated data in zones where no event was present on the unmigrated data (see for instance inline I-200, the zone below 3 s). As a result of 3-D migration, events move in and out of cross-sections in the inline and crossline directions. Lack of events after 3-D migration should not be alarming — they just moved to another location outside the plane of the cross-section of migrated data under consideration.

Figure 7.2-19 shows selected time slices from the unmigrated 3-D DMO-stack volume of data. Compare with the time slices from the 3-D poststack time-migrated volume of data shown in Figure 7.2-20. Note, for instance from the 2200-ms time slice of the migrated data, the delineation of the thrusted fault zone. Also, recall from migration principles that anticlines appear bigger than their actual sizes on unmigrated data. As such, compare the 2200-ms time slices from the unmigrated data (Figure 7.2-19) and the migrated data (Figure 7.2-20) and note that the contours in the latter bounded by the fault trend are narrowed, indicating an anticlinal, imbricate structure.

### Trace interpolation

A typical 3-D survey has a trace spacing in the crossline direction that normally is coarser than the trace spacing in the inline direction, in some cases, by as much as four times. This coarse spacing can cause spatial aliasing in the crossline direction. The problem of spatial aliasing was discussed extensively in the 2-D Fourier transform and further aspects of migration in practice within the context of 2-D Fourier transform and 2-D migration, respectively.

Recall that the coarser the trace spacing and the steeper the event dip of interest, the lower the threshold frequency at which spatial aliasing begins to take effect. In this section, we shall review trace interpolation as a means to circumvent the adverse effect of spatial aliasing in the crossline direction in 3-D migration. Trace interpolation does not create data, it merely unwraps the f − k spectrum of the input data so that aliased frequencies are mapped to the correct quadrant in the f − k domain. Finally, data do not necessarily need to be trace-interpolated in the crossline direction down to the trace spacing of the inline direction. Instead, signal bandwidth and subsurface dip can be taken into consideration to compute the optimum trace spacing to avoid spatial aliasing (equation 1-7).

During the first decade of 3-D seismic exploration, various methods have been developed for trace interpolation [20] [21]. Presently, trace interpolation usually is done using one-step complex spatial prediction filters [22] [23] in the frequency-space domain [24]. Consider a 2-D CMP-stacked data P(x, t) to be trace interpolated such that the trace spacing of the interpolated data is half of the trace spacing of the original data. Assume that P(x, t) is made up of a discrete set of dipping linear events and that the waveform for a given event is invariant from trace to trace. Also assume that the data to be interpolated do not contain any random noise. Under these assumptions, Pk(ω), the discrete Fourier transform of P(x, t) in the time direction, may be represented by the following combination of amplitude and phase spectra in the frequency-space domain [24]:

 ${\displaystyle P_{k}(\omega )=\sum _{j=1}^{N}A_{j}(\omega )\exp(-i\omega k\Delta \tau _{j}),k=0,1,2,\ldots ,M,}$ (11)

where Aj(ω) is the amplitude spectrum of the wavelet associated with the jth dipping event, Δτj is the time shift along the jth dipping event from trace to trace, k is the trace index, M is the number of traces, and N is the number of dipping events.

A one-step complex prediction filter is designed as a trace interpolation operator for each frequency component in the frequency-space domain and applied to the input data associated with the frequency component with twice the frequency [24]. The following are the steps involved in trace interpolation using one-step prediction filters. Details of the design and application of prediction filters for trace interpolation are provided in Section G.5.

1. Start with a 3-D volume stacked data P(x, y, t) that is to be migrated, and assume that the data volume is adequately sampled in the inline direction, but needs to be interpolated in the crossline direction before migration. Apply Fourier transform in the time direction, P(x, y, ω).
2. Sort the data into crosslines.
3. Then sort each crossline complex matrix of data P(y, ω) into complex arrays P(y) for each frequency component ω.
4. Design a one-step prediction filter from the data array P(y) of frequency ω/2 and apply it to P(y) of frequency ω to obtain the interpolated array Q(y) (Section G.5).
5. Interlace the original data array P(y) with the interpolated data array Q(y) to obtain the output array R(y) with twice the number of elements as in the input array P(y).
6. Repeat steps (d) and (e) for all frequencies. Then, combine the complex arrays R(y) and sort them into their crossline complex matrix R(y, ω).
7. Apply inverse Fourier transform to obtain the interpolated crossline section R(y, t).
8. Repeat steps (b) through (g) for all the crossline sections.

Figure 7.2-21 shows a CMP-stacked section with trace spacing of 10 m and the same stack with every other trace dropped so as to make the trace spacing 20 m, 40 m, and 80 m. Note the increasingly less distinctive character of the steep flanks of the diffractions and steeply dipping reflections at coarser trace spacing. Migrations of the four stacked sections in Figure 7.2-21 are shown in Figure 7.2-22. Note that migration of the stacked section with sufficiently small trace spacing provides a clear image, albeit in time, of both gently dipping and steeply dipping reflectors. With coarser trace spacing, however, the image quality of the steeply dipping events at the upper half of the stacked section begins to deteriorate rapidly. Spatial aliasing at coarse trace spacing has caused the aliased frequency components of the dipping events to move in the wrong direction. As a result, the migrated section in Figure 7.2-22d with very coarse trace spacing (80 m) is very noisy, especially in the upper half, and void of coherent reflections that are present in the unmigrated stacked section (Figure 7.2-21d).

Figure 7.2-23a shows the stacked section in Figure 7.2-21b with trace spacing of 20 m after interpolation to trace spacing of 10 m. Compare this interpolated section with the stacked section in Figure 7.2-21a and subtract the two sections to obtain the difference section shown in Figure 7.2-24a. The difference section contains alternating traces with zero sample values since trace interpolation retains the original traces in the input section with 20-m trace spacing (Figure 7.2-21b) while creating extra traces in between the original traces to produce a section with 10-m trace spacing (Figure 7.2-23a). Note that the traces with live samples in the difference section carry negligibly small energy associated with steeply dipping events within the upper half of the section.

Figure 7.2-23b shows the stacked section in Figure 7.2-21c with a trace spacing of 40 m after interpolation to trace spacing of 20 m. Compare this interpolated section with the stacked section in Figure 7.2-21b and subtract the two sections to obtain the difference section shown in Figure 7.2-24b. Again, the difference section contains alternating traces with zero samples interlaced with traces that contain nonzero samples corresponding to errors in trace interpolation.

Finally, Figure 7.2-23c shows the stacked section in Figure 7.2-21d with trace spacing of 80 m after interpolation to a trace spacing of 40 m. Compare this interpolated section with the stacked section in Figure 7.2-21c and subtract the two sections to obtain the difference section shown in Figure 7.2-24c. This difference section indicates that interpolation from a coarse trace spacing to a finer trace spacing may not be sufficiently accurate, and as a result, may not reproduce the signal coherency that is present in the original stacked section (Figure 7.2-21c). Hence, there is a limit as to how coarse the input trace spacing can be for trace interpolation to provide acceptable results. The degree of accuracy in trace interpolation has a direct impact on the fidelity of the image obtained from the migration of the interploated data. Compare the migrated sections in Figures 7.2-25a,b,c with those in Figures 7.2-22a,b,c. Note that the image quality is comparable when interpolation is performed using input data with a reasonably coarse trace spacing — from a 20-m to 10-m interval or from a 40-m to 20-m interval. The image quality however is not faithfully preserved when interpolation is performed using input data with very coarse trace spacing — from an 80-m to 40-m interval.

Figure 7.2-27  Difference sections: (a) difference between the stacked sections in Figures 7.2-21b and 7.2-26a; (b) difference between the stacked sections in Figures 7.2-21a and 7.2-26b.

Cascaded interpolation to increasingly finer trace spacing causes accumulation of interpolation errors. Figure 7.2-26a shows the stacked section with a 20-m trace spacing after trace interpolation applied to the already interpolated stacked section in Figure 7.2-23c with a 40-m trace spacing. Compare this section with that in Figure 7.2-21b and subtract the two sections to obtain the difference section shown in Figure 7.2-27a. Note the significant energy content of the difference section indicating that cascaded interpolation produces unfavorable results. A third cascade in the interpolation yields the stacked section in Figure 7.2-26b with a 10-m trace spacing from the already interpolated section in Figure 7.2-26a with a 20-m trace spacing. Compare this section with that in Figure 7.2-21a and subtract the two sections to obtain the difference section shown in Figure 7.2-27b. The energy content of the difference section indicates further deterioration of the signal quality with an increased number of cascades in interpolation.

## 7.3 3-D poststack migration

To understand 3-D migration, consider a point scatterer that is buried in a constant-velocity medium. The zero-offset traveltime curve in two dimensions is a hyperbola given by equation (4-4). The curvature of the hyperbolic trajectory for amplitude summation is governed by the velocity function. The equation for this trajectory derived from the geometry of Figure 4.1-15 is rewritten here as

 ${\displaystyle t^{2}=\tau ^{2}+{\frac {4x^{2}}{v_{rms}^{2}}}.}$ (12)

Refer to Figure 4.1-15 to review the 2-D migration based on diffraction summation. Given the output time τ, which coincides with the the apex of the hyperbola, compute the input time t at sample location B, and deposit the amplitude at that input location on the output section at location A, corresponding to the output time τ at the apex of the hyperbola. Assuming a horizontally layered earth model, the velocity function used to compute the traveltime trajectory described by equation (12) is the rms velocity at the apex of the hyperbola at time τ (normal moveout).

We can imagine that the zero-offset response of a point scatterer in three dimensions is a hyperboloid of revolution described by the equation

 ${\displaystyle t^{2}=\tau ^{2}+{\frac {4(x^{2}+y^{2})}{v_{rms}^{2}}},}$ (13)

where x and y are the inline and crossline coordinates of the input sample at two-way zero-offset time t. Migration in three dimensions amounts to summing amplitudes at times t over the surface of the hyperboloid and placing the result at time τ that coincides with the apex of the hyperboloid.

In migration principles, we learned that the simple diffraction summation technique for migration can be improved by making the appropriate amplitude and phase corrections based on the far-field term of the Kirchhoff integral solution to the scalar wave equation before summation. As for the 2-D case (equation 4-5), the output 3-D image Pout(x0, y0, z = /2, t = 0) from Kirchhoff summation at a subsurface location (x0, y0, z) is computed from the 3-D zero-offset wavefield Pin(x, y, z = 0, t), which is measured at the surface (z = 0), by a summation over an areal aperture in the inline and crossline directions given by

 ${\displaystyle P_{out}={\frac {\Delta x\Delta y}{4\pi }}\sum _{y}\sum _{x}\left[{\frac {\cos \theta }{v_{rms}r}}\rho (t)\ast P_{in}\right],}$ (14)

where vrms is the rms velocity at the output point (x0, y0, z), and ${\displaystyle r={\sqrt {(x-x_{0})^{2}+(y-y_{0})^{2}+z^{2}}},}$ which is the distance between the input (x, y, z = 0) and the output (x0, y0, z) points. The output image Pout is computed at (x0, y0, z = /2, t = 0) using the input wavefield Pin at (x, y, z = 0, t − r/v). For a formal mathematical treatise of the Kirchhoff integral solution to the 3-D scalar wave equation, see Section H.1.

Equation (14) can be used to compute the wavefield at any depth z. To obtain the migrated section at an output time τ, equation (14) must be evaluated at z = /2 and the imaging principle must be invoked by mapping amplitudes of the resulting wavefield at t = 0 onto the migrated section at output time τ. The complete migrated section is obtained by performing the summation in equation (14) over an areal aperture and setting t = 0 for each output location. For 2-D migration, as many as 300 traces may be included in the summation. In three dimensions, this implies that as many as 70,000 traces may need to be included in the summation. Among the early works in 3-D migration are French [25] who provides an excellent and physically intuitive description of the process based on wavefield modeling, and Schneider [26] who provides a detailed theoretical discussion on 3-D Kirchhoff migration with some field data examples.

The rho filter ρ(t) in equation (14) corresponds to the time derivative of the measured wavefield, which yields the 90-degree phase shift and adjustment of the amplitude spectrum by the ramp function ω of frequency (Table A-1, Section A.1). Since the rho filter is independent of the spatial variables, it actually can be applied to the output of the summation in equation (14). Finally, the far-field term in equation (14) is proportional to the cosine of the angle of propagation (the directivity term or the obliquity factor) and is inversely proportional to vr (the spherical spreading term) in three dimensions.

### Separation versus splitting

We now discuss some practical alternatives to summing over the surface of the hyperboloid. The diffraction hyperboloid, which is strictly associated with a constant-velocity medium, has a hyperbolic cross-section in any azimuthal direction. A two-stage summation over the surface of the hyperboloid can be performed as follows:

1. Sum along the hyperbolic cross-sections in the inline direction and place the summed amplitudes at the local apexes of these hyperbolas. The hyperboloid now is collapsed to a hyperbola that is in the plane perpendicular to the direction of the first summation. This hyperbola comprises the first summed amplitudes at the local apexes and is contained in the plane of the crossline direction.
2. Next, sum the energy along this hyperbola and place it to its apex, which is also the apex of the original hyperboloid and is where the image should be placed.

Gardner et al. (1978) proposed this two-pass approach to do 3-D migration that involves successive 2-D migrations of the inlines followed by 2-D migrations of the crosslines. A rigorous theoretical treatise of the two-stage Kirchhoff summation for 3-D zero-offset migration was later given by Jacubowicz and Levin [27].

The underlying mathematical notion for the two-pass approach is full separation of extrapolation operators in the inline and crossline directions [28]. Another method of 3-D migration is based on the numerical procedure called splitting [28] [29]. Here, extrapolation operators are independently applied to the data in both the inline and crossline directions at each downward-continuation step (Section G.1). In contrast, full separation calls for complete migration in one direction — for instance, the inline direction, before any operation is applied in the orthogonal direction.

Figure 7.3-1  Algorithms for 3-D migration based on the splitting and separation techniques of implementing migration operators.

The sequence of numerical operations in the splitting and separation methods is shown in Figure 7.3-1. Note that the number of operations in both schemes is the same, only that the order of the operations differ. It turns out that if you have strong dependency of velocities on spatial variables in the inline and crossline directions, then splitting is more accurate than separation. So splitting would be the required scheme for depth migration, whereas separation may be acceptable for time migration. As noted earlier, 3-D migration based on full separation is known as the two-pass approach. Three-dimensional migration based on splitting is considered a one-pass approach.

The two-pass or one-pass strategy for doing 3-D poststack migration should not be confused with the type of algorithm — finite-difference, Kirchhoff summation, or frequency-wavenumber used in migration. The two-pass approach based on separation is now outmoded by the one-pass approach based on the splitting of the 3-D extrapolation operator with various improvements in accuracy and efficiency. Li [30] introduced a correction term to the one-pass implicit scheme to compensate for the accumulated errors resulting from splitting. Li’s correction was adapted later to explicit schemes by Etgen and Nichols [31] (Section G.2). One-pass algorithms based on stable explicit extrapolation operators designed in frequency-space [32] [33],[34] [35] [36] [37] [38] and frequency-wavenumber domains [39] [40] have now become more widely used than those based on implicit schemes. Before we discuss the more modern one-pass explicit schemes, however, for their historical significance, we shall first discuss the two-pass finite-difference implicit scheme based on separation and one-pass implicit finite-difference scheme based on splitting of the 3-D extrapolation operator without any additional correction term.

As for 2-D migration (migration principles), a convenient domain to do 3-D time or depth migration is the frequency-space domain. Since the need for 3-D imaging becomes especially significant in areas with steep dips and lateral velocity variations, the 2-D frequency-space algorithm (Section D.4) adapted for one-pass 3-D migration is used in the examples shown in this section. Extrapolation operators based on rational approximations to the exact dispersion relation are applied to the 3-D stacked volume of data in the inline and crossline directions. A concise theoretical description of 3-D migration based on the one-pass approach is provided in Section G.1.

### Impulse response of the one-pass implicit finite-difference 3-D migration

The ideal impulse response of a 3-D migration operator is a hollow hemisphere, with semicircular cross-sections in the vertical planes and circular symmetry on time slices. Figures 7.3-2 and 7.3-3 show the 3-D 15-degree and 45-degree operator impulse responses for the one-pass approach. The corresponding 2-D operators are shown in Figure 4.2-1. Selected lines and time slices are shown to better illustrate the shape of the responses. There are 101 inlines and 101 crosslines.

Figure 7.3-2  A 15-degree split 3-D migration operator impulse response. See Figure 4.2-1 for the 2-D equivalent response. Except for the amplitudes, the vertical cross-section at center line (51) is the same as the 2-D impulse response.

The center line 51 exhibits the elliptical shape of the 15-degree 2-D operator. As we go to lines away from the center line, the ellipse naturally gets smaller in size. The time slices of the 3-D 15-degree operator retain the circular symmetry. But the traveltimes are not quite correct — it causes undermigration in all directions.

Examine the time slices from the 15-degree 3-D impulse response (Figure 7.3-2) and note that the response is circularly symmetric. In contrast, the time slices of the 3-D 45-degree operator (Figure 7.3-3) do not have the circular symmetry; instead, they exhibit a diamond shape. The 15-degree operator causes the same amount of undermigration in all azimuthal directions. The 3-D 45-degree operator causes less undermigration, but the amount of migration error depends on the azimuth — the largest error being in the diagonal directions.

Figure 7.3-3  A 45-degree split 3-D migration operator impulse response. See Figure 4.2-1 for the 2-D equivalent response. Except for the amplitudes, the vertical cross-section at center line (51) is the same as the 2-D impulse response.

Although the one-pass algorithm is appropriate to handle lateral velocity variations — hence, suitable for depth migration as well as for time migration — it has problems of its own. The 15-degree 3-D operator is circularly symmetric (Figure 7.3-2), but when you ask for more steep-dip accuracy, as with the 45-degree 3-D operator, you end up with azimuthal variations in migration error that is manifested in the form of diamond-shaped time slices in the impulse response (Figure 7.3-3).

What is the remedy for the diamond-shaped behavior on the time slices of the one-pass 3-D 45-degree operator (Figure 7.3-3)? Ristow [41] discusses further splitting of the 3-D migration operator in the two diagonal directions, in addition to the inline and crossline directions, as a way to achieve a more circular time slice. Ristow [41] also gives a least-squares approximation (for a specified dip angle) to the 3-D migration operator by two 2-D operators. The problem with this technique is that data need to be sorted and operated upon in the two diagonal directions. This not only increases the computational cost, but also causes the extrapolation of smaller and smaller spatial data arrays in the proximity of the corners of the survey area.

Li [30] examined the difference between the exact 3-D phase-shift extrapolation operator and the one-pass extrapolation operator, and derived a residual extrapolation term to shape the diamond response of the 45-degree operator (Figure 7.3-3) to a circular response of a desired 3-D operator.

Certainly, the best way to achieve circular symmetry is by way of the 3-D phase-shift extrapolation operator (Section G.3). To accommodate lateral velocity variations, a residual correction can be applied at each depth level before extrapolating down to the next depth level. This algorithm, known as phase-shift-plus-correction scheme [39] [40], will be discussed later in the section.

Another way to circumvent the diamond shape is to Fourier transform the wavefield in the direction with the milder velocity variations, say in the crossline direction, and perform 2-D migration for each of the crossline wavenumbers [42]. Again, this would require the constant-velocity assumption in the direction of Fourier transformation, something contrary to the reason for doing 3-D surveys and doing 3-D migration.

One other suggested solution [43] [44] is to compute the exact 3-D phase shift operator in the Fourier transform domain, inverse Fourier transform in the inline and crossline directions, and then truncate to get a finite-size complex operator in the frequency-space domain. This would require the constant-velocity assumption over the spatial extent of the operator in the inline and crossline directions. Nevertheless, you would not have the errors in difference approximations to differential operators as with the finite-difference algorithms, since you are computing the derivatives of the wavefield in the inline and crossline directions in the Fourier transform domain. For the 3-D wave extrapolation from one depth step to another, the 2-D operator in the inline-crossline domain is then convolved with the complex wavefield for each frequency component.

While this scheme based on direct convolution of the extrapolation operator with the wavefield at one depth level is more accurate than the one-pass scheme based on operator splitting, it is computationally expensive. It can only be made practical if truncation errors are minimized, the 3-D extrapolation operator can be made spatially compact, and is applied, efficiently. Indeed, 2-D stable explicit extrapolation operators can be designed with minimal truncation errors within a specified dip range [32] [33] (Section D.5). Blacquiere [45] extended the explicit method by Holberg [32] to 3-D poststack depth migration. This extension includes computing and tabulating the coefficients of 2-D explicit extrapolation operators in the frequency-space domain. Again, the operator is applied to the wavefield at each depth level by direct convolution. Finally, Hale [34] devised an algorithm for 3-D poststack migration that circumvents the cost of computing and tabulating 2-D operators, and applying them by direct convolution. As will be shown later in the section, a 2-D explicit operator can be used to create a steep-dip 3-D extrapolation operator that preserves near-circular symmetry by way of the McClellan transform [34]. This technique only requires computing 1-D explicit extrapolation filter coefficients and uses an efficient recursive filtering scheme based on Chebychev polynomials to perform the convolution of the extrapolation operator with the input wavefield.

### Two-pass versus one-pass implicit finite-difference 3-D migration in practice

The two-pass approach for 3-D migration is valid strictly for a constant-velocity medium [41] [27]. Since we never have a constant-velocity medium, it is reasonable to question the practical value of the two-pass 3-D migration. Now consider the case of the vertically varying velocity model shown in Figure 5.1-17.

Compute the 3-D zero-offset response of the three point scatterers buried in the subsurface below the center midpoint (32) of the center inline (32) at the layer boundaries (denoted by the asterisks in the velocity model in Figure 5.1-17). Figure 7.3-4a shows four selected inlines from this 3-D zero-offset data that consist of 63 inlines and 63 crosslines. Obviously the traveltime responses are circularly symmetric, the shallow one being an exact hyperboloid, while the two deeper ones are approximately hyperboloid.

One-pass 3-D migration collapses the energy to the apexes of the three diffraction traveltime surfaces at the center midpoint (32) of the center line (32) (Figure 7.3-4b). Now consider two-pass 3-D migration. First migrate in the inline direction (Figure 7.3-4c). The center inline (32) is migrated properly, but the lines farther from the center line are increasingly overmigrated. Moreover, within one line, say, line 2, the shallower the diffraction, the more the overmigration. This is easy to understand once we refer to the velocity model (Figure 5.1-17).

On lines farther from the center line, the 3-D diffraction events arrive at increasingly late times; for example, on line 2 in Figure 7.3-4a. Because the velocity varies with depth, diffractions on line 2 get migrated with velocities higher than the true velocities associated with these diffractions. The diffractions on line 2 are sideswipes from the point scatterers situated beneath the center line (32). Because velocities corresponding to late times are used to migrate energy generated by shallow diffractors, these sideswipes are overmigrated. Now sort the data into crosslines (Figure 7.3-4d). It is necessary to display only those lines closer to the center line. Note that most of the energy has collapsed to the hyperbolas on the center line (32); however, because of the error induced by the first pass caused by the effect of v(z), some energy still remains on the neighboring lines. Also note that the error is more significant at shallow diffraction apexes. Finally, migration in the crossline direction collapses most of the energy to the center midpoint (32) (Figure 7.3-4e). Sort the data back to inlines (Figure 7.3-4f) and compare them with the one-pass result (Figure 7.3-4b).

Energy spreading to the neighboring midpoints results from the error of the two-pass approach in the presence of vertical velocity variations. The test results shown in Figure 7.3-4 also suggest that the larger the vertical velocity gradient, the more error induced by the two-pass approach. Another important point is that with the two-pass approach and for a v(z) function, it makes a difference in which direction migration is done first. Smearing is greatest in the direction in which the first pass of 3-D migration is performed.

Figure 7.3-5 shows selected time slices from the unmigrated volume of data (Figure 7.3-4a), and from the output of one-pass (Figure 7.3-4b) and two-pass (Figure 7.3-4e) migrations. The time slices at 400, 600, and 800 ms coincide with the apex times of the zero-offset hyperboloidal traveltime surfaces at the center line 32 (Figure 7.3-4b). Note that the time slices also show that much of the energy has collapsed to the apexes of the hyperboloids. The residual energy left behind however appears to be more on the two-pass result than on the one-pass result.

Theoretically, it is true that the two-pass method only is appropriate for velocities that vary slowly in the vertical direction and do not vary laterally. There is also the question of migration algorithm accuracy versus accuracy in the velocity field used for migration. In many practical situations, error induced by the two-pass approach probably is less than the possible errors associated with the uncertainty in the velocity field. Hence, it may be difficult on real data to tell whether imperfect imaging results from the velocity errors or the two-pass approach. In practice, we usually find that the two-pass approach yields acceptable results in areas in which dips are small, vertical velocity variations are moderate, and lateral velocity variations are within the limits of time migration.

Now consider the field data example in Figure 7.3-6. The top left is an inline and the top right is a crossline stacked section from a land 3-D survey. Migration in the inline direction yields the sections shown in the center row. Note the partially collapsed diffraction energy on the inline section (center left) and a pronounced enhancement of the diffraction energy on the crossline section (center right). Sorting the data and migrating in the crossline direction yields the results shown in the bottom row. These are the two-pass 3-D migrated sections.

A note of theoretical significance is that the crossline section after 2-D migration in Figure 7.3-6 can be considered as a true 2-D wavefield within the bounds of the zero-offset assumption. This is because it does not have sideswipe energy contained in the original unmigrated stacked section, which is a 2-D cross-section of a true 3-D wavefield, again within the bounds of the zero-offset assumption.

Figure 7.3-7 shows the comparison between the two-pass (center row) and one-pass (bottom row) 3-D migrations of the same inline and crossline sections (top row). The apparent overmigration (just below location A) seen on the two-pass result may be attributed to the error of the two-pass approach resulting from large vertical velocity gradients associated with these data.

Note that the above analysis is entirely within the framework of time migration. Two-pass 3-D migration algorithms cannot even be considered approximately accurate in situations where lateral velocity variations are large enough to warrant depth migration.

We now summarize the significant differences between one-pass and two-pass 3-D migrations based on these experiments with synthetic and field data. Also, we shall make reference to the results of the quantitative work by Dickinson [46] with regard to the two-pass approach. We remind ourselves that, here, we refer to the two-pass approach based on separation and the one-pass approach based on splitting of the 3-D extrapolation operator without any correction term to compensate for errors in extrapolation such as that by Li [30].

1. Two-pass 3-D migration causes overmigration, and one-pass 3-D migration causes undermigration of steeply dipping events in a medium with significant vertical velocity gradients. This is illustrated using a circularly symmetric salt dome model in Figure 7.3-8.
2. Error in two-pass 3-D migration is negligible for small dips.
3. For steep dips, the error is the same order of magnitude as the error caused by a lack of precise knowledge of migration velocities.
4. Largest overmigration with the two-pass approach usually occurs at around the 45-degree azimuth. No error occurs when the dip is entirely in the inline or crossline directions — the case of 2-D rather than 3-D (Figure 7.3-8).
5. Largest undermigration with the one-pass approach also occurs at around the 45-degree azimuth (Figure 7.3-8).
6. The two-pass approach translates the dipping event out of its true position. After the two-pass 3-D migration, the dip is correct in contrast to overmigration caused by erroneously too high velocities which not only mispositions the event but also changes its dip.
7. Finally, the two-pass 3-D migration results depend on which direction the first-pass migration has been performed — if it is the inline direction, and if the velocity gradient is large enough for the two-pass approach to cause overmigration, then you get more overmigration in that direction.

The undermigration effect of one-pass 3-D migration illustrated in Figure 7.3-8 can be observed in the field data example shown in Figures 7.3-9 and 7.3-10. The time slices from the stacked volume of data infer a nearly circular salt dome (left column in Figure 7.3-9). The time slices from the desired image volume obtained from 3-D poststack phase-shift migration also exhibit the near-circular shape of the salt dome (right column in Figure 7.3-9). Note how the size of the circular events gets smaller after migration on the same time slice as was illustrated in Figure 7.3-8. Figure 7.3-10 shows a comparison of the time slices from the image volumes obtained from one-pass 3-D migration using the frequency-space 45-degree implicit scheme in split mode (left column) and 3-D phase-shift migration (right column). It may be claimed that the near-circular shape especially evident in the 1600-ms and 1700-ms time slices of the image volume from 3-D phase-shift migration is not quite preserved by the one-pass 3-D migration. Instead, we see a distorted circle much the same as sketched in Figure 7.3-8b for the one-pass 45-degree scheme. This distortion suggests undermigration mostly in the two diagonal directions caused by the one-pass schemes. At deeper time slices (Figure 7.3-10), the differences between the one-pass and phase-shift schemes are increasingly more subtle because the dip of the salt flank decreases at late times.

Note the practical significance of the last aspect on the list. For instance, if you decide to do two-pass 3-D migration in an overthrust area, then you may want to migrate first in the strike direction where the velocity variations may generally be changing less rapidly than the dip direction to minimize the overmigration effect of the two-pass approach. Also, once you complete the migration in the direction of mild velocity variations, you can afford to revise the second pass in the direction of strong velocity variations without going back to the original stacked data.

### Explicit schemes combined with the mcclellan transform

One-pass and two-pass 3-D poststack migration algorithms based on implicit finite-difference schemes are now scarcely used. The primary reason is the azimuth-dependent positioning errors they incur (Figure 7.3-8). Improvements to the one-pass scheme to achieve circular symmetry in the impulse response [41] [30] are computationally involved. Instead, explicit schemes have gained a wide acceptance since the mid-90s. They provide circular symmetry and efficiency in the design and application of the extrapolation operators [33],[34] [35] [36] [37] [38]. As with the implicit methods, explicit methods can be formulated to perform both time and depth migration of 3-D CMP stacked data.

There are three aspects to the design and application of 3-D explicit schemes.

1. Design of a 2-D explicit operator in the frequency-space domain,
2. Transformation of this 2-D operator to a 3-D wave extrapolator, and
3. Application of the 3-D wave extrpolator to migrate 3-D CMP-stacked data.

Here, we shall review an explicit method based on the McClellan transform [34]. An alternative design of 2-D explicit operators based on the Remez exchange algorithm (Section D.5) and the Laplace transform to convert the 2-D operator into a 3-D wave extrapolator is given by Soubaras [47]. Table 7-2 provides a list of combinations of schemes used in the design of 2-D explicit operators, transforming them into 3-D wave extrapolators and the application of the 3-D extrapolators to migrate 3-D CMP-stacked data. Mathematical details of the design and application of 3-D explicit extrapolation operators are left to Section G.2.

As for the one-pass implicit schemes, explicit schemes that use the McClellan transform for 3-D poststack migration are implemented in the frequency-space (ω, x, y) domain. This enables application of the extrapolation filter at each depth level to each frequency component of the wavefield, independently.

Wave extrapolation in 3-D is described in the frequency-wavenumber domain by

 ${\displaystyle P(k_{x},k_{y},z+\Delta z,\omega )=P(k_{x},k_{y},z,\omega )\exp(-ik_{z}\Delta z),}$ (15)

where, kx and ky are inline and crossline wavenumbers, respectively, kz is the vertical wavenumber, and ω is frequency in units of radians per unit time, and P is the wavefield that is being extrapolated in depth z using a depth step size Δz. The desired extrapolation operator D(k) = exp(−ikzΔz) for a specific ratio of 2ω/v is given by

 ${\displaystyle D(k)=\exp \left[-i{\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk}{2\omega }}\right)^{2}}}\Delta z\right],}$ (16)

where

 ${\displaystyle k={\sqrt {k_{x}^{2}+k_{y}^{2}}}.}$ (17)
 2-D design Least-squares Remez exchange 3-D design McClellan transform Laplace transform 3-D application Chebychev recursion Chebychev recursion

We want to compute an extrapolation filter hn in the frequency-space (ω, x, y) domain with a Fourier transform H(k, ω, v)

 ${\displaystyle H(k)=h_{0}+2\sum _{n=1}^{N}h_{n}\cos(nk),}$ (18)

that best approximates the desired transform given by equation (16) (Section D.5). Since D(k) is symmetric with respect to k = 0, the complex filter coefficients hn are even — h−n = hn, and the number of filter coefficients 2N + 1 is odd. As a result, we only need to compute N coefficients.

Methods to obtain hn such that the actual Fourier transform H(k) given by equation (18) closely approximates the actual Fourier transform D(k) given by equation (16) is described by Holberg [32] and Hale [33]. The former is based on least-squares and the latter is based on a modified Taylor series. The Hale method is described in Section D.5. Since the desired amplitude spectrum is |D(k)| = 1 for all k, then we require the actual amplitude spectrum to be |H(k)| ≤ 1 for all k. This will ensure that the operator is stable — that is, amplitudes of the extrapolated wavefield will not grow from one depth level to another.

Following the determination of the 1-D extrapolation filter for specific ω and v values, they are convolved with the McClellan transform filter template (Tables G-1 and G-2) to translate it to a 2-D filter in the (x, y) domain for each frequency ω. The filter is then applied to the appropriate frequency component of the wavefield P(x, y, ω, z) in the (ω, x, y) domain at depth z to extrapolate it to depth z + Δz by way of an efficient recursive scheme that is based on Chebychev polynomials. At each depth, the image P(x, y, z + Δz, t = 0) is obtained by invoking the imaging principle that is equivalent to summing the extrapolated wave components over frequency.

In practice, the usual implementation of the McClellan transform method requires equal spatial sampling intervals, Δx and Δy, in inline and crossline directions, respectively. When the two sampling intervals are not equal, as in most marine 3-D surveys, then trace interpolation needs to be done prior to migration to obtain data with equal sampling intervals in inline and crossline directions.

Figure 7.3-11 shows the desired impulse response of a 3-D poststack migration operator. There are 101 inlines and 101 crosslines, with line 51 being the center line. This impulse response was generated using the 3-D Stolt migration with constant velocity (Section G.4). Hence, it is accurate for all dips and frequencies. Note from the time slices the perfect circular symmetry. We also achieved the circular symmetry with the implicit 15-degree split operator (Figure 7.3-2), but failed with the implicit 45-degree split operator (Figure 7.3-3). In fact, using a steeper-dip implicit split operator also fails to respond with circular symmetry (Figure 7.3-12). In the case of Figure 7.3-12, the steep-dip implicit scheme was designed by the cascaded application (twice) of the 65-degree operator (Section D.5). Note from the vertical cross-sections the dispersive noise inherent to implicit schemes. Specifically, implicit schemes cause evanescent energy that behaves like propagating energy, rather than letting it vanish with depth. The time slices exhibit a distinctive diamond shape meaning that positioning errors in the form of undermigration with the splitting methods vary with azimuth. The largest error occurs along the two diagonals, and the least error occurs along the inline and crossline directions. Both dispersive noise and azimuthal asymmetry make the impulse response in Figure 7.3-12 undesirable.

Now, examine the impulse response of the 3-D migration operator based on the explicit scheme discussed in Section D.5 and the McClellan transform discussed in Section G.2 (Figure 7.3-13). The 1-D explicit operator hn has 2N + 1 = 39 filter coefficients, and the McClellan filter template is 3 × 3 in size (Table G-1). The stable explicit operator hn provides cleaner vertical cross-sections, free of dispersive noise compared to the implicit scheme (Figure 7.3-12). The explicit scheme does not allow the evanescent energy to turn into propagating energy. The dip accuracy of the explicit scheme is governed by the length of the filter hn — the longer the filter the steeper the dip the algorithm can handle. The McClellan transform provides the more circular time slice compared to the splitting method (Figure 7.3-12).

Figure 7.3-14 shows the vertical cross-sections and time slices of the impulse response of another 3-D migration operator based on the McClellan transform. The 1-D explicit operator hn has 2N + 1 = 39 filter coefficients and the McClellan filter template is 5 × 5 in size (Table G-2). Because we have the same number of coefficients (39) in the explicit 1-D filter hn, the steep-dip accuracy is the same in both Figures 7.3-13 and 7.3-14. However, the McClellan filter template used in the impulse response of Figure 7.3-14 provides a more circular behavior at higher cost than that used in the impulse response of Figure 7.3-13.

We now perform 3-D migration of the 3-D zero-offset synthetic data set in Figure 7.3-4a. The velocity field is based on a horizontally layered earth model shown in Figure 5.1-17. Figure 7.3-15 shows selected vertical sections of the migrated output from the explicit scheme combined with the McClellan transform. The wave extrapolation in this case was performed in true depth z. The time slices from the input 3-D zero-offset wavefield and the depth slices from the 3-D image output of 3-D migration are shown in Figure 7.3-16. Note that energy associated with the three point scatterers has collapsed onto the apexes of the traveltime surfaces, which coincide with the depths where the scatterers are located.

### The phase-shift-plus-correction method

The phase-shift method (Section G.3) is designed to accommodate vertically varying velocity field for migration, only. Nevertheless, it can be extended to accommodate lateral velocity variations [39] [48]. Consider a spatially varying 3-D migration velocity field v(x, y, z) and a laterally averaged, but vertically varying velocity function ${\displaystyle {\bar {v}}(z)}$. The basic idea is to first extrapolate in depth with the phase-shift extrapolator (equation G-22) using the vertically varying velocity function ${\displaystyle {\bar {v}}(z)}$. This is followed by the application of a convolutional operator as in the explicit schemes (Section G.2) to the extrapolated wavefield at the same depth. This second operation is fundamentally a residual extrapolation to account for the difference between the laterally varying velocity field v(x, y, z) and the vertically varying velocity function ${\displaystyle {\bar {v}}(z)}$. The migration algorithm that incorporates such a correction term into the phase-shift algorithm has been known as phase-shift-plus-correction (PSPC) method. Such a splitting of wave extrapolation within one depth step may appear to be computationally awkward. Nevertheless, numerically efficient schemes have been devised to apply the convolutional operator in the residual extrapolation step [39]. Depending on the degree of lateral velocity variations, the PSPC method can be used either as a time or depth migration algorithm.

Figure 7.3-17 shows selected crosslines from the 3-D migrated volumes of data associated with the 3-D survey of Figure 7.2-1 using three different migration algorithms — the implicit one-pass on top, the explicit McClellan transform in the middle, the PSPC at the bottom. Note that both of the explicit schemes based on the McClellan transform and the PSPC methods have produced comparable and better images of the subsurface in the vicinity of the overthrust between inline locations 200 and 300. The same region of the subsurface appears to be undermigrated by the implicit one-pass scheme. Incorrect positioning by the implicit one-pass scheme also is pronounced in the inline sections shown in Figure 7.3-18. Note, for instance, poor reflector continuity between 1.5-2 s and crossline locations 200 and 400 — events have not been moved completely into this section from the neighboring sections. Moreover, note the wobbly character of these events — this is caused by the circularly asymmetric impulse response of the implicit one-pass operator (Figure 7.3-12). Notice also the generally noisy character of the results from the implicit one-pass scheme — this may be attributed to the dispersive noise produced by the algorithm. Specifically, it treats evanescent energy as if it is propagating energy and generates a noisy impulse response at steep dips (Figure 7.3-12).

Figure 7.3-19 shows selected time slices from the results of migration using the three different migration algorithms — the implicit one-pass on the left, the explicit McClellan transform in the middle, the explicit phase-shift-plus-correction on the right. Note the wobbly behavior of the contours on the time slices from the implicit one-pass scheme — this is related to the azimuthal asymmetry of its impulse response (Figure 7.3-12). The differences are exhibited more clearly on the enlarged view in Figure 7.3-20. The two explicit schemes — the McClellan transform and the phase-shift-plus-correction, produce comparable images.

Study the crosslines (Figure 7.3-17), the inlines (Figure 7.3-18), and the time slices in Figure 7.3-19, and note the differences in imaging by the implicit and explicit schemes. It is almost certain that in the future, explicit schemes, because of their ease of design and implementation, will become standard. Increased computer power will further encourage use of the explicit schemes.

## 7.4 3-D prestack time migration

In prestack time migration, we discussed 2-D prestack time migration as the rigorous solution to the problem of conflicting dips with different stacking velocities. When the reflector geometries that give rise to this problem have a 3-D behavior, it is necessary to image the subsurface using 3-D prestack time migration. Fault-plane reflections associated with rotated fault blocks are one such case that requires imaging in 3-D and before stack. Aside from solving the problem of conflicting dips with different stacking velocities and thus producing an improved migrated stack volume, 3-D prestack time migration produces common-reflection-point (CRP) gathers which can be used for amplitude variation with offset analysis.

As for the 2-D case, the robust alternative to 3-D prestack time migration in practice is to apply NMO and 3-D DMO corrections followed by 3-D poststack time migration. We shall begin this section by presenting a widely accepted procedure for 3-D prestack time migration that is based on a modification to this robust alternative. Specifically, it is assumed that following NMO and 3-D DMO corrections, each of the common-offset volumes of data is equivalent to a 3-D zero-offset wavefield, and thus can be migrated using a 3-D poststack time migration algorithm, individually. Also in this section, we shall review crossline migration as a process that reduces the imaging problem from three dimensions to two dimensions [49] [50].

Figure 7.4-1 shows three inline sections from a 3-D CMP-stacked data volume associated with a marine 3-D survey with no 3-D DMO correction. Four time slices associated with the 3-D CMP-stacked data volume as in Figure 7.4-1 are shown in Figure 7.4-2. From the inline sections and the time slices, note that the subsurface structural setting involves highly complicated fault patterns. The cross-sections of the image volume from 3-D poststack time migration along the same inline traverses as in Figure 7.4-1 are shown in Figure 7.4-3. Although the fault planes themselves are not delineated, the presence of the complicated fault patterns is evident from the cross-sections (Figure 7.4-3) and the time slices (Figure 7.4-4) of the image volume.

A subset of the velocity spectra used to pick stacking velocities is shown in Figure 7.4-5. Figure 7.4-6 shows the cross-sections of the 3-D migration velocity field along three inline and three crossline traverses. This velocity field was derived by interpolating the vertical functions picked from the velocity spectra as in Figure 7.4-5 and spatially smoothing the resulting velocity volume.

Now apply 3-D DMO correction to preserve the fault-plane reflections. Figure 7.4-7 shows the three inline sections and Figure 7.4-8 shows the four time slices from the 3-D DMO-stack data volume. The cross-sections of the image volume from 3-D poststack time migration along the same inline traverses as in Figure 7.4-7 are shown in Figure 7.4-9 and the time slices from the same volume are shown in Figure 7.4-10. Note that 3-D DMO correction preserves the fault-plane reflections on stacked data and the subsequent 3-D poststack time migration delineates the fault blocks with improved imaging (compare Figure 7.4-3 with Figure 7.4-9).

A subset of the velocity spectra used to pick velocities after 3-D DMO correction is shown in Figure 7.4-11. Figure 7.4-12 shows the cross-sections of the 3-D migration velocity field along three inline and three crossline traverses. This velocity field was derived by spatially interpolating the vertical functions picked from the velocity spectra as in Figure 7.4-11.

The kinematics of 3-D prestack time migration can be formulated as an extension of the kinematics of 2-D prestack time migration (migration velocity analysis and E.5) in the same manner as for 3-D DMO correction. In fact, Figure 7.2-11, which describes the 3-D DMO process, can also be used to describe 3-D prestack time migration. Consider a trace from a common-cell gather with no NMO nor DMO correction and map the amplitude at time sample A to neighboring cells which are coincident with the source-receiver azimuthal direction associated with that input trace along the semi-elliptical trajectory that describes the kinematics of 3-D prestack time migration. The nonzero-offset migration ellipse is described by equation (5-33). Repeat the process for all the traces from the same common-cell gather and map the amplitudes in the same manner.

For the hypothetical recording geometry with a common-cell gather coincident with a common-midpoint gather and traces in the gather covering a 360-degree source-receiver azimuthal range, but having the same source-receiver separation, the elliptical trajectories associated with all the traces constitute an ellipsoid of revolution as shown in Figure 7.2-11a. This ellipsoid describes the kinematics of the impulse response of a 3-D prestack migration operator.

In migration principles, we discussed semicircle superposition and diffraction summation concepts for 2-D zero-offset migration. Applying those concepts to the case of 3-D zero-offset migration, the latter can be conceptualized as the spreading of amplitudes on each input stacked trace in the (x, y, t) volume over the surface of a hollow hemisphere (Figure 7.3-11). Superposition of the resulting hemispherical surfaces yields the (x, y, z) image volume. Alternatively, for a given output sample of a trace in the (x, y, z) image volume, amplitudes over the surface of the hyperboloid of revolution in the (x, y, t) volume of the input 3-D zero-offset wavefield can be summed and placed on that output sample location. The Kirchhoff summation technique for migration incorporates the amplitude and phase factors described in Section H.1.

Similarly, for an input trace with a specific source-receiver azimuth and offset as depicted in Figure 7.2-11a, 3-D prestack time migration can be conceptualized either by way of a semi-elliptical superposition using equation (5-33) or a diffraction summation over the traveltime trajectory described by equation (5-32). The recording geometry of Figure 7.2-11a, however, never exists in reality. Instead, 3-D recording geometries give rise to nonuniform source-receiver azimuthal and offset coverage, and midpoint scattering over the survey area. 3-D DMO correction implicitly regularizes the spatial sampling of the prestack 3-D data. As a result, the data can be decoupled and sorted into common-offset volumes each of which is considered a replica of a 3-D zero-offset wavefield. This then enables us to adopt the robust approach for 2-D prestack time migration based on DMO correction and common-offset migration (prestack time migration) to develop an efficient workflow for 3-D prestack time migration as described below.

### 3-D DMO correction combined with 3-D common-offset migration

In prestack time migration, we learned a prestack time migration strategy based on DMO correction combined with common-offset migration [51]. Aside from prestack imaging, this strategy provides the opportunity to derive a velocity field associated with events in their migrated positions. We shall now develop a workflow for 3-D prestack time migration based on 3-D DMO correction combined with 3-D common-offset migration. This workflow yields common-reflection-point (CRP) gathers which can be stacked to produce the image volume from 3-D prestack time migration. The CRP gathers can also be used to derive a 3-D rms velocity field associated with the migrated data [52].

1. Starting with input prestack data, apply NMO correction using flat-event velocities. These are picked from velocity spectra computed over a sparse grid as in Figure 7.4-5.
2. Apply 3-D DMO correction and sort the data into common-offset volumes. Figures 7.4-13, 7.4-14, and 7.4-15 show three common-offset sections with 400-m, 1000-m, and 1800-m offsets after 3-D DMO correction along the inline traverses as in Figure 7.4-7. These common-offset sections exhibit steeply dipping reflections associated with the fault planes which conflict with the gently dipping reflections associated with the sedimentary strata.
3. Following NMO and 3-D DMO correction, each common-offset volume is assumed to be a replica of a 3-D zero-offset section, and thus, can be migrated using a 3-D zero-offset migration algorithm. A spatially averaged but vertically varying velocity function can be used to perform the migrations of the common-offset volumes of data. The same common-offset sections as in Figure 7.4-13, 7.4-14, and 7.4-15 after 3-D common-offset migration are shown in Figures 7.4-16, 7.4-17, and 7.4-18.
4. Sort the migrated common-offset volumes into CRP gathers.
5. Apply inverse NMO correction using the velocity field as in step (a) and repeat the velocity analysis. A subset of the velocity spectra used to pick stacking velocities after 3-D common-offset migration is shown in Figure 7.4-19.
6. Create a 3-D velocity field by interpolating the vertical functions picked from the velocity spectra as in Figure 7.4-19 and spatially smoothing the resulting velocity volume. The cross-sections from the 3-D migration velocity field along three inline and three crossline traverses are shown in Figure 7.4-20.
7. Apply NMO correction to the common-reflection-point (CRP) gathers using the postmigration velocity field from step (e) (Figure 7.4-20). Selected CRP gathers are shown in Figure 7.4-21. The CRP stacks constitute the cross-sections of the image volume from 3-D prestack time migration. Three such sections along the same line traverses as in Figure 7.4-7 are shown in Figure 7.4-22, and the time slices from the image volume are shown in Figure 7.4-23. Compare the cross-sections with those from the two image volumes based on 3-D poststack time migration of the conventional CMP-stacked data volume (Figure 7.4-3) and that of the 3-D DMO stack volume (Figure 7.4-9), and note the significant improvement in the imaging of the complex fault blocks.
8. To obtain the migrated volume with the updated velocity field (Figure 7.4-20), first, perform 3-D zero-offset inverse migration (equivalent to 3-D zero-offset forward modeling) of the resulting image volume from step (g) using the same velocity function used to perform the 3-D common-offset migrations as in step (c). The cross-sections of the resulting 3-D zero-offset wavefield along the same inline traverses as in Figure 7.4-7 are shown in Figure 7.4-24. Selected time slices from this modeled data volume are shown in Figure 7.4-25 which can be compared with those from the CMP stack (Figure 7.4-2) and DMO stack (Figure 7.4-8) volumes. The modeled volume of data in Figure 7.4-22 may be treated as equivalent to the unmigrated stack volume.
9. The final step involves 3-D zero-offset migration of the modeled volume of data using the 3-D migration velocity field as in Figure 7.4-20. The cross-sections along the inline traverses as in Figure 7.4-22 from the 3-D prestack time migration sequence described above are shown in Figure 7.4-26, and the selected time slices from the image volume are shown in Figure 7.4-27. The same cross-sections as in Figure 7.4-26 with reverse polarity are displayed in Figure 7.4-28.

Compare the results of 3-D prestack time migration (Figure 7.4-26) with the results from 3-D poststack time migration of the data with (Figure 7.4-9) and without (Figure 7.4-3) 3-D DMO correction. Note the improvement in imaging the fault planes with 3-D prestack time migration.

### Crossline migration

In 3-D poststack migration, we saw that 3-D poststack migration can be performed by a two-stage summation over the surface of the 3-D zero-offset diffraction hyperboloid. In practice, the two-pass 3-D poststack migration involves successive 2-D poststack migrations of the inlines followed by 2-D poststack migrations of the crosslines. Similarly, 3-D prestack migration also may be performed in two stages. While two-pass 3-D migration is only an approximation to a one-pass 3-D migration (3-D poststack migration), the computational benefits of the former may encourage its use in practice under certain circumstances.

The 3-D nonzero-offset traveltime t (Section G.6) for a zero source-receiver azimuth represents an ellipsoid (Figure G-2) given by

 ${\displaystyle {\frac {x^{2}}{(vt/2)^{2}}}+{\frac {y^{2}}{(vt/2)^{2}-h^{2}}}+{\frac {z^{2}}{(vt/2)^{2}-h^{2}}}=1,}$ (19)

where x and y are inline and crossline directions, respectively, v is the medium velocity and h is the half-offset.

The horizontal cross-section of this ellipsoid at a depth z is an ellipse (Figure 7.4-29). For a 3-D recording geometry with an arbitrary source-receiver azimuth defined by angle θ, the inline-crossline Cartesian coordinates (x, y) in terms of the local Cartesian coordinates (x′, y′) are given by

 ${\displaystyle x=x^{\prime }\cos \theta -y^{\prime }\sin \theta }$ (20a)

and

 ${\displaystyle y=x^{\prime }\sin \theta +y^{\prime }\cos \theta .}$ (20b)

By substituting equations (20a,20b) into equation (19), it can be shown that the vertical cross-section of the ellipsoid along an arbitrary source-receiver azimuth θ is also an ellipse [49]. This characteristic of the traveltime ellipsoid makes it possible to perform 2-D migration of 3-D prestack data in one direction only, so as to create 2-D prestack data regularized in offset and azimuth [49] [50] [53]. This one direction is most likely to be the crossline direction for most cases that warrant the two-pass migration strategy in practice — thus the term ''crossline migration''. Because of the irregular offsets and source-receiver azimuths of 3-D prestack data, an appropriate choice for crossline migration algorithm is based on Kirchhoff summation.

The output gathers from crossline migration along one inline can be associated with a hypothetical 2-D line in that inline direction with no sideswipe energy. These gathers can be processed subsequently, and a 2-D migration in the inline direction can be performed by one of the following strategies:

1. 2-D prestack time migration based on a workflow that includes 2-D DMO correction and 2-D common-offset migration (prestack time migration) to generate common-reflection-point (CRP) gathers for amplitude variation with offset (AVO) analysis (analysis of amplitude variation with offset), or
2. 2-D earth modeling in depth (earth modeling in depth) to derive an accurate image in depth (earth imaging in depth) by 2-D prestack depth migration in the inline direction that is assumed to coincide with the dominant structural dip direction.

Shown in Figure 7.4-30a is an inline section from a 3-D unmigrated zero-offset volume of data. The latter was created as part of the 3-D common-offset migration strategy for 3-D prestack time migration discussed earlier in this section. Specifically, following 3-D DMO correction and 3-D common-offset migration, data were stacked and unmigrated using the velocity field that was used in 3-D common-offset migration and an algorithm applicable to 3-D zero-offset wavefields to obtain a 3-D zero-offset wavefield. The section in Figure 7.4-30a is, therefore, a 2-D inline cross-section from the 3-D zero-offset wavefield. The 3-D zero-offset wavefield is subsequently migrated using the velocity field derived from a velocity analysis of the data after 3-D common-offset migration, and again, an algorithm applicable to 3-D zero-offset wavefields. The 2-D cross-section from the remigrated data volume along the same inline traverse is shown in Figure 7.4-30b.

Figure 7.4-29  (a) The 3-D nonzero-offset traveltime surface associated with 3-D prestack migration is an ellipsoid given by equation (19), and the vertical cross-section of the ellipsoid along an arbitrary source-receiver azimuth AA′ as denoted in (b) is an ellipse. Adapted from [54].

The 2-D zero-offset section derived from crossline migration along the same inline traverse is shown in Figure 7.4-30c, and the image from 2-D prestack time migration of the gathers from the crossline migration is shown in Figure 7.4-30d. The unmigrated section in Figure 7.4-30c represents a 2-D zero-offset wavefield and, as such, it does not contain any sideswipe energy. Whereas the unmigrated section in Figure 7.4-30a represents a cross-section from a 3-D zero-offset wavefield, and therefore, it contains sideswipe energy. The image quality of the migrated sections in Figures 7.4-30b and 7.4-30d is comparable.

A crossline migration workflow includes the following steps:

1. Perform velocity analysis of unmigrated 3-D prestack data to derive a 3-D velocity field that is appropriate for migration in the crossline direction. Often, a single, vertically varying velocity function is adequate to use in crossline migration.
2. Perform crossline migration and create common-cell gathers along selected inline traverses. The resulting gathers will have traces with regular offset distribution and uniform source-receiver azimuth that is coincident with the inline direction. From here onward, apply the 2-D prestack time migration sequence described in prestack time migration.
3. First, perform velocity analysis at selected locations along the inline traverses for which data have been migrated in the crossline direction. Figure 7.4-31 shows the common-cell gathers and the computed velocity spectra at four analysis locations along an inline traverse. The gathers are from crossline-migrated data, and as such, have traces with regular offset distribution and uniform azimuth.
4. Apply NMO and 2-D DMO corrections to the gathers from step (b) using the velocity field derived from the analysis in step (c).
5. Perform 2-D common-offset migration using an appropriately smoothed form of the velocity field derived from the analysis in step (c).
6. Apply inverse NMO correction using the same velocity field as in step (d) and repeat the velocity analysis to derive a velocity field from the data which now have been first crossline migrated then 2-D prestack migrated in the inline direction. Figure 7.4-32 shows the common-cell gathers and the computed velocity spectra at four analysis locations along the same inline traverse as in Figure 7.4-31 after 2-D prestack time migration.
7. Apply NMO correction using the velocity field from step (f), stack the data and unmigrate using the same velocity field that was used in common-offset migration in step (e). The resulting sections shown in Figure 7.4-33 represent 2-D zero-offset wave-fields, and thus contain no sideswipe energy.
8. Use the velocity field from step (f) to remigrate the stacked sections. The resulting migrated sections shown in Figure 7.4-34 represent the final product from a 3-D prestack time migration based on the two-pass strategy outlined above.

Aside from being an integral component of a two-pass 3-D prestack time migration workflow described above, crossline migration has two other applications.

1. By restricting the aperture for Kirchhoff summation to crossline trace spacing, crossline migration can be used to perform azimuth moveout (AMO) correction of 3-D prestack data (processing of 3-D seismic data). The resulting prestack data still require DMO correction, but only as a 2-D process. As a bonus byproduct of this application, crossline migration yields prestack data with regular offset distribution which can be used as input to a 3-D prestack depth migration that requires uniformly sampled data (3-D prestack depth migration).
2. Crossline migration can be used to re-orient 3-D prestack data along a desired direction. This application is useful in merging data from a number of neighboring or partially overlapping 3-D surveys which have been conducted using different recording directions. Specifically, data from one survey can be crossline migrated to generate gathers along the inline direction associated with another survey.

### 3-D migration velocity analysis

Earlier in this section we learned that a workflow based on 3-D DMO correction combined with 3-D common-offset migration can include a step for estimating a 3-D rms velocity field from the migrated prestack data. Now, we shall develop an alternative workflow for 3-D migration velocity analysis based on a 3-D extension of the velocity-independent imaging technique by Fowler [55] described in migration velocity analysis.

1. Start with 3-D prestack data P(x, y, h, t) in coordinates of inline x, offset 2h, crossline y, and event time t in the unmigrated position, and apply NMO correction, 3-D DMO correction followed by inverse NMO correction.
2. Create a constant-velocity stack (CVS) volume P(x, y, t0; vDMO) using a constant velocity vDMO, where t0 is the zero-offset event time after 3-D DMO correction. Shown in Figure 7.4-35 is a CVS volume created by using a velocity of 2400 m/s. The prestack data are from the same survey as that of the data shown in Figure 7.4-33.
3. Migrate the constant-velocity stack (CVS) volume P(x, y, t0; vDMO) to create a constant-velocity migration (CVM) volume P(x, y, τ; vmig), where τ is the event time after migration, using the constant velocity associated with the CVS volume itself and a 3-D constant-velocity Stolt algorithm (Section G.4). Shown in Figure 7.4-36 is the CVM volume associated with the CVS volume of Figure 7.4-35 with a velocity of 2400 m/s.
4. Repeat steps (b) and (c) for a range of constant velocities. In the present example, a total of 50 CVS volumes and corresponding 50 CVM volumes were created using a velocity range of 1200-3650 m/s with an increment of 50 m/s.
5. Extract the time slices from the CVM volumes for a specified time τj and combine them to form a velocity cube P(x, y, vmig; τj). Figure 7.4-37a shows such a velocity cube for a time of 2000 ms. Repeat for a set of N times τj, j = 1, 2, …, N. In the present example, a total of 26 velocity cubes were created for times from 600 to 3100 ms with an increment of 100 ms.
6. For 3-D visualization and interpretation, form a super volume that is a composite of the 26 velocity cubes from step (e) placed one on top of the other (Figure 7.4-38).
7. Based on the maximum-event-amplitude picking strategy described in migration velocity analysis, interpret each of the velocity cubes and create velocity strands associated with the constant times τj, j = 1, 2, …, N, along selected crossline (or inline) traverses (Figure 7.4-37b).
8. Combine the velocity strands for a specific time τj and create an rms velocity map associated with that time. Repeat for all times τj, j = 1, 2, …, N, and obtain a set of rms velocity maps (Figure 7.4-39).
9. Create a 3-D rms velocity field using the rms velocity maps from step (h). Figures 7.4-40 and 7.4-41 show selected inline and crossline sections from the rms velocity volume. Compare with the corresponding image sections from time migration shown in Figure 7.4-34, and note that the 3-D rms velocity field honors the structural characteristics of the subsurface. Specifically, the inline sections in Figure 7.4-40 manifest the basin boundary indicated by the change in color from green (the basin interior) to orange-red. The crossline sections in Figure 7.4-41 manifest the presence of young sediments with low velocities (the shallow purple zone) and older sediments with relatively higher velocities (the green) disrupted by an erosional channel between 1-2 s (the dark green).

The 3-D rms velocity field derived from the workflow outlined above is associated with events in their migrated positions. Compared to the 3-D DMO velocity field associated with events in their yet unmigrated positions, it is the preferred velocity field with which you would want to migrate your data using a desired 3-D post- or prestack time migration algorithm (Appendix G). If you wish to extend your analysis from time to depth domain, the 3-D rms velocity field associated with time-migrated data also is the preferred velocity field to derive a 3-D interval velocity field that can then be used to obtain an image in depth (3-D structural inversion applied to seismic data from the northeast China).

### Aspects of 3-D prestack time migration — a summary

The improvement in imaging with 3-D prestack time migration may sometimes be marginal compared to imaging with 3-D poststack time migration of 3-D DMO-stacked volume of data. Nevertheless, the benefits of 3-D prestack time migration are not limited to the improved image of the subsurface as listed below.

1. 3-D prestack time migration is the appropriate strategy for imaging conflicting dips with different stacking velocities, such as reflections from steeply dipping fault planes with 3-D geometry and gently dipping strata.
2. The 3-D image volume derived from 3-D prestack time migration is used as an input to 3-D zero-offset amplitude inversion to estimate an acoustic impedance model of the earth (acoustic impedance estimation).
3. The CRP gathers from 3-D prestack time migration are used to perform prestack amplitude inversion to derive amplitude variation with offset (AVO) attributes (analysis of amplitude variation with offset).
4. The 3-D image volume can be inverse migrated by way of 3-D zero-offset wavefield modeling to derive an unmigrated 3-D zero-offset data volume.
5. The 3-D rms velocity field estimated from the 3-D prestack time-migrated data, as from step (f) of the workflow described above, can be used to derive a 3-D interval velocity field by Dix conversion.
6. Finally, the 3-D interval velocity field and the 3-D zero-offset wavefield can be used to derive an earth image in depth by 3-D poststack depth migration. The earth image volume in depth can then be interpreted to delineate a set of reflector geometries associated with key geological markers. The interval velocity field combined with the reflector geometries may be used to build an initial earth model in depth (earth modeling in depth).

## 7.5 Interpretation of 3-D seismic data

After 3-D migration, a 3-D data volume is suitable for derivation of the 3-D subsurface geological model. Because of the completeness of data within that volume, the interpreter has more information than he would from the analysis of a 2-D seismic data set. Although more information leads to less uncertainty in deriving the geological model, interpretation of abundant data can be exhausting. The use of an interactive interpretation workstation is a common way to deal with large 3-D data volumes. Moreover, an interactive environment is versatile in viewing the 3-D data volume. For example, we can examine vertical sections in inline, crossline, or any arbitrary direction, as well as horizontal sections, namely, time slices. An interactive environment also can provide the capability to improve interpretation. For example, horizon flattening, correlation of marker horizons across faults, and some image processing tools to enhance certain features within the data volume can be implemented.

First-generation 3-D interpretation systems were based on the principle of line-by-line interpretation. Contemporary 3-D interpretation systems are based on volume visualization and manipulation of amplitudes within the image volume to enhance structural and stratigraphic features of interest. With the availability of powerful hardware and software visualization tools, interpretation no longer is picking a set of horizons and discarding the rest of the image volume. Instead, interpretation involves using both traveltime and amplitude information contained in the entire image volume. Specifically, structural interpretation is primarily based on picking traveltimes that are coincident with geological layer boundaries, and stratigraphic interpretation is based on manipulation of seismic amplitudes to enhance subtle features associated with depositional environment and sedimentology.

### Time slices

A time slice contains events from more than one reflection horizon at the same time level. A spatially high-frequency event on a time slices is either a steeply dipping event or a high-frequency event in time. Thus, from the time slices in Figure 7.5-1, we can infer steep dip at location H from the high-frequency character of the event, and gentle dip at location L from the low-frequency character. Additionally, contours associated with a reflection horizon can be traced from time slices. If contours narrow between a shallow and a deeper time slice, then the feature is a structural low (Figure 7.5-1). Conversely, if contours widen between a shallow and a deeper time slice, then the feature is a structural high (Figure 7.5-2).

Time slices can be used to generate structural contour maps. Some 2-D smoothing can be applied to time slices to ease contouring (Figure 7.5-3b). Edge enhancement can be used to better delineate zero crossings (Figure 7.5-3c). Some image processing tools also can be used to detect and enhance subtle structural features on time slices. Illumination of a time slice by a pseudosun is an example of such a technique (Figure 7.5-3d). The angle of illumination and its direction can be manipulated for optimum detection. The main structure in Figure 7.5-3d is a salt dome. Note the presence of tensional faults associated with salt diapirism. Bone [56], Brown [57], and Brown [58] have made pioneering use of of time slices in 3-D interpretation.

One common artifact that is observable on time slices is the presence of horizontal striations. The pseudo-sun illumination has enhanced the streaks in Figure 7.5-3d. These striations typically are oriented along the shooting direction. Although there are several causes for time-slice striations, one cause is the positioning errors present in the navigation data. Schultz and Lau [59] identify these striations as time shifts from one inline to another (crossline statics). They also suggest a poststack procedure to estimate and eliminate these time shifts in the inline-crossline wavenumber domain.

Besides contouring, time slices also are useful in quality control. Figure 7.5-4 is a time slice before and after residual statics corrections were applied to the data. Note the improvement in the signal-to-noise ratio and the continuity of events in Figure 7.5-4b.

Time slices can be used to check for consistency in picking time horizons along inlines and crosslines. Actually, in some areas with very complex structures, time slices can sometimes be used to trace faults and horizon contours. Figure 7.5-5a shows a time slice that exhibits a major fault from left to right across the survey area and a series of oblique fault blocks adjoining the major fault zone. A time slice contributes one contour level for a horizon as illustrated in Figure 7.5-5b. Such a display can be used for the quality control of structural interpretation.

### 3-D visualization

Picking horizons and faults and then abandoning the 3-D image volume no longer is practiced in interpreting 3-D seismic data. Horizon picking amounts to interpreting seismic traveltimes, and as such, it yields the knowledge for structural geology, only. The knowledge for stratigraphy, depositional environment and history, and even lithology can only be inferred by the combined interpretation of traveltimes and amplitudes contained in the seismic image volume.

3-D visualization facilitates combined interpretation of multiple volumes of data — image volumes, velocity volumes, and attribute volumes such as those derived from amplitude variation with offset (AVO) analysis and acoustic impedance estimation. Rapid detection of amplitude anomalies, and understanding and correlating the geometries of the various layer boundaries and faults are made possible by the power of 3-D visualization. Combining visualization with interpretation tools such as opacity removal, seed detection, and edge detection enable the interpreter to get the most out of the data volumes.

3-D visualization of seismic amplitudes is based on the concept of representing each sample in the data volume by a 3-D object called voxel. The voxel representation actually is an extension of a 2-D representation of a sample by a pixel. Each voxel is color coded by the amplitude of the sample associated with it.

Figure 7.5-6a shows a single depositional unit carved out from an image volume. Note the distinctive pattern of the basin-edge faults. Representation of these faults as a set of surfaces in the absence of the depositional unit would only provide a limited understanding of the structural geology of the subsurface. Figure 7.5-6b shows an example of enhancing structural features such as the intensive fault pattern over the survey area by introducing opacity to the time slab from the image volume. The fault pattern also can be seen on a map view as in Figure 7.5-6c.

### Removal of opacity

Figure 7.5-7a shows a horizontal slab of about 100 ms in thickness isolated from a 3-D poststack time-migrated volume of data. The amplitudes within the slab were actually manipulated to enhance the structural and de-positional features. Amplitude manipulation involves suppressing a range of positive and negative amplitudes below a specified threshold value so as to remove the opacity effect caused by the suppressed amplitudes. The same slab with labeled depositional features is shown in Figure 7.5-7b. Among the features identified are a carbonate shelf, the shelf edge, a channel system, and two submarine fan lobes. The color indicates the seismic amplitudes retained in the slab. Figure 7.5-7b is a product of what interpretive work means — a description of the depositional environment.

Figure 7.5-2  Part 1: Selected time slices from a land 3-D survey starting at 1000 ms and moving down to 1660 ms at 60-ms intervals. (Data courtesy Nederlandse Aardolie Maatschappij B.V.)

To demonstrate that interpretation is not just picking a set of time horizons, but it is more what you do after picking, consider the data set in Figure 7.5-8. We see a sample time slice (Figure 7.5-8a), a crossline (Figure 7.5-8b) and an inline (Figure 7.5-8c). The objective is to identify sand dunes which are confined between horizons T6 and T8. None of the conventional displays in the form of cross-sections of the image volume as shown in Figure 7.5-8 can provide a rapid identification of the sand dunes.

Figure 7.5-2  Part 2: Selected time slices from a land 3-D survey starting at 1000 ms and moving down to 1660 ms at 60-ms intervals. (Data courtesy Nederlandse Aardolie Maatschappij B.V.)

Start with the entire image volume and isolate a time slab bounded by the time horizons T6 and T8 as labeled in Figures 7.5-8b,c. Then, remove opacity to enhance the series of sand dunes as shown in Figure 7.5-9a. A vertical cross-section of the slab with the individual sand dunes is shown in Figure 7.5-9b.

An interpretation of a deltaic sequence is made possible by 3-D visualization and amplitude manipulation. Figure 7.5-10a shows the sequence under consideration highlighted on a vertical cross-section from a time-migrated volume of data. Figure 7.5-10b shows the isolated deltaic sequence with the individual members of the sequence identified by removing opacity from the subvolume that includes the deltaic sequence highlighted in Figure 7.5-10a.

### Seed detection

Figure 7.5-11a shows a horizontal slab isolated from a 3-D poststack time-migrated volume of data. The objective is to isolate a bright spot located within the slab. First, remove opacity to obtain the image shown in Figure 7.5-11b. Note that the bright spot now can be easily detected. Next, label one voxel within the interior of the bright spot as if you are implanting a seed. Then, search for all neighboring voxels that represent amplitudes within a specified range of the seed amplitude. This process results in isolating a set of voxels that are contiguous in the form of a subvolume as highlighted in Figure 7.5-11c. The detected subvolume associated with the bright spot is displayed in this figure within the context of the surrounding depositional environment.

Figure 7.5-3  Processing time-slice data for improved interpretation: (a) A time slice, (b) after 2-D smoothing, and (c) edge enhancement; (d) another time slice with pseudo-sun illumination applied to enhance the intensive fracturing associated with salt diapirism. (Data in (a) courtesy Shell Oil Company and Esso, and data in (d) BP Amoco.)

The seed detection can be repeated by adjusting the threshold amplitude associated with the seed voxel to produce a series of subvolumes with increasing size as shown in Figures 7.5-12a through f. It is evident that the case shown in Figure 7.5-12b best defines the size of the bright spot. The subvolume that represents the bright spot can be isolated completely from the entire image volume as shown in Figure 7.5-12g. By proper color coding of amplitudes within the subvolume, the potential hydrocarbon-bearing feature labeled as sweet sand can be identified. Also, observe the improvised structural contours around the sand body.

Figure 7.5-4  A time slice (a) before and (b) after residual statics corrections. (Data courtesy Nederlandse Aardolie Maatschappij B.V.)

Seed detection can be used to interpret horizons with complex geometry without actually picking traveltimes. Based on a crude definition of the horizon geometry, a horizon-consistent time slab is extracted from the data volume used in interpretation. Figure 7.5-13a shows a color-coded time horizon associated with the top boundary of a karstic limestone layer interpreted from a time-migrated volume of data using the seed detection technique applied to the time slab. The same horizon in map view is shown Figure 7.5-13b.

### Structural interpretation

A 3-D structural interpretation session may begin with viewing selected inline and crossline sections to acquire a regional understanding of the subsurface geology. Other orientations, such as vertical sections along a dominant dip direction, also may be needed to determine the structural pattern. Time slices then are studied to check the structural pattern. These previews may be made dynamic in an interactive environment; vertical or horizontal sections can be viewed in rapid succession as one would a film strip from a motion picture. Any change in structure in space and time can thus be grasped with ease.

Figure 7.5-14 shows an image volume derived from 3-D prestack time migration of data from a marine 3-D survey. There are 296 inlines and 1300 crosslines in the data volume. Figures 7.5-15 and 7.5-16 show snapshots of the image volume as it is being visualized in the inline and crossline directions, respectively. The time slices shown in Figure 7.5-17 also represent snapshots from the 3-D visualization session. Note the complex structural pattern represented by a principal fault across the survey area from left to right and a series of fault blocks on both sides of this fault.

The water bottom and six deeper time horizons were picked using a combination of seed detection and line-based interpretation.

1. First, the seismic event for each horizon was identified within the image volume. Then, at locations with good continuity and signal-to-noise ratio, seed points were placed at the seismic event. For consistency, all seed points for a given horizon were assigned identically — either at trough or peak. By using the seed detection method, the event was tracked away from each seed point as far as possible laterally in all directions. The result of seed detection is a series of surface patches associated with the horizon under consideration (Figure 7.5-18). Figure 7.5-19 shows the surface patches for the five horizons for which seed detection was applied. The water-bottom horizon is not shown since it is almost flat. The seed detection was not suitable for the deepest horizon H6. Note that for horizons H1-H5, seed detection was largely successful away from the fault zones.
2. Where seed detection failed, control points were picked along a grid of selected inlines and crosslines (Figure 7.5-20). In fact, the deepest horizon H6 was picked entirely using the line-based interpretation. Figure 7.5-21 shows the surface patches derived from seed detection and line-based picks for all six horizons.

The combined result of seed-based and line-based interpretation of all six horizons is shown in Figure 7.5-22. Horizons H1-H6 of Figure 7.5-22 are ordered from top to bottom. Specifically, interpretation has produced a set of picks for each horizon. These picks are used in grid calculations to create complete surfaces for all the horizons over the entire survey area. Figure 7.5-23 shows the gridded surface associated with horizon H3 using the picks from seed-based and line-based interpretations shown in Figure 7.5-20. The gridded surfaces for all six horizons based on the picks shown in Figure 7.5-21 are shown individually in Figure 7.5-24. The 3-D view that exhibits the final result of picking shown in Figure 7.5-22 is repeated for the final result of gridding as shown in Figure 7.5-25.

During and after an interpretation session, it is imperative to check the consistency of the picked horizons with the seismic data volume that was used in the interpretation. Figure 7.5-26 shows combined views of the time horizons and seismic cross-sections from the image volume as in Figure 7.5-18. The final step in structural interpretation is time-to-depth conversion of time horizons interpreted from the time-migrated volume of data. This step requires knowledge of interval velocities for each layer, and as such, it implies earth modeling in depth. We shall therefore refer to the discussion on time-to-depth conversion and calibration of the results of depth conversion to well data in earth modeling in depth.

A structural interpretation session normally would include explicit delineation of fault framework. Nevertheless, the objective in the present case study is the investigation of stratigraphic features within the reservoir zone, and as such, the time horizons shown in Figure 7.5-22 were interpreted without the explicit definition of the fault patterns. Instead, the fault surfaces were interpreted as part of the horizons themselves.

### Stratigraphic interpretation

A session for 3-D stratigraphic interpretation may begin with splitting the image volume as in Figure 7.5-18 into subvolumes that correspond to individual depositional units bounded by the time horizons derived from structural interpretation. Figure 7.5-27 shows the image volume stripped down to each time horizon as in Figure 7.5-22, while Figure 7.5-28 shows the subvolumes that represent the depositional unit for each layer.

The stratigraphic interpretation involves removal of opacity and seed detection.

1. Removal of opacity is applied to either a horizontal time slab with a specified thickness, typically a few to tens of time samples, or to a depositional unit bounded by the time horizons derived from structural interpretation. Figure 7.5-29 shows opacity removal applied to a thin horizontal slab that includes the water bottom, the horizontal slab slightly deeper than the water bottom, and the depositional unit labeled as H1 in Figure 7.5-28. Note the enhanced images of a complex channel system at the water bottom, and the intensive fracture system that begins to develop immediately below the water bottom and increases in complexity as we go deeper in the image volume. When we reach horizon H3 as labeled in Figure 7.5-28, we observe a highly complex fault system (Figure 7.5-30).
2. Now we look inside a specific depositional unit, in this case the unit bounded by horizon H3 on top as labeled in Figure 7.5-28. The subvolume defined by a horizon-consistent slab associated with this unit is shown also in Figure 7.5-31a. Examine selected vertical cross-sections from this slab as shown in Figure 7.5-31b and identify stratigraphic features of interest. In this case, we note the presence of a buried channel represented by a strong reflection event with a lateral extent up to 400 m (Figure 7.5-31c). Remove the opacity within the subvolume under consideration and uncover a remnant of an ancient channel system as shown in Figure 7.5-32. Note that the fault system, which has been formed by the subsequent tectonic disturbance, disrupts the channel system at several locations.
3. Apply seed detection to the event that represents the channel within the subvolume and isolate it completely from the rest of the subvolume. Shown in Figure 7.5-33 is the channel with the surrounding flood plain represented by the white speckles. The color within the channel represents the seismic amplitudes. Also shown in Figure 7.5-33 is horizon H4 which defines the base of the reservoir unit that contains the channel. A close-up view of the structural and stratigraphic interpretation result in Figure 7.5-33 is shown in Figure 7.5-34a.

As for structural interpretation, the result of stratigraphic interpretation needs to be checked for consistency with the image volume of data used in the interpretation. Figure 7.5-34b shows a views of the channel, horizon H4, and an inline section from the image volume derive from 3-D prestack time migration. The color within the channel in this view represents the elevation; note that it varies considerably because of the disturbance of the channel by the subsequent faulting.

## Exercises

Exercise 7-1. From the time slices in Figure 7.E-1, infer the structural pattern in the subsurface.

Exercise 7-2. Approximately at what angle is the pseudosun illuminating the time slice in Figure 7.5-3?

Exercise 7-3. What are the advantages and disadvantages of dip-line and strike-line shooting in terms of (a) crossline smear, (b) DMO, (c) spatial aliasing, and (d) velocity estimation?

Exercise 7-4. Suppose you did a 2-D survey over the salt dome structure shown in Figure 7.3-9 and applied 2-D imaging to your data. Then you did a 3-D survey and 3-D imaging. Consider mapping the top of the salt. Which one, 2-D or 3-D seismic, will imply a salt dome structure with a larger spatial extent and structural closure?

Exercise 7-5. If you use the 90-degree 3-D phase-shift algorithm (Section G.3) in 3-D time migration, do you have to do it in one pass?

Exercise 7-6. Consider a dual-boat operation in a 3-D survey such that line spacings alternate between 112.5 and 37.5 m. Is this equivalent to a single-boat operation with 75-m line spacing?

Exercise 7-7. How would you obtain 3-D migration of a crossline section without having to do 3-D migration of the entire 3-D data set?

## Appendix G: Mathematical foundation of 3-D migration

### G.1 Implicit methods

The 2-D scalar wave equation (D-1) given in Section D.1 can be adapted to three dimensions as

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

This equation describes propagation of a 3-D compressional zero-offset wavefield P(x, y, z, t) in a medium with constant material density and compressional wave velocity v(x, y, z), where x is the horizontal spatial axis in the inline direction, y is the horizontal spatial axis in the crossline direction, z is the depth axis (positive downward), and t is time. Given the upcoming seismic wavefield P(x, y, z = 0, t), which is recorded at the surface, we want to determine reflectivity P(x, y, z, t = 0). This requires extrapolating the surface wavefield to depth z, then collecting it at t = 0.

The solution to the 2-D scalar wave equation (D-1) is given by equation (D-7), which can be adapted to extrapolate a 3-D zero-offset wavefield in depth

 ${\displaystyle P(k_{x},k_{y},z,\omega )=P(k_{x},k_{y},z,\omega )\ \ \exp(-ik_{z}z),}$ (G-2)

where

 ${\displaystyle k_{z}={\frac {2\omega }{v}}{\sqrt {1-X^{2}-Y^{2}}},}$ (G-3)

with

 ${\displaystyle X={\frac {vk_{x}}{2\omega }},}$ (G-4a)

and

 ${\displaystyle Y={\frac {vk_{y}}{2\omega }}.}$ (G-4b)

In equations (G-4a,G-4b), kx and ky are the wavenumbers in the inline and crossline directions, respectively, and ω is the temporal frequency in units of radians per unit time. Note that the inline and crossline normalized wavenumbers X and Y are coupled in equation (G-3).

First, consider the constant-velocity case and rewrite the dispersion relation of equation (G-3) as for the 2-D case (equation D-49a) to get

 ${\displaystyle \omega _{\tau }={\sqrt {\omega ^{2}-{\frac {v^{2}k_{x}^{2}}{4}}-{\frac {v^{2}k_{y}^{2}}{4}}}},}$ (G-5)

where ωτ = vkz/2 is the output temporal frequency.

Suppose we perform 2-D migration on all the inlines by mapping the input frequency ω to an output frequency ω1 using the relation

 ${\displaystyle \omega _{1}={\sqrt {\omega ^{2}-{\frac {v^{2}k_{x}^{2}}{4}}}}.}$ (G-6)

We then sort the already migrated inlines in the crossline direction. Next, suppose that we perform 2-D migration on all the crosslines. This second-pass migration is performed by mapping ω1 of equation (G-5) to an output frequency ω2 using the relation

 ${\displaystyle \omega _{2}={\sqrt {\omega _{1}^{2}-{\frac {v^{2}k_{y}^{2}}{4}}}}.}$ (G-7)

When equation (G-7) is substituted into equation (G-6), we obtain

 ${\displaystyle \omega _{2}={\sqrt {\omega ^{2}-{\frac {v^{2}k_{x}^{2}}{4}}-{\frac {v^{2}k_{y}^{2}}{4}}}},}$ (G-8)

which is the same as equation (G-5). This means that 3-D migration can be performed in two passes: 2-D migration in the inline direction followed by 2-D migration in the crossline direction. The above derivation of a 90-degree 3-D migration operator [27] is based on the constant-velocity assumption.

Mathematically, the two-pass approach is based on the idea of full separation [28] of migration operators in the inline and crossline directions. The separation is strictly valid for the constant-velocity medium. When velocities vary spatially, splitting [28] of the inline and crossline operators is necessary (Figure 7.3-1). A rigorous mathematical treatment of separation and splitting of extrapolation operators can be found in a paper by Brown [29].

To implement separation or splitting, the inline and crossline wavenumbers need to be decoupled in the 3-D dispersion relation (equation G-3). By Taylor series expansion of the square root and by keeping the first three terms, we get the 15-degree dispersion relation (Section D.3)

 ${\displaystyle k_{z}\approx {\frac {2\omega }{v}}\left(1-{\frac {X^{2}}{2}}-{\frac {Y^{2}}{2}}\right).}$ (G-9)

When the cross-dip component is zero (Y = 0), equation (G-9) reduces to its 2-D form (equation D-54). The inline and crossline terms in equation (G-9) are separable when velocity is slowly variable in z or independent of x and y [29].

The impulse response of the 2-D 15-degree migration operator with its dispersion relation given by equation (D-42a) is an ellipse (Figure 4.3-1). The impulse response of the 3-D 15-degree migration operator with its dispersion relation given by equation (G-9) is an ellipsoid (Figure 7.3-2). Except for the amplitudes, the vertical cross-section at center line (51) is the same as the 2-D impulse response. Within typical bandwidth of seismic data, the 3-D 15-degree impulse response departs from the ideal 3-D impulse response (a hollow hemisphere) at dips greater than approximately 30 degrees. The time slices of the impulse response for equation (G-9) are circular in shape.

Additional problems arise when more accuracy is desired; the X and Y terms in equation (G-3) no longer are decoupled. For example, the 45-degree approximation (Section D.4) to equation (G-3) introduces cross terms. Another approximation is to represent the single square root in equation (G-3) by two square roots [41] as

 ${\displaystyle k_{z}\approx {\frac {2\omega }{v}}\left[{\sqrt {1-X^{2}}}+{\sqrt {1-Y^{2}}}-1\right].}$ (G-10)

Note that equation (G-10) has the same form as equation (D-15). By Taylor series expansion of the two square roots, equation (G-10) reduces to equation (G-9). Assuming a small cross-dip component, more accuracy in the inline direction is achieved by making the 45-degree approximation given by

 ${\displaystyle k_{z}\approx {\frac {2\omega }{v}}\left(1-{\frac {X^{2}}{2-{\frac {X^{2}}{2}}}}-{\frac {Y^{2}}{2}}\right).}$ (G-11)

When more accuracy is required in both directions, we can apply the 45-degree expansion to equation (G-10) in both the inline and crossline directions given by

 ${\displaystyle k_{z}\approx {\frac {2\omega }{v}}\left(1-{\frac {X^{2}}{2-{\frac {X^{2}}{2}}}}-{\frac {Y^{2}}{2-{\frac {Y^{2}}{2}}}}\right).}$ (G-12)
Figure G-1  An algorithmic description of one-pass implicit 3-D poststack migration based on 3-D wave extrapolation and imaging in the frequency-space domain (Section G.1).

The impulse response of the 2-D 45-degree migration operator with its dispersion relation given by equation (D-58) is heart shaped (Figure 4.3-1). The impulse response of the 3-D 45-degree migration operator with its dispersion relation given by equation (G-12) is shown in Figure 7.3-3. Except for the amplitudes, the vertical cross-section at center line (51) is the same as the 2-D impulse response. Within the typical bandwidth of seismic data, the 3-D 45-degree impulse response departs from the ideal 3-D impulse response (a hollow hemisphere) at dips greater than approximately 65 degrees. The time slices of the impulse response for equation (G-12) are diamond shaped. Note that the response is accurate in both the inline and crossline directions, and that the largest error in the form of undermigration occurs in the two diagonal directions.

In summary, to implement separation or splitting, the inline and crossline terms of the 3-D migration operator need to be decoupled. This can be achieved by making rational approximations (equation G-9 or G-12) to the square root in equation (G-3). Separation works for 90-degree, constant-velocity (i.e., Stolt migration) or 15-degree, slowly varying v(z) cases. When there are large vertical velocity gradients or strong lateral velocity variations, splitting must be used. The algorithm for a one-pass implicit frequency-space 3-D poststack migration that uses the splitting method is identical to the 2-D algorithm outlined in Figure 4.1-23, except for the additional application of the diffraction term in the crossline direction (Figure G-1).

### G.2 Explicit methods

Rewrite equation (G-2) for extrapolating a 3-D zero-offset wavefield in discrete depth steps of Δz

 ${\displaystyle P(k_{x},k_{y},z+\Delta z,\omega )=P(k_{x},k_{y},z,\omega )\ \ \exp(-ik_{z}\Delta z).}$ (G-13a)

Combine equations (G-3) and (G-4) to express the vertical wavenumber kz in terms of the other transform variables as

 ${\displaystyle k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk_{x}}{2\omega }}\right)^{2}-\left({\frac {vk_{y}}{2\omega }}\right)^{2}}},}$ (G-13b)

where, v is the extrapolation velocity.

Substitute

 ${\displaystyle k={\sqrt {k_{x}^{2}+k_{y}^{2}}}}$ (G-14)

into equation (G-13b), and obtain

 ${\displaystyle k_{z}={\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk}{2\omega }}\right)^{2}}}.}$ (G-15)

With the vertical wavenumber kz given by equation (G-15), we now write the desired extrapolation operator D(k) given by the exponential term in equation (G-13a) for a specific ratio of 2ω/v as

 ${\displaystyle D(k)=\exp \left[-i{\frac {2\omega }{v}}{\sqrt {1-\left({\frac {vk}{2\omega }}\right)^{2}}}\Delta z\right].}$ (G-16)

This operator can be, in principle, inverse transformed back to the frequency-space domain, and thus, extrapolation described by equation (G-13) can be achieved by convolution of the wavefield with the extrapolation filter in the frequency-space domain.

Here, we describe a practical scheme [33],[34] which involves the design of a stable 2-D extrapolation filter, use of the McClellan transform to obtain a 3-D extrapolation operator without actually computing a 3-D operator, and replacing convolution for filter application by an efficient recursive scheme based on Chebychev polynomials.

Start with the extrapolation operator in the Fourier transform domain given by equation (G-16). To accommodate lateral velocity variations, we want to compute an extrapolation filter hn in the frequency-space (ω, x, y) domain with a Fourier transform ${\displaystyle H(k={\sqrt {k_{x}^{2}+k_{y}^{2}}},\omega ,v)}$ that best approximates the desired transform given by equation (G-16). Since D(k) is symmetric with respect to k = 0, the complex filter coefficients hn are even — h−n = hn, and the number of filter coefficients 2N + 1 is odd. As a result, we only need to compute N coefficients. The actual Fourier transform of hn is expressed as

 ${\displaystyle H(k)=h_{0}+2\sum _{n=1}^{N}h_{n}\cos(nk).}$ (G-17)

Methods to obtain hn such that the actual Fourier transform H(k) given by equation (G-17) closely approximates the actual Fourier transform D(k) given by equation (G-16) is described by Holberg [32] and Hale [33]. The former is based on least-squares and the latter is based on a modified Taylor series. The Hale method is described in Section D.5. Since the desired amplitude spectrum is |D(k)| = 1 for all k, then we require the actual amplitude spectrum to be |H(k)| ≤ 1 for all k. This will ensure that the operator is stable — that is, amplitudes of the extrapolated wavefield will not grow from one depth level to another.

By using a special case of the recursive formula for the Chebychev polynomials [34], equation (G-17) can be rewritten in terms of cos k as

 ${\displaystyle H(k)=h_{0}+2h_{1}\cos k+2h_{2}(2\cos k\cos k-1)+\cdots ,}$ (G-18)

so that H(k) is expressed entirely in terms of hn and cos k. As a result, convolution with the filter hn can be performed by a recursive Chebychev filter structure described by Hale [34]. This is essential for an efficient application of a 2-D explicit extrapolation filter, since direct convolution in two dimensions is very costly.

 1/8 1/4 1/8 1/4 −1/2 1/4 1/8 1/4 1/8

For the design of the 3-D explicit operator, define the transform G(k) = cos k and note from equation (G-14) that

 ${\displaystyle G(k_{x},k_{y})=\cos {\sqrt {k_{x}^{2}+k_{y}^{2}}}.}$ (G-19)

For computational efficiency, G(kx, ky) given by equation (G-19) may be approximated as [34]

 ${\displaystyle G(k_{x},k_{y})\approx -1+{\frac {1}{2}}(1+\cos k_{x})(1+\cos k_{y}).}$ (G-20)

This approximation is exact for kx = 0 and ky = 0, best for small kx and ky, and deteriorates at large kx and ky.

G(kx, ky) of equation (G-20) is the Fourier transform of the 2-D McClellan transformation filter given by Table G-1 in the x − y domain, independent of ω and v [34].

Using the McClellan filter coefficients given by Table G-1, the 1-D extrapolation filter hn is transformed into a 2-D extrapolation filter with approximate circular symmetry [34]. Again, by way of a Chebychev filter structure, application of this filter is done more efficiently as compared with direct convolution. The impulse response of the 3-D explicit migration operator that uses the 3 × 3 McClellan transform template of Table G-1 is shown in Figure 7.3-13.

The circular symmetry of the McClellan transform filter can be improved by an alternative approximation to G(kx, ky) of equation (G-19) by [34]

 ${\displaystyle G(k_{x},k_{y})=-1+{\frac {1}{2}}(1+\cos k_{x})(1+\cos k_{y})-{\frac {c}{2}}(1-\cos 2k_{x})(1-\cos 2k_{y}),}$ (G-21)

where c is an adjustable scalar to make the Fourier transform of the 2-D filter based on the McClellan transform match the desired Fourier transform exactly (equation G-16) at some ${\displaystyle k={\sqrt {k_{x}^{2}+k_{y}^{2}}}.}$ A typical value for c = 0.0255. G(kx, ky) of equation (G-21) is the Fourier transform of the 2-D McClellan transformation filter given by Table G-2 in the x − y domain, again, independent of ω and v.

 −c/8 0 c/4 0 −c/8 0 1/8 1/4 1/8 0 c/4 1/4 −(1 + c)/2 1/4 c/4 0 1/8 1/4 1/8 0 −c/8 0 c/4 0 −c/8

The impulse response of the 3-D explicit migration operator that uses the 5 × 5 McClellan transform template of Table G-2 is shown in Figure 7.3-14.

Aside from the application of the McClellan transform, the circular symmetry of the 3-D impulse response can also be attained by applying the correction term by Li [30] to the extrapolated wavefield, intermittently, at some depth steps. The objective is to compensate for the accumulated errors caused by the approximations made to the exact 3-D extrapolation operator given by equation (G-16). Li [30] devised a correction term originally to compensate for the errors caused by the 45-degree split form of the dispersion relation given by equation (G-12). The error incurred during wave extrapolation using an implicit scheme arises from the difference between the exact form of the dispersion relation given by equation (G-13b) and an approximate form such as that given by equation (G-12).

Etgen and Nichols [31] later adapted Li’s correction to compensate for the errors caused by an approximation made to the explicit operator given by equation (G-16). The error incurred during wave extrapolation using an explicit scheme can be attributed to the difference between the exact form of the wavenumber k given by equation (G-14) and an approximate form ${\displaystyle {\tilde {k}}}$ such as that given by equation (G-19)

 ${\displaystyle {\tilde {k}}=\cos ^{-1}G(k_{x},k_{y}),}$ (G-22)

with G(kx, ky) given by the approximate form in equation (G-21).

Write equation (G-16) for the approximate form of the wavenumber k given by equation (G-22) as

 ${\displaystyle {\tilde {D}}(k)=\exp \left[-i{\frac {2\omega }{v}}{\sqrt {1-\left({\frac {v{\tilde {k}}}{2\omega }}\right)^{2}}}\Delta z\right].}$ (G-23)

By using the exact (equation G-16) and the approximate (equation G-23) forms of the extrapolation operators, the error in the extrapolation accumulated after nΔz depth intervals is then given by

 ${\displaystyle E(k_{x},k_{y},\omega )=\exp \left\{-i{\frac {2\omega }{v}}\left[{\sqrt {1-\left({\frac {vk}{2\omega }}\right)^{2}}}-{\sqrt {1-\left({\frac {v{\tilde {k}}}{2\omega }}\right)^{2}}}\right]n\Delta z\right\}.}$ (G-24)

Therefore, the extrapolated wavefield P(kx, ky, nΔz, ω) given by equation (G-13a) is compensated for the errors incurred by the approximate form of G(kx, ky) given by equation (G-21) by applying the conjugate of E(kx, ky, ω) of equation (G-24) to the extrapolated wavefield

 ${\displaystyle P(k_{x},k_{y},n\Delta z,\omega )=E^{\ast }(k_{x},k_{y},\omega )P(k_{x},k_{y},n\Delta z,\omega ).}$ (G-25)

Following Li’s correction adapted to the explicit schemes [31] as given by equation (G-25), the wave extrapolation is continued with the approximate form (equation G-23) until more error is accumulated.

### G.3 3-D phase-shift migration

We start with the solution of the 3-D scalar wave equation to extrapolate a 3-D zero-offset wavefield P(x, y, z = 0, t) as given by equation (G-2) in the frequency-wavenumber domain. We also assume a horizontally layered earth model associated with a vertically varying velocity function v(z). By inverse Fourier transforming equation (G-2), we have

 ${\displaystyle P(x,y,z,t)=\int \int P(k_{x},k_{y},0,\omega )\ \exp(-ik_{z}z)\ \exp(-ik_{x}x-ik_{y}y+i\omega t)\ dk_{x}\ dk_{y}\ d\omega ,}$ (G-26)

where kz is defined by equation (G-13) adapted to the exploding reflectors model by replacing v with v/2.

The imaging principle t = 0 then is applied to equation (G-26) to get the migrated section P(x, y, z, t = 0) as

 ${\displaystyle P(x,y,z,t=0)=\int \int \int P(k_{x},k_{y},0,\omega )\ \exp(-ik_{x}x-ik_{y}y-ik_{z}z)\ dk_{x}\ dk_{y}\ d\omega ,}$ (G-27)

This is the equation for the phase-shift method [60] adapted to 3-D migration. Equation (G-27) involves integration over frequency and inverse Fourier transformation along inline and crossline axes at each depth step.

### G.4 3-D stolt migration

We now consider the special case of constant velocity v. As described in Section D.7, Stolt [61] devised a migration technique that involves an efficient mapping in the Fourier transform domain from temporal frequency ω to vertical wavenumber kz. We rewrite equation (G-13) to get an explicit expression for ω as

 ${\displaystyle \omega ={\frac {v}{2}}{\sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}.}$ (G-28)

By keeping the horizontal wavenumbers kx and ky unchanged and differentiating, we get

 ${\displaystyle d\omega ={\frac {v}{2}}{\frac {k_{z}}{\sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}}dk_{z}.}$ (G-29)

When equations (G-28) and (G-29) are substituted into equation (G-27), we obtain

 ${\displaystyle P(x,y,z,t=0)=\int \int \int \left[{\frac {v}{2}}{\frac {k_{z}}{\sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}}\right]\times P\left[k_{x},0,{\frac {v}{2}}{\sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}\right]\exp(-ik_{x}x-ik_{y}y-ik_{z}z)\ dk_{x}\ dk_{y}\ dk_{z}.}$ (G-30)

This is the equation for constant-velocity 3-D Stolt migration. It involves two operations in the f − k domain. First, the temporal frequency ω is mapped onto the vertical wavenumber kz via equation (G-28). Second, the amplitudes are scaled by the quantity S:

 ${\displaystyle S={\frac {v}{2}}{\frac {k_{z}}{\sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}}},}$ (G-31)

which is equivalent to the obliquity factor associated with Kirchhoff migration (Section H.1).

### G.5 Trace interpolation

In Section F.6 we developed the theory of the spatial prediction filter. The objective at the time was to attenuate random noise by applying the spatial prediction filter to data in the frequency-space domain [22].

Now, we want to design a complex Wiener prediction filter to interpolate data in the spatial direction. Consider a CMP-stacked data set P(x, t), where x is the CMP axis and t is the two-way zero-offset time axis. Apply Fourier transform in the t direction to decompose this 2-D data set to its frequency components P(x, ω). For each frequency component, define a complex array P : P(x, ω) in the x direction. Specifically, we want a filter F : F(x) such that, when applied to the input data array P : P(x, ω), it yields an estimate of the input array D : P(x + Δx/2, ω), at x + Δx/2, where D is the desired output array and Δx/2 is the prediction lag [24].

The spatial prediction filtering is expressed by the following convolutional relation

 ${\displaystyle P(x,\omega )\ast F(x)=Y(x),}$ (G-32)

where Y(x) represents the actual output from prediction filtering. Consider the discrete form of equation (G-32), with F(x) represented by the m-length complex array F : (F0, F1, F2, …, Fm−1), P(x, ω) represented by the n-length complex array P : (P0, P1, P2, …, Pn−1), and Y(x) represented by the (m + n − 1)-length complex array Y : (Y0, Y1, Y2, …, Ym+n−1). Equation (G-32) can then be expressed in matrix form

 ${\displaystyle {\begin{pmatrix}P_{0}&0&\ldots &0\\P_{1}&P_{0}&\ldots &0\\P_{2}&P_{1}&P_{0}&\vdots \\\vdots &P_{2}&P_{1}&P_{0}\\P_{m-1}&\vdots &P_{2}&P_{1}\\0&P_{m-1}&\vdots &P_{2}\\\vdots &\ldots &P_{m-1}&\vdots \\0&0&\ldots &P_{m-1}\end{pmatrix}}{\begin{pmatrix}F_{0}\\F_{1}\\F_{2}\\\vdots \\F_{n-1}\end{pmatrix}}={\begin{pmatrix}Y_{0}\\Y_{1}\\Y_{2}\\\vdots \\Y_{m+n-1}\end{pmatrix}}.}$ (G-33)

To develop the theory for trace interpolation using spatial prediction filters, consider the special case of a three-point one-step prediction filter (F0, F1, F2) and a six-point input data array (P0, P1, P2, P3, P4, P5). We also want to specialize the output array Y(x) to be the input array P(x) one unit of distance ahead, (P1, P2, P3, P4, P5). Equation (G-33) for this special case takes the form

 ${\displaystyle {\begin{pmatrix}P_{0}&0&0\\P_{1}&P_{0}&0\\P_{2}&P_{1}&P_{0}\\P_{3}&P_{2}&P_{1}\\P_{4}&P_{3}&P_{2}\\P_{5}&P_{4}&P_{3}\\0&P_{5}&P_{4}\end{pmatrix}}{\begin{pmatrix}F_{0}\\F_{1}\\F_{2}\end{pmatrix}}={\begin{pmatrix}P_{1}\\P_{2}\\P_{3}\\P_{4}\\P_{5}\\0\\0\end{pmatrix}}.}$ (G-34)

Augment the column array on the right-hand side to the matrix on the left-hand side to obtain

 ${\displaystyle {\begin{pmatrix}P_{1}&P_{0}&0&0\\P_{2}&P_{1}&P_{0}&0\\P_{3}&P_{2}&P_{1}&P_{0}\\P_{4}&P_{3}&P_{2}&P_{1}\\P_{5}&P_{4}&P_{3}&P_{2}\\0&P_{5}&P_{4}&P_{3}\\0&0&P_{5}&P_{4}\end{pmatrix}}{\begin{pmatrix}1\\-F_{0}\\-F_{1}\\-F_{2}\\\end{pmatrix}}={\begin{pmatrix}0\\0\\0\\0\\0\\0\\0\end{pmatrix}}.}$ (G-35)

Write out explicitly the equations associated with rows 3, 4, and 5 which contain terms with all three filter coefficients (F0, F1, F2) as

 ${\displaystyle P_{3}-P_{2}F_{0}-P_{1}F_{1}-P_{0}F_{2}=0,}$ (G-36a)

 ${\displaystyle P_{4}-P_{3}F_{0}-P_{2}F_{1}-P_{1}F_{2}=0,}$ (G-36b)

 ${\displaystyle P_{5}-P_{4}F_{0}-P_{3}F_{1}-P_{2}F_{2}=0.}$ (G-36c)

Next, split each of these three equations into two parts — first part has terms that contain the data array elements (P0, P2, P4) and second part has terms that contain the data array elements (P1, P3, P5). Then, write the two parts in matrix form

 ${\displaystyle {\begin{pmatrix}-F_{2}&-F_{0}&0\\0&-F_{1}&1\\0&-F_{2}&-F_{0}\end{pmatrix}}{\begin{pmatrix}P_{0}\\P_{2}\\P_{4}\end{pmatrix}}-{\begin{pmatrix}F_{1}&-1&0\\F_{2}&F_{0}&0\\0&F_{1}&-1\end{pmatrix}}{\begin{pmatrix}P_{1}\\P_{3}\\P_{5}\end{pmatrix}}=0.}$ (G-37)

This equation suggests that given the prediction filter (F0, F1, F2), the data array (P1, P3, P5) can be calculated from the data array (P0, P2, P4). The known data array (P0, P2, P4) is associated with three consecutive traces in a CMP-stacked data set. The unknown array to be calculated (P1, P3, P5) would then be associated with the traces halfway between the traces that correspond to the known data array (P0, P2, P4). Hence, we note that equation (G-37) can be used to perform spatial interpolation.

Write equation (G-37) in matrix notation

 ${\displaystyle {\mathbf {F} }_{d}{\mathbf {P} }_{d}-{\mathbf {F} }_{u}{\mathbf {P} }_{u}=0,}$ (G-38)

where Fd and Fu are the coefficient matrices with their elements in terms of the prediction filter coefficients (F0, F1, F2), and Pd : (P0, P2, P4) and Pu : (P1, P3, P5) are the known input data array and the interpolated data array, respectively.

Equation (G-38) is of the form d′ = Lu as for a generalized linear inverse problem (Section F.3). The least-squares solution u = (LT*L)−1LT*d is adopted for the interpolated data array Pu. In the present case, the matrix L corresponds to Fu, and the array d corresponds to FdPd, so that

 ${\displaystyle {\mathbf {P} }_{u}=({\mathbf {F} }_{u}^{T\ast }{\mathbf {F} }_{u})^{-1}{\mathbf {F} }_{u}^{T\ast }{\mathbf {F} }_{d}{\mathbf {P} }_{d},}$ (G-39)

where T denotes matrix transpose and the asterisk denotes complex conjugate.

The question remains as to how to calculate the prediction filter from the known data (P0, P2, P4). We shall make the assumption that the input data set is made up of M traces that comprise a set of N dipping events represented by the following combination of amplitude and phase spectra in the frequency-space domain [24]:

 ${\displaystyle P_{k}(\omega )=\sum _{j=1}^{N}A_{j}(\omega )\ \exp(-i\omega k\Delta \tau _{j}),k=0,1,2,\ldots ,M,}$ (G-40)

where Aj(ω) is the amplitude spectrum of the jth dipping event, Δτj is the time shift along the jth dipping event from trace to trace, and k is the trace index. Define

 ${\displaystyle z_{j}^{k}(\omega )=\exp(-i\omega k\Delta \tau _{j}),}$ (G-41)

so that equation (G-40) takes the form

 ${\displaystyle P_{k}=\sum _{j=1}^{N}A_{j}z_{j}^{k},k=0,1,2,\ldots ,M,}$ (G-42)

where the variable ω is omitted for convenience.

Now, given the data array (P0, P2, P4) and the one-step prediction filter (F0, F1, F2), we want to compute the next element of the data array P6. Refer to equation (G-36a) and replace the elements of the data array (P0, P1, P2) with the elements of the data array (P0, P2, P4). Also, replace the next element P3 of the data array (P0, P1, P2) with the next element P6 of the data array (P0, P2, P4) to get

 ${\displaystyle P_{4}F_{0}+P_{2}F_{1}+P_{0}F_{2}=P_{6}.}$ (G-43)

Write this equation using the form given by equation (G-42)

 ${\displaystyle \left(\sum _{j=1}^{N}A_{j}z_{j}^{4}\right)F_{0}+\left(\sum _{j=1}^{N}A_{j}z_{j}^{2}\right)+F_{1}\left(\sum _{j=1}^{N}A_{j}z_{j}^{0}\right)F_{2}=\sum _{j=1}^{N}A_{j}z_{j}^{6}.}$ (G-44)

Consider a case of three dipping events. Set N = 3 in equation (G-44) to get

 ${\displaystyle z_{1}^{4}F_{0}+z_{1}^{2}F_{1}+F_{2}=z_{1}^{6},}$ (G-45a)

 ${\displaystyle z_{2}^{4}F_{0}+z_{2}^{2}F_{1}+F_{2}=z_{2}^{6},}$ (G-45b)

 ${\displaystyle z_{3}^{4}F_{0}+z_{3}^{2}F_{1}+F_{2}=z_{3}^{6},}$ (G-45c)

which can be written in matrix form, while explicitly incorporating the variable ω [62],

 ${\displaystyle {\begin{pmatrix}z_{1}^{4}(\omega )&z_{1}^{2}(\omega )&1\\z_{2}^{4}(\omega )&z_{2}^{2}(\omega )&1\\z_{3}^{4}(\omega )&z_{3}^{2}(\omega )&1\end{pmatrix}}{\begin{pmatrix}F_{0}(\omega )\\F_{1}(\omega )\\F_{2}(\omega )\\\end{pmatrix}}={\begin{pmatrix}z_{1}^{6}(\omega )\\z_{2}^{6}(\omega )\\z_{3}^{6}(\omega )\\\end{pmatrix}}.}$ (G-46)

Equation (G-46) can be solved for the prediction filter (F0, F1, F2) for each frequency component of the input data array (P0, P2, P4). But this filter, when applied to the data array, yields the next element P6 (equation G-43). We need a filter that yields (P1, P3, P5, …), — the elements halfway between the elements of the data array (P0, P2, P4, …) equivalent to the output from interpolation.

To expedite a determination of this interpolation filter, derive the prediction filter ${\displaystyle (F_{0}^{\prime },\ F_{1}^{\prime },\ F_{2}^{\prime })}$ associated with the data array that would have been created by interpolation, (P0, P1, P2). Write equation (G-43) for the case of the data array (P0, P1, P2) to compute the next element P3 using the one-step prediction filter ${\displaystyle (F_{0}^{\prime },\ F_{1}^{\prime },\ F_{2}^{\prime })}$ as

 ${\displaystyle P_{2}F_{0}^{\prime }+P_{1}F_{1}^{\prime }+P_{0}F_{2}^{\prime }=P_{3}.}$ (G-47)

Write this equation using the form given by equation (G-42)

 ${\displaystyle \left({\sum \limits _{j=1}^{N}{{A_{j}}z_{j}^{2}}}\right){F_{0}^{\prime }}+\left({\sum \limits _{j=1}^{N}{{A_{j}}z_{j}^{1}}}\right){F_{1}^{\prime }}+\left({\sum \limits _{j=1}^{N}{{A_{j}}z_{j}^{0}}}\right){F_{2}^{\prime }}=\sum \limits _{j=1}^{N}{{A_{j}}z_{j}^{3}}.}$ (G-48)

For the case of three dipping events, setting N = 3 in equation (G-48), we get

 ${\displaystyle z_{1}^{2}F_{0}^{\prime }+z_{1}^{1}F_{1}^{\prime }+F_{2}^{\prime }=z_{1}^{3},}$ (G-49a)

 ${\displaystyle z_{2}^{2}F_{0}^{\prime }+z_{2}^{1}F_{1}^{\prime }+F_{2}^{\prime }=z_{2}^{3},}$ (G-49b)

 ${\displaystyle z_{3}^{2}F_{0}^{\prime }+z_{3}^{1}F_{1}^{\prime }+F_{2}^{\prime }=z_{3}^{3},}$ (G-49c)

which can be written in matrix form, while explicitly incorporating the variable ω,

 ${\displaystyle {\begin{pmatrix}z_{1}^{2}(\omega )&z_{1}^{1}(\omega )&1\\z_{2}^{2}(\omega )&z_{2}^{1}(\omega )&1\\z_{3}^{2}(\omega )&z_{3}^{1}(\omega )&1\\\end{pmatrix}}{\begin{pmatrix}F_{0}^{\prime }(\omega )\\F_{1}^{\prime }(\omega )\\F_{2}^{\prime }(\omega )\\\end{pmatrix}}={\begin{pmatrix}z_{1}^{3}(\omega )\\z_{2}^{3}(\omega )\\z_{3}^{3}(\omega )\\\end{pmatrix}}.}$ (G-50)

Now, note from the definition expressed by equation (G-41) that

 ${\displaystyle z_{j}^{2k}(\omega /2)=z_{j}^{k}(\omega ).}$ (G-51)

With this definition, rewrite the matrix equation (G-50) as

 ${\displaystyle {\begin{pmatrix}z_{1}^{4}(\omega /2)&z_{1}^{2}(\omega /2)&1\\z_{2}^{4}(\omega /2)&z_{2}^{2}(\omega /2)&1\\z_{3}^{4}(\omega /2)&z_{3}^{2}(\omega /2)&1\\\end{pmatrix}}{\begin{pmatrix}F_{0}^{\prime }(\omega )\\F_{1}^{\prime }(\omega )\\F_{2}^{\prime }(\omega )\end{pmatrix}}={\begin{pmatrix}z_{1}^{6}(\omega /2)\\z_{2}^{6}(\omega /2)\\z_{3}^{6}(\omega /2)\end{pmatrix}}.}$ (G-52)

Additionally, rewrite equation (G-46) explicitly in terms of ω/2 as

 ${\displaystyle {\begin{pmatrix}z_{1}^{4}(\omega /2)&z_{1}^{2}(\omega /2)&1\\z_{2}^{4}(\omega /2)&z_{2}^{2}(\omega /2)&1\\z_{3}^{4}(\omega /2)&z_{3}^{2}(\omega /2)&1\\\end{pmatrix}}{\begin{pmatrix}F_{0}(\omega /2)\\F_{1}(\omega /2)\\F_{2}(\omega /2)\end{pmatrix}}={\begin{pmatrix}z_{1}^{6}(\omega /2)\\z_{2}^{6}(\omega /2)\\z_{3}^{6}(\omega /2)\end{pmatrix}}.}$ (G-53)

Finally, compare equations (G-52) and (G-53) and note that the prediction filter ${\displaystyle [F_{0}^{\prime }(\omega ),F_{1}^{\prime }(\omega ),F_{2}^{\prime }(\omega )]}$ that we need to do trace interpolation of the data component with frequency ω actually is the prediction filter [F0(ω/2), F1(ω/2), F2(ω/2)] that is computed directly from the uninterpolated data component with half the frequency ω/2.

Figure G-2  Rayapth geometry for 3-D nonzero-offset traveltime equation (G-59) derived in Section G.6.

### G.6 3-D nonzero-offset traveltime equation

Refer to the 3-D raypath geometry shown Figure G-2. For convenience, we shall consider a 3-D recording parallel to the inline direction x with zero source-receiver azimuth. Specifically, consider a midpoint M : (x, y, 0) on crossline y that is associated with a source S : (x − h, y, 0) and a receiver R : (x + h, y, 0), where h is the half-offset. The distance along the raypath SD from source S : (x − h, y, 0) to a point D(0, 0, z) in the subsurface is given by

 ${\displaystyle SD={\sqrt {(x-h)^{2}+y^{2}+z^{2}}},}$ (G-54a)

and the distance along the raypath DR from point D : (0, 0, z) to a receiver R : (x + h, y, 0) is given by

 ${\displaystyle DR={\sqrt {(x+h)^{2}+y^{2}+z^{2}}}.}$ (G-54b)

The traveltime t associated with the raypath SDR then is given by

 ${\displaystyle vt={\sqrt {(x+h)^{2}+y^{2}+z^{2}}}+{\sqrt {(x-h)^{2}+y^{2}+z^{2}}},}$ (G-55)

where v is the medium velocity.

We shall show that equation (G-55) represents an ellipsoid (Figure G-2b). Square both sides to get

 ${\displaystyle {\begin{array}{l}v^{2}t^{2}=2{\sqrt {[(x+h)^{2}+y^{2}+z^{2}][(x-h)^{2}+y^{2}+z^{2}]}}\\\quad \quad \quad +[(x+h)^{2}+y^{2}+z^{2}]+[(x-h)^{2}+y^{2}+z^{2}].\end{array}}}$ (G-56)

Combine the second and third terms on the right-hand side and simplify the terms inside the square root

 ${\displaystyle v^{2}t^{2}=2{\sqrt {(x^{2}-h^{2})^{2}+2(y^{2}+z^{2})(x^{2}+h^{2})+z^{4}}}+2(x^{2}+h^{2}+y^{2}+z^{2}).}$ (G-57)

Perform further algebraic manipulation to collect the terms in x, y, and z

 ${\displaystyle (v^{2}t^{2}-4h^{2})x^{2}+v^{2}t^{2}y^{2}+v^{2}t^{2}z^{2}={\frac {v^{4}t^{4}}{4}}-v^{2}t^{2}h^{2}.}$ (G-58)

Finally, normalize by the terms on the right-hand side and rearrange the terms in the denominators to obtain

 ${\displaystyle {\frac {x^{2}}{(vt/2)^{2}}}+{\frac {y^{2}}{(vt/2)^{2}-h^{2}}}+{\frac {z^{2}}{(vt/2)^{2}-h^{2}}}=1.}$ (G-59)

Equation (G-59) represents an ellipsoid. The horizontal cross-section of this ellipsoid at a depth z is an ellipse. At z = 0, the ellipse has the following parameters:

1. Semi-major axis in inline x direction: a = vt/2.
2. Semi-minor axis in crossline y direction: ${\displaystyle b={\sqrt {(vt/2)^{2}-h^{2}}}}$.
3. Distance from center to either focus: ${\displaystyle {\sqrt {a^{2}-b^{2}}}=h}$.

The ellipsoid of equation (G-59) in the x − y − z volume describes the kinematics of the impulse response of a 3-D nonzero-offset migration operator applied to 3-D prestack data.

When equation (G-59) is specialized to the zero-offset case, h = 0, we get

 ${\displaystyle {\frac {x^{2}}{(vt/2)^{2}}}+{\frac {y^{2}}{(vt/2)^{2}}}+{\frac {z^{2}}{(vt/2)^{2}}}=1,}$ (G-60a)

which describes a hollow hemisphere in the (x, y, z) volume for a constant t with a radius vt/2. This hemisphere represents the kinematics of the impulse response of a 3-D zero-offset migration operator applied to 3-D poststack data.

When equation (G-55) is specialized to the 3-D zero-offset case, h = 0, we get

 ${\displaystyle vt=2{\sqrt {x^{2}+y^{2}+z^{2}}},}$ (G-60b)

which describes the diffraction hyperboloid of revolution in the (x, y, t) volume for a constant z.[63] [64] [65]

## References

1. Levin, 1983, Levin, F. K., 1983, The effects of streamer feathering on stacking: Geophysics, 48, 1165–1171.
2. Levin 1984, Levin, F. K., 1984, The effect of binning on data from a feathered streamer: Geophysics, 49, 1386–1387.
3. Bentley and Yang, 1982, Bentley, L. and Yang, M., 1982, Scatter of midpoints grouped in cells and its effects on stacking: unpublished technical document, Western Geophysical Company.
4. Larner and Ng, 1984, Larner, K. and Ng, P., 1984, Strike versus dip shooting: 53th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 256–260.
5. Levin (1971), Levin, F. K., 1971, Apparent velocities from dipping interface reflections: Geophysics, 36, 510–516.
6. Lehmann and Houba (1985), Lehmann, H. J. and Houba, W., 1985, Practical aspects of determination of 3-D stacking velocities: Geophys. Prosp., 33, 155–163.
7. Deregowski and Rocca, 1981, Deregowski, S. and Rocca, F., 1981, Geometrical optics and wave theory for constant-offset sections in layered media: Geophys. Prosp., 29, 374–387.
8. Forel and Gardner, 1988, Forel, D. and Gardner, G. H. F., 1988, A three-dimensional perspective on two-dimensional dip-moveout: Geophysics, 53, 604–610.
9. Meinardus and Schleicher, 1993, Meinardus, H. A. and Schleicher, K. L., 1993, 3-D time-variant dip moveout by the f − k method: Geophysics, 58, 1030–1041.
10. Ferber, 1991, Ferber, R., 1991, A filter, delay and spread technique for 3-D DMO: Geophys. Prosp., 39, 737–755.
11. Beasley and Klotz, 1992, Beasley, C. J. and Klotz, R., 1992, Equalization of DMO for irregular spatial sampling: 58th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 970–973.
12. Goodway and Ragan (1995), Goodway, W. N. and Ragan, B., 1995, Focused 3-D: consequences of midpoint scatter and spatial sampling in acquisition design, processing and interpretation: Presented at the Can. Soc. Expl. Geophys. Mtg.
13. Brink and Ronen (1998), Brink, M. and Ronen, S., 1998, Coverage and binning issues for marine seismic surveys: The Leading Edge, 17, 1600–1605.
14. Beasley and Mobley, 1998, Beasley, C. J. and Mobley, E., 1998, Spatial dealiasing of 3-D DMO: The Leading Edge, 17, 1590–1594.
15. Ronen et al., 1991, Ronen, S., Sorin, V., and Bale, R., 1991, Spatial dealiasing of 3-D seismic reflection data: Geophys. J. Internat., 105, 503–511.
16. Ronen, 1985, Ronen, J., 1985, Multichannel inversion in reflection seismology: Ph. D. thesis, Stanford University.
17. Ronen, 1994, Ronen, S., 1994, Handling irregular geometry: Equalized DMO and beyond: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1545–1548.
18. Ronen et al., 1995, Ronen, S., Nichols, D., Bale, R. and Ferber, R., 1995, Several field data examples of dealiasing DMO: 66th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1168–1171.
19. Bale et al., 1996, Bale, R., Ferber, R., Nichols, D., Stolte, C. and Ronen, S., 1996, Several field data examples of dealiasing DMO: 66th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1168–1171.
20. Rothman et al., 1981, Rothman, D., Larner, K. L., and Chambers, R., 1981, Trace interpolation and design of 3-D surveys: Presented at the 39th Ann. Mtg. Eur. Assoc. Expl. Geophys.
21. Ronen, 1987, Ronen, J., 1987, Wave-equation trace interpolation: Geophysics, 52, 973–984.
22. Canales, 1984, Canales, L., 1984, Random noise reduction: 54th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 525–527.
23. Gulunay, 1986, Gulunay, N., 1986, F − X decon and complex Wiener prediction filter: 56th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 279–281.
24. Spitz, 1991, Spitz, S., 1991, Seismic trace interpolation in the f − x domain: Geophysics, 56, 785–794.
25. French (1974), French, W. S., 1974, Two-dimensional and three-dimensional migration of model-experiment reflection profiles: Geophysics, 39, 265–277.
26. Schneider (1978), Schneider, W. A., 1978, Integral formulation for migration in two and three dimensions: Geophysics, 43, 49–76.
27. Jacubowicz and Levin (1983), Jakubowicz, H. and Levin, S., 1983, A simple exact method of 3-D migration — theory: Geophys. Prosp., 31, 34–56.
28. Claerbout, 1985, Claerbout, J. F., 1985, Imaging the earth’s interior: Blackwell Scientific Publications.
29. Brown, 1983, Brown, D. L., 1983, Applications of operator separation in reflection seismology: Geophysics, 48, 288–294.
30. Li (1991), Li., Z, 1991, Compensating finite-difference errors in 3-D migration and modeling: Geophysics, 56, 1650–1660.
31. Etgen and Nichols (1999), Etgen, J. T. and Nichols, D., 1999, Application of the Li correction to explicit depth extrapolation methods: 69th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1366–1369.
32. Holberg, 1988, Holberg, O., 1988, Towards optimum one-way wave propagation: Geophys. Prosp., 36, 99–114.
33. Hale, 1991a, Hale, D., 1991a, Stable explicit depth extrapolation of seismic wavefields: Geophysics, 56, 1770-1777.
34. Hale, 1991b, Hale, D., 1991b, 3-D depth migration via McClellan transforms: Geophysics, 56, 1778–1785.
35. Blacquiere, 1991, Blacquiere, G., 1991, Optimized McClellan transformation filters applied to one-pass 3-D depth extrapolation filters: 61st Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1126-1129.
36. Biondi and Palacharla, 1994, Biondi, B. and Palacharla, G., 1994, 3-D migration using rotated McClellan filters: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1278–1281.
37. Gaiser, 1994, Gaiser, J. E., 1994, Optimum transformation and filter structure for explicit finite-difference 3-D migration: Presented at the 56th Mtg. Eur. Assoc. Expl. Geophys.
38. Notfors, 1995, Notfors, C. D. B., 1995, Accurate and efficient explicit 3-D migration: 65th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1224–1227.
39. Kosloff and Kessler, 1987, Kosloff, D. and Kessler, D., 1987, Accurate depth migration by a generalized phase-shift method: Geophysics, 52, 1074–1084.
40. Reshef and Kessler, 1989, Reshef, M. and Kessler, D., 1989, Practical implementation of three-dimensional poststack depth migration: Geophysics, 54, 309–318.
41. Ristow (1980), Ristow, D., 1980, 3-D downward extrapolation of seismic data in particular by finite-difference methods: PhD thesis, University of Utrecht, The Netherlands.
42. Black et al., 1987, Black, J. L. and Leong, T. K., 1987, A flexible, accurate approach to one-pass 3-D migration: 57th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 559–560.
43. Kitchenside and Jacubowicz, 1987, Kitchenside, P., and Jacubowicz, H., 1987, Operator design for 3-D depth migration: 57th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 556–558.
44. Kitchenside, 1988, Kitchenside, P., 1988, Steep-dip 3-D migration: Some issues and examples: 58th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 976–978.
45. Blacquiere et al. (1989), Blacquiere, G., Debeye, H. W. J., Wapenaar, C. P. A., and Berkhout, A. J., 1989, 3-D table-driven migration; Geophys., Prosp., 37, 925–958.
46. Dickinson (1988), Dickinson, J. A., 1988, Evaluation of two-pass three-dimensional migration: Geophysics, 53, 32–49.
47. Soubaras (1996), Soubaras, R., 1996, Explicit 3-D migration using equiripple polynomial expansion and Laplacian synthesis: Geophysics, 61, 1386–1393.
48. Pai, 1988, Pai, D. M., 1988, Generalized f − k (frequency-wavenumber) migration in arbitrarily varying media: Geophysics, 53, 1547–1555.
49. Berryhill, 1991, Berryhill, J. R., 1991, Kinematics of crossline prestack migration: Geophysics, 56, 1674–1676.
50. Canning and Gardner, 1996, Canning, A. and Gardner, G. H. F., 1996, A two-pass approximation to 3-D prestack depth migration: Geophysics, 61, 409–421.
51. Marcoux et al., 1987, Marcoux, M. O., Godfrey, R. J., and Notfors, C. D., 1987, Migration for optimum velocity evaluation and stacking: Presented at the 49th Ann. Mtg. European Ass. Expl. Geophys.
52. Ferber, 1994, Ferber, R., 1994, Migration to multiple offset and velocity analysis: Geophys. Prosp., 42, 99–112.
53. Devaux et al., 1996, Devaux, V., Gardner, G. H. F., and Rampersad, T., 1996, 3-D prestack depth migration by Kirchhoff operator splitting: 66th Ann. Internal. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 455–458.
54. Baysal and Kosloff, 1998, Baysal, E. and Kosloff, D., 1998, Constructing 2-D prestack data from 3-D volume via two-pass 3-D Kirchhoff migration: Open-file technical document, Paradigm Geophysical.
55. Fowler (1984), Fowler, P., 1984, Velocity-independent imaging of seismic reflectors: 54th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 383–385.
56. Bone et al. (1975), Bone, M. R., Giles, B. F., and Tegland, E. R., 1975, Analysis of seismic data using horizontal cross-sections: Presented at the 45th Ann. Internat. Mtg., Soc. Expl. Geophys.
57. Brown (1978), Brown, A. R., 1978, 3-D seismic interpretation methods: Presented at the 48th Ann. Internat. Soc. Expl. Geophys. Mtg.
58. Brown et al. (1982), Brown, A. R., Graebner, R. J., and Dahm, C. G., 1982, Use of horizontal seismic sections to identify subtle traps, in The Deliberate Search for the Subtle Trap: Am. Assoc. Petr. Geol. Memoir, 32, 47–56.
59. Schultz and Lau (1984), Schultz, P. S. and Lau, A., 1984, Poststack estimation of three-dimensional crossline statics: Geophysics, 49, 227–236.
60. Gazdag, 1978, Gazdag, J., 1978, Wave-equation migration by phase shift: Geophysics, 43, 1342–1351.
61. Stolt (1978), Stolt, R. H., 1978, Migration by Fourier transform: Geophysics, 43, 23–48.
62. Porsani, 1997, Porsani, M. J., 1997, Seismic trace interpolation using half-step prediction filters: Preprint.
63. Berkhout, A. J., 1985, 3-D seismic processing with an eye to the future: World Oil, 201, No. 5, 91–94.
64. Biondi, B., Fomel, S., and Chemingui, N., 1998, Azimuth moveout for 3-D prestack imaging: Geophysics, 63, 574–588.
65. Chemingui, N. and Biondi, B., 1997, Equalization of irregular data by iterative inversion: 67th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1115–1118.