Encyclopedia Of Geology

GEOPHYSICS, VOL. 66, NO. 6 (NOVEMBER-DECEMBER 2001); P. 1699–1713, 17 FIGS.

Case History A unified 3-D seismic workflow

¨ Yilmaz∗ , Irfan Tanir‡ , and Cyril Gregory‡ Oz

ABSTRACT

recorded data and finite cable length used in recording further compound the problem. Additionally, traveltime picking that is needed for most velocity estimation techniques and time-to-depth conversion as well as picking depth horizons from depth-migrated data to delineate reflector geometries are all adversely affected by noise present in the data. All things considered, we can only expect to do our best in estimating what may be called an initial model, and update this model to get an acceptable final model. The objective behind the design of the seismic workflow described in this paper is to attain the best estimate of a structurally consistent initial model based on rms velocities associated with migrated data, so as to minimize the work required to update the model. The unified workflow involves analysis of seismic data both in time and depth, and follows a pathway that starts with the application of 3-D dip-moveout correction and 3-D prestack time migration to derive an rms velocity field. This is followed by estimation of an accurate, structurally consistent initial model by Dix conversion of rms velocities and interpretation of a set of depth horizons from 3-D poststack depth migration. To update the initial model, the image gathers derived from 3-D prestack depth migration are analyzed for residual moveout. The resulting final model is then used to perform 3-D prestack depth migration to obtain an image volume in depth. The final phase of the workflow includes structural and stratigraphic interpretation of the image volume with the ultimate objective of obtaining a seismically derived reservoir model.

A geophysicist who practices seismic data analysis for earth modeling and imaging in depth is overwhelmed by the prolific number of inversion methods to estimate layer velocities and delineate reflector geometries—the two constituents of a seismically defined earth model in depth. Given a specific type of structural play, the key to estimating an accurate earth model in depth, however, is a workflow that is based on a judicious combination of inversion methods appropriately selected for their robustness. We present a unified workflow for processing, inversion, and interpretation of 3-D seismic data that is applicable to low-relief structures and complex structures associated with extensional and compressional tectonics. With some modifications, the workflow also is applicable to complex overburden structures associated with salt and overthrust tectonics. Although doing it right the first time is most desirable, there is never a situation where this is possible when estimating an earth model in depth. A fundamental problem with inversion for earth modeling is velocity-depth ambiguity. This means that an error in layer velocity can be indistinguishable from an error in reflector geometry. The velocity-depth ambiguity that is inherent to seismic inversion makes it very difficult to obtain the right answer (an adequate representation of the true geological model), let alone do it the first time. Limitations in the resolving power of the methods to estimate layer velocities that arise from the band-limited nature of the

Manuscript received by the Editor July 26, 2000; revised manuscript received February 20, 2001. ∗ Anatolian Geophysical, 28 Lalebayiri, Kemer Country, Kemerburgaz, Istanbul 34993, Turkey. E-mail: [email protected]. ‡Paradigm Geophysical, Chobham House, Christ Church Way, Woking, Surrey GU21 1JF, United Kingdom. E-mail: [email protected]; [email protected].  c 2001 Society of Exploration Geophysicists. All rights reserved. 1699

1700

Yilmaz et al. INTRODUCTION

We describe the unified seismic workflow by a 3-D structural inversion case study from the onshore northeast China. The case study involves seismic data from a 3-D survey that covers an area of nearly 150 km2 . The surface area is roughly a rectangle with dimensions of approximately 13 km in the in-line direction and 12 km in the cross-line direction. Survey statistics include more than 1.8 million recorded traces and nearly 120,000 stacked traces. The average fold is 15, and the in-line and cross-line trace spacings are 25 and 50 m, respectively. There are 19 swaths of data recorded using four cables that are 100 m apart, each with 60 receiver groups at 50-m interval. The shots were placed in the direction perpendicular to the swaths. Prior to the execution of the unified seismic workflow, a study of the fold of coverage was made to examine its variation with respect to source-receiver offset and azimuth ranges. Nonuniform fold, irregular offset, azimuth, and midpoint distributions within common-cell gathers are issues that need to be addressed in 3-D data analysis. Lack of uniformity in the fold 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 migration, and amplitude variation with offset (AVO) analysis are adversely affected by the acquisition footprint characterized by nonuniform fold, irregular spatial sampling, and undersampling of the recorded data. Although these issues can be addressed to some extent during data analysis (Beasley and Klotz, 1992; Ronen, 1994; Ronen et al., 1995; Beasley and Mobley, 1998), they must be avoided during data recording as much as possible. During the acquisition of the data used in this case study, exceptional effort was made to attain a fairly uniform fold over 90% of the survey area. Uniformity in fold also extends to source-receiver azimuthal coverage and offset distribution in common-cell gathers. We execute the unified seismic workflow in five phases, as outlined in Figure 1. A brief description of each phase is as follows: 1) 3-D DMO processing with deliverables that include a set of 3-D DMO gathers, a 3-D DMO stack volume, a 3-D DMO velocity field, and an image volume derived from 3-D poststack time migration. 2) 3-D prestack time migration with deliverables that include a set of common-reflection-point (CRP) gathers, an image volume from 3-D prestack time migration, a 3-D rms velocity field, and a 3-D zero-offset wavefield modeled from the image volume. The latter may be treated as the equivalent of an unmigrated stack volume associated with the 3-D prestack time-migrated data. 3) Stratigraphic inversion with deliverables that include the AVO attribute volumes derived from prestack amplitude inversion and an acoustic impedance volume derived from poststack amplitude inversion. In this phase, we also derive a 3-D interval velocity field by Dix conversion of the 3-D rms velocity field from phase 2 and use it to create an image volume from 3-D poststack depth migration. 4) Structural inversion that involves estimation of an initial 3-D velocity-depth model and model updating with deliverables that include a final 3-D velocity-depth model

and an image volume derived from 3-D prestack depth migration. 5) Structural and stratigraphic interpretation, and calibration to well data with the ultimate objective of building a reservoir model that is fundamentally a petrophysical representation of an earth model in depth. In the present case study, we execute the entire unified workflow outlined above, except for stratigraphic inversion in phase 3 and reservoir modeling in phase 5. 3-D DMO PROCESSING

Shown in Figure 2a is a shot record with 240 traces that correspond to the four 60-channel receiver cables. The quality of the recorded data in general is remarkably good, except for some coherent, dispersive low-frequency and high-amplitude noise at near offsets that may be associated with the ground-roll energy. Prestack signal processing of the data included geometric spreading correction using a t 2 -scaling function, spiking deconvolution, and time-variant spectral whitening (Figure 2b). This parsimonious sequence sufficed to attain a desired broad and flat spectrum within the passband (compare Figures 2c and 2d). Note that the coherent linear noise that dominates the raw shot record (Figure 2a) does not stand out in the shot gather after spectral balancing (Figure 2b). Following the signal processing, the data are ready for DMO processing. The steps of the process are as follows: 1) Sort the common-shot gathers to common-cell gathers and perform stacking velocity analysis over a sparse grid of 2 × 2 km. 2) Apply NMO correction using flat-event velocities followed by 3-D DMO correction to remove source-receiver azimuth and dip effects on the stacking velocities. 3) Apply inverse NMO correction and repeat the velocity analysis over a finer grid of 0.5 × 0.5 km.

FIG. 1. A unified seismic workflow which begins with field data after DMO processing on top left and ends with an estimate of a reservoir model at bottom right. The arrows indicate the data flow, and the numbers on the left refer to different phases in the workflow: (1) 3-D DMO processing, (2) 3-D prestack time migration, (3) stratigraphic inversion, (4) structural inversion, and (5) structural and stratigraphic interpretation.

A Unified 3-D Seismic Workflow

4) Create a 3-D DMO velocity field from the vertical functions picked in step 3. 5) Apply NMO correction to the 3-D DMO-corrected data from step 3 using the velocity field from step 4 and stack the data. 6) Perform 3-D poststack time migration using a laterally smoothed form of the 3-D DMO velocity field from step 4. To circumvent the adverse effect of spatial aliasing on migration, especially in the cross-line direction, trace interpolation before migration was considered. However,

1701

this process deteriorated the stacked data in zones of intensive faulting associated with extensional tectonism. Therefore, in lieu of trace interpolation, the unmigrated data in step 5 were filtered down to a passband of 6–36 Hz. The deliverables from phase 1 (3-D DMO processing) include a set of 3-D DMO-corrected gathers, a volume of 3-D DMO-stacked data, a volume of 3-D DMO velocity field, and an image volume derived from 3-D poststack time migration. 3-D PRESTACK TIME MIGRATION

While the interpreter may have a sneak preview of the subsurface image using the volume of 3-D poststack time-migrated data from phase 1, we move on to phase 2 for 3-D prestack time migration. The improvement in imaging with 3-D prestack time migration may sometimes be marginal compared to imaging with 3-D poststack time migration of a 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 follows: 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 CRP gathers from 3-D prestack time migration are used to perform prestack amplitude inversion to derive the AVO attributes. 3) The 3-D image volume derived from 3-D prestack time migration can be used as input to 3-D zero-offset amplitude inversion to estimate an acoustic impedance model of the earth. 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 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 zerooffset 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.

FIG. 2. (a) A common-shot gather recorded by using four cables each with 60 receiver groups, (b) the same gather after t 2 -scaling, spiking deconvolution, and time-variant spectral whitening, (c) time-variant amplitude spectrum of the record in (a), and (d) time-variant amplitude spectrum of the record in (b). In (c) and (d), the horizontal axis is frequency and the vertical axis is the two-way time as in (a) and (b), and the graphs below the time-variant spectra are the amplitude spectra averaged over the time axis.

We follow a prestack time-migration strategy based on DMO correction combined with common-offset migration (Marcoux et al., 1987). Aside from prestack imaging, this strategy provides the opportunity to derive a velocity field associated with events in their migrated positions. The process that underlies this strategy yields 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. The process is as follows: 1) Sort the 3-D DMO- and NMO-corrected data from phase 1 to common-offset volumes. 2) Assume that each of the common-offset volumes of data is a replica of a 3-D zero-offset wavefield and perform

1702

3)

4)

5)

6)

7)

Yilmaz et al.

3-D zero-offset time migration using a regionally averaged, but vertically varying velocity function derived from the 3-D DMO velocity field from phase 1. Again, to circumvent the adverse effect of spatial aliasing, the common-offset volumes of data were filtered down to a passband of 6–36 Hz prior to migration. If, instead, trace interpolation was desired to circumvent spatial aliasing, it must preserve relative amplitudes of the CRP gathers from prestack time migration and not distort AVO. An interpolation scheme that meets this requirement is based on frequency-space Wiener prediction filtering (Spitz, 1991). Sort the migrated common-offset volumes of data back to common-cell gathers, apply inverse NMO correction using the 3-D DMO velocity field from phase 1, and repeat the velocity analysis over a grid of 0.5 × 0.5 km. Create a 3-D rms velocity field associated with the migrated data from the vertical functions picked in step 3. Figure 3a shows the cross-section of the 3-D rms velocity volume along an in-line traverse. Apply NMO correction to the common-cell gathers from step 3 using the rms velocity field from step 4, and stack the data. Model a 3-D zero-offset wavefield from the image volume derived from the 3-D common-offset migration and stack in step 5 using the same vertically varying velocity function as in step 2. Figure 3b shows an inline section from the modeled 3-D zero-offset wavefield volume. The final step in phase 2 involves remigration of the stacked data. Specifically, perform 3-D poststack time migration for which the input data volume is the 3-D zero-

offset modeled wavefield from step 6 and the migration velocity field is the 3-D rms velocity field from step 4. Figure 4 shows an in-line and a cross-line section from the image volume derived from 3-D prestack time migration based on the above sequence. Whereas the benefits of 3-D prestack time migration are not limited to the improved image of the subsurface as listed above, note from these two sections that the basin boundaries and the internal structural characteristics associated with extensional and compressional tectonism have been delineated accurately. The deliverables from phase 2 (3-D prestack time migration) include a set of CRP gathers, a volume of the 3-D zero-offset wavefield which may be considered as equivalent to an unmigrated stack volume, a volume of the 3-D rms velocity field, and an image volume derived from 3-D prestack time migration. 3-D MIGRATION VELOCITY ANALYSIS

In the previous section, we demonstrated that a strategy based on 3-D DMO correction combined with 3-D commonoffset migration can include a step for estimating a 3-D rms velocity field from the migrated prestack data in phase 2 of the unified workflow shown in Figure 1. Now, we develop an alternative strategy for 3-D migration velocity analysis to derive the 3-D rms velocity field based on a 3-D extension of the velocity-independent imaging technique by Fowler (1984): 1) Start with 3-D prestack data P(x, y, h, t) in coordinates of in-line x, offset 2h, cross-line y, and event time t in

FIG. 3. Phase 2 deliverables: (a) cross-section of the 3-D rms velocity field along an in-line traverse, and (b) cross-section of the 3-D zero-offset wave-field along the same traverse.

A Unified 3-D Seismic Workflow

2)

3)

4)

5)

6)

7)

the unmigrated position, and apply NMO correction, 3-D DMO correction followed by inverse NMO correction. 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 5a is a CVS volume created by using a velocity of 2400 m/s. Migrate the 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 CVM algorithm. Shown in Figure 5b is the CVM volume associated with the CVS volume of Figure 5a with a velocity of 2400 m/s. Repeat steps 2 and 3 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. Extract the time slices from the CVM volumes for a specified time τ j and combine them to form a velocity slab P(x, y, vmig ; τ j ). Figure 6a shows such a velocity slab 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 slabs was created for times from 600 to 3100 ms with an increment of 100 ms. For 3-D visualization and interpretation, form a supervolume that is a composite of the 26 velocity slabs from step 5 placed one on top of the other (Figure 7). Just as you would pick the peak semblance value at an event time on a conventional velocity spectrum, interpret each of the velocity slabs by picking maximum

1703

event amplitudes along selected cross-line (or in-line) traverses and create velocity strands associated with the constant times τ j , j = 1, 2, . . . , N (Figure 6b). 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 8). 9) Create a 3-D rms velocity field using the rms velocity maps from step 8. Figure 9 shows an in-line and a crossline section from the rms velocity volume. Compare with the corresponding image sections from time migration shown in Figure 4 and note that the 3-D rms velocity field honors the structural characteristics of the subsurface. Specifically, the in-line section in Figure 9a manifests the basin boundary indicated by the change in color from green (the basin interior) to orange red. The cross-line section in Figure 9b manifests 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 strategy outlined above is associated with events in their migrated positions. Compared with 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. If you wish to extend your analysis from the time to the depth domain, the 3-D rms velocity field associated with the timemigrated data also is the preferred velocity field to derive a

FIG. 4. Two cross-sections from the image volume derived from 3-D prestack time migration in phase 2: (a) an in-line section and (b) a cross-line section.

1704

Yilmaz et al.

3-D interval velocity field that can then be used to obtain an image in depth. FROM RMS TO INTERVAL VELOCITIES

representation of the reflection amplitudes at vertical incidence and zero offset. Alternatively, poststack amplitude inversion can be applied to the image volume from 3-D prestack time migration derived in phase 2.

As we note from the list of deliverables from phase 2 (Figure 1), we are not implementing phase 2 (3-D prestack time migration) just for imaging the subsurface. We make use of the 3-D rms velocity field and the 3-D zero-offset wavefield from phase 2 to move from the time to the depth domain in the analysis. Although it will not be considered in the present case study, phase 3 of the generalized workflow also includes amplitude inversion that uses the CRP gathers from 3-D prestack time migration and the image volume itself (Figure 1). The steps in phase 3 are as follows: 1) Apply prestack amplitude inversion to the CRP gathers and derive the AVO attribute volumes. These include the intercept and gradient volumes, P-wave and S-wave reflectivity and impedance volumes, and pseudo-Poisson and fluid-factor volumes. To satisfy the underlying assumptions for AVO analysis, and specifically the Zoeppritz equations that describe the partitioning of energy associated with an incident compressional plane wave at a horizontal layer boundary into its reflected and refracted compressional- and shear-wave components, we want to use the CRP gathers from 3-D prestack time migration that contain events in their migrated positions in prestack amplitude inversion rather than the DMO gathers that contain events in their unmigrated positions. 2) Apply poststack amplitude inversion to the AVO intercept volume to derive the acoustic impedance volume. It is assumed that the AVO intercept attribute is a close

FIG. 6. (a) A velocity slab formed by combining the time slices from a set of CVM volumes as in Figure 5b at the same time, and (b) the velocity strands interpreted from the velocity slab in (a). The event amplitudes in (a) are colorcoded: black is positive, red is negative, and yellow is the vicinity of zero crossing.

FIG. 5. (a) A CVS volume, and (b) a CVM volume. See text for details.

A Unified 3-D Seismic Workflow

3) Apply Dix conversion to the 3-D rms velocity field (Figure 3a) from phase 2 to derive a 3-D interval velocity field. Figure 10a shows an in-line cross-section from the interval velocity volume. To satisfy the vertical incidence assumption for Dix conversion, we want to use the 3-D rms velocity field derived in conjunction with 3-D prestack time migration of the data in phase 2 that is associated with events in their migrated positions rather than the DMO velocity field in conjunction with 3-D DMO correction of the data in phase 1 that is associated with events in their unmigrated positions. 4) Perform 3-D poststack depth migration of the 3-D zerooffset wavefield from phase 2 using the 3-D interval velocity field from step 3. The deliverables from phase 3 include a volume of 3-D interval velocity field and an image volume derived from 3-D poststack depth migration. STRUCTURAL INVERSION

We want to construct a structurally consistent 3-D velocitydepth model using the deliverables from phase 3 and update it to obtain a final 3-D velocity-depth model. We then want to use this earth model in depth to create an earth image in depth from 3-D prestack depth migration. The steps in phase 4 (Figure 1) are as follows:

FIG. 7. A portion of the super volume composed of the velocity slabs as in Figure 6a.

1705

1) Interpret a set of depth horizons from the image volume derived from 3-D poststack depth migration in phase 3 (Figure 10b). These depth horizons correspond to layer boundaries with significant velocity contrast. 2) To preserve the vertical and lateral velocity variations in the 3-D interval velocity field derived from Dix conversion (Figure 10a), create a set of phantom depth horizons by subdividing each of the layers bounded by the depth horizons (Figure 10b) interpreted from the 3-D poststack depth-migrated volume of data into a set of 10 sublayers. The phantom horizons can be conformed either to the top boundary or the base boundary or, as in the present case, to both the top and base boundary of the layer under consideration. Figure 10c shows the principal depth horizons as in Figure 10b interleaved with the phantom horizons in green. 3) Intersect the 3-D interval velocity volume (Figure 10a) from phase 3 with the principal and phantom depth horizons (Figure 10c) from step 2 and extract horizonconsistent interval velocity surfaces. Then, combine these layer velocities with the reflector geometries represented by the depth horizons to create an initial, structurally consistent 3-D velocity-depth model. Figure 10d shows an in-line cross-section from the initial 3-D velocity-depth model. 4) Perform 3-D prestack depth migration and generate a set of image gathers along every tenth in-line traverse,

FIG. 8. A set of rms velocity maps created by gridfitting the velocity strands as in Figure 6b. Each map represents the cross-section of the 3-D rms velocity field at a constant time associated with the prestack time-migrated data. The horizontal axis is in the in-line direction and the vertical axis is in the cross-line direction. Each map is colorcoded independently: purple represents low velocities, red to dark brown represents high velocities, and green represents medium-range velocities.

1706

Yilmaz et al.

and compute the residual moveout semblance spectra for model updating. To explain the basis for residualmoveout analysis for model updating, consider conventional stacking velocity analysis of a common-midpoint (CMP) gather. Compute the velocity spectrum from the CMP gather and pick a velocity function. Following the NMO correction of the CMP gather using this velocity function, events should look flat if the velocity function had been picked correctly. If the picking was done incorrectly, then you would observe events with residual moveout. In principle, this residual moveout can be computed and used to update the initially picked velocity function. Analogous to the conventional stacking velocity analysis, if the initial velocity-depth model has been estimated with sufficient accuracy, then image gathers derived from prestack depth migration using this model should exhibit flat events associated with the layer boundaries included in the model. Any errors in layer velocities and/or reflector geometries, on the other hand, should give rise to residual moveout along those events on the image gathers. Assume that the residual moveout is parabolic and compute the semblance spectra. Figures 11a and 11b show selected panels of residual-moveout analysis. Each panel comprises the image gather at the analysis location and the residual moveout semblance spectrum computed from it. The horizontal axis of the semblance plane represents the depth error, and the vertical axis represents the depth of the event. The semblance spectrum has two quadrants that correspond to positive and negative residual moveouts or, equivalently, to positive and negative depth errors. A flat event would yield a semblance peak that coincides with the vertical axis that coincides with

zero depth error, whereas an event with residual moveout would yield a semblance peak situated either in the left or right quadrant depending on the sign of the depth error. While many of the events in the image gathers show very small or no residual moveout (Figure 11a), note that some events exhibit significant residual moveout errors (Figure 11b). Perform residual moveout corrections (Deregowski, 1990) and update the initial velocity-depth model. 5) Repeat the residual moveout analysis a few times until the residual moveout errors have been reduced significantly. The residual moveout analysis from the last iteration of model updating shown in Figures 11c and 11d indicates significant reduction of the residual moveout errors (compare, for instance, Figure 11b with Figure 11d). We remind ourselves that the spirit with which the seismic workflow described in this paper is designed, as stated in the introductory section, is to attain the best estimate of a structurally consistent initial model based on rms velocities associated with migrated data, so as to minimize the work required to update the model. An otherwise erroneous initial estimate often causes excessive number of iterations which fails to minimize the residual moveout errors. Figures 12a and 12b show in-line and cross-line cross-sections, respectively, from the final velocity-depth model. 6) Perform 3-D prestack depth migration using the final 3-D velocity-depth model from step 5. The input to 3-D prestack depth migration is the prestack data set with signal processing applied as in Figure 2b. Figures 12c and 12d show selected in-line and cross-line sections, respectively, from the image volume in depth derived from 3-D

FIG. 9. Two cross-sections from the 3-D rms velocity volume derived from the rms velocity maps as in Figure 8: (a) an in-line section and (b) a cross-line section. Each section is colorcoded independently: purple represents low velocities, red to dark brown represents high velocities, and green represents medium-range velocities. Compare (a) with the cross-section shown in Figure 3a.

A Unified 3-D Seismic Workflow

prestack depth migration. Compare Figures 12a and 12b with the image sections in Figures 12c and 12d, and note the fine-grained, yet structurally consistent, lateral and vertical variations in velocities. Moreover, compare the image in depth in Figure 12c with the image in time in Figure 4, and note that the imaging of the fault planes within the basin and the basin boundary itself has been greatly improved by 3-D prestack depth migration. The deliverables from phase 4 (structural inversion) include a volume of structurally consistent 3-D interval velocity field and an image volume derived from 3-D prestack depth migration.

1707

3-D PRESTACK DEPTH MIGRATION

In the section on 3-D prestack time migration, we defined a strategy based on the application of 3-D zero-offset migration theory to 3-D common-offset data. Specifically, following NMO and 3-D DMO corrections, prestack data will have been corrected for dip and source-receiver azimuth effects. As a result, the common-offset volumes of data are decoupled representations of a 3-D zero-offset wavefield; thus, each common-offset volume can be migrated independently using a 3-D zero-offset time migration algorithm. Decoupling of common-offset data volumes, in a strict theoretical sense, is only possible for events with hyperbolic

FIG. 10. (a) A cross-section along an in-line traverse as in Figure 3a from the interval velocity volume derived from Dix conversion of the rms velocity volume, (b) depth horizons interpreted from the image volume derived from 3-D poststack depth migration of the 3-D zero-offset wavefield as in Figure 3b using the 3-D interval velocity field as in (a), (c) phantom horizons (in green) inserted between the principal horizons as in (b), and (d) the same interval velocity field as in (a) but transcribed into a structurally consistent form using the horizons in (c).

1708

Yilmaz et al.

moveout (Yilmaz, 2001). For events with nonhyperbolic moveout that calls for depth migration, each of the resulting common-offset volumes following the applications of NMO and 3-D DMO corrections would not be a representation of a 3-D zero-offset wavefield. As such, and again in theory, decoupling of common-offset data volumes is not acceptable for prestack depth migration. Despite this theoretical mandate, it is surprisingly pleasant to see that indeed, in practice, the common-offset strategy for 3-D prestack time migration sometimes can also be applicable to 3-D prestack depth migration. Whereas following NMO and 3-D DMO corrections, each of the common-offset volumes is time-migrated using a 3-D rms velocity field, depth migration of the common-offset volumes is done using a 3-D velocity-depth model. There are two reasons why you may choose to apply the common-offset strategy to 3-D prestack depth migration. (1) You may wish to use a 3-D zero-offset depth migration

algorithm other than the Kirchhoff summation, and thus avoid the troublesome task of computing 3-D nonzero-offset traveltimes needed for the latter. (2) For line output or for selected image gathers, Kirchhoff summation is the appropriate algorithm. For volume output from 3-D prestack depth migration, other algorithms can be more efficiently applied to commonoffset data. For instance, a 3-D zero-offset depth migration based on the McLellan transform (Hale, 1991) or phase-shiftplus-correction (PSPC) (Kosloff and Kessler, 1987) can be used to migrate each of the 3-D DMO-corrected common-offset data volumes, independently. The final image volume from 3-D prestack depth migration is then produced by summing the resulting common-offset image volumes (Figures 12c, d). A choice of algorithms for 3-D prestack depth migration to create an image volume is summarized as follows: 1) You may wish to first apply azimuth-moveout (AMO) correction (Biondi et al., 1998) to regularize the pres tack data, then use a frequency-wavenumber algorithm (Biondi and Palacharla, 1994) to perform the 3-D prestack depth migration. The method involves design and application of a wavefield extrapolation operator on the whole of the 3-D prestack data without decoupling the common-offset volumes. 2) Alternatively, you may wish to perform 3-D prestack depth migration using the Kirchhoff summation technique (Vidale, 1988; Podvin and Lecomte, 1991; Reshef, 1991; Lecomte, 1999). The Kirchhoff summation implicitly handles geometry irregularities associated with the prestack data. An integral part of the Kirchhoff summation technique is calculation of traveltimes associated with source and receiver locations at the surface and reflection points in the subsurface (Reshef and Kosloff, 1986; Moser, 1991; Meshbey et al., 1993; Nichols, 1996; Vesnaver, 1996; Sethian and Popovici, 1999). 3) A third alternative, which is suggested in this paper, is based on depth migration of the decoupled 3-D commonoffset volumes of data that have been corrected for geometry irregularities and reduced to zero offset by way of 3-D DMO correction (Yilmaz, 2001). Again, the preferred choice for the algorithm to perform the 3-D common-offset depth migration is based on wave extrapolation. The common-offset strategy for depth migration will not work strictly for a case of complex overburden structure that gives rise to complex, nonhyperbolic moveout. A general guideline for the choice of strategy for 3-D prestack depth migration is as follows:

FIG. 11. Two image gathers and the residual moveout semblance spectra computed from them before [(a) and (b)] and after [(c) and (d)] residual moveout corrections for model update. The arrows indicate the vertical axis in black that corresponds to the zero residual moveout. The white curves in the semblance panels represent the interval velocity function at the analysis location before the application of the residual moveout corrections for model update.

1) If you are dealing with a low-relief structure and moderate lateral velocity variations, use 3-D prestack depth migration for conducting image-gather analysis for model verification and update, but not necessarily for imaging in depth. For the purpose of model verification and update, you need to work with only a sparse grid of image gathers; the Kirchhoff summation is the suitable algorithm to produce such a set of image gathers without creating a whole image volume. As for imaging in depth, you may be content with 3-D poststack depth migration; a finitedifference or frequency-wavenumber algorithm that uses wave extrapolation are the suitable algorithms.

A Unified 3-D Seismic Workflow

2) If you are dealing with a complex structure and moderateto- strong lateral velocity variations, for imaging use a 3-D prestack depth migration algorithm based on decoupling of common-offset volumes of data, as suggested in this paper. 3) Finally, if you are dealing with a complex overburden structure, for imaging the use of a 3-D prestack depth mi-

1709

gration algorithm without decoupling the common-offset data is imperative. STRUCTURAL AND STRATIGRAPHIC INTERPRETATION

We now interpret the image volume from 3-D prestack depth migration, a structural intepretation based on picking a set of

FIG. 12. Two cross-sections from the structurally consistent, final 3-D interval velocity volume following the model update based on residual moveout analysis as in Figure 11: (a) an in-line section and (b) a cross-line section. Also, two cross-sections from the image volume derived from 3-D prestack depth migration using the final model represented by the 3-D interval velocity field as in (a) and (b): (c) an in-line section and (d) a cross-line section.

1710

Yilmaz et al.

depth horizons and stratigraphic interpretation based on amplitude manipulation (phase 5 as in Figure 1). Figure 13 shows six depth horizons that coincide with the geologic markers in the area. The top three horizons, D H 1, D H 2, and D H 3, were delineated largely by using the seed detection technique that involves connecting neighboring voxels (3-D pixels) with am-

plitudes that are within a specified range. The bottom three horizons, D H 4, D H 5, and D H 6, were delineated by line-based interpretation that involves creating a set of horizon strands along selected inlines and crosslines. The effect of extensional tectonism is evident on the top three depth horizons. Note in particular the series of faults

FIG. 13. Depth horizons interpreted from the image volume derived from 3-D prestack depth migration as in Figures 12c, d: (a) is the shallowest horizon and (f) is the deepest horizon.

A Unified 3-D Seismic Workflow

across the lowlands indicated by the green. Examine the depth horizons in reverse order, starting with D H 6 and ending with D H 1, to reconstruct the structural evolution of the area. In the beginning, there was an erosional surface cut by a deep channel, as seen on horizon D H 6. There may have been a period of compressional tectonism that gave rise to the gently folded surfaces in the lowlands portion of horizons D H 5 and D H 4. Then, a reversal from compressional tectonism to extensional tectonism began to cause subsidence of the basin, which in turn gave rise to the formation of the series of faults parallel to the ancient shoreline. Throughout its geologic times, the highlands indicated by the orange on all horizons remained as a hiatus until the geologic time approximately coincident with the age of the overburden above horizon D H 1. Now split the image volume into subvolumes that represent individual depositional units as shown in Figures 14a, b. Observe the extensional fault patterns on the top surfaces of these subvolumes. To investigate the interior of each of these depositional units, cut them into thin slices. Then, apply transparency to each of the thin slices to identify depositional features of interest. Shown in Figure 14c is the map view of a thin slice from the depositional unit in Figure 14a that exhibits the stream traverse disrupted by the faults. As indicated in Figure 1, the seismic pathway we followed from phases 1 through 5 in this case study should lead us to constructing a model for the reservoir. Begin with a detailed delineation of the top and base of the deltaic sequence that represents the reservoir unit shown in Figure 15a. Structural interpretation of the image volume derived from 3-D prestack depth migration (Figures 12c, d) yields the top and base surfaces shown in Figure 15b and 15c. Then, define the interior geometry of the reservoir by splitting the deltaic sequence bounded by the two surfaces into a set of thin layers, as shown in Figure 15a. Cross-sections of the reservoir unit along selected in-line traverses after sublayering are shown in Figure 16. Next, extract the subvolume associated with the reservoir unit from within the image volume derived from 3-D prestack depth migration (Figure 17). You may want to use the subvolume to estimate a set of seismic attributes. The sublayers (Figure 16) then are populated (based on any available well data and the results of the seismic attribute analysis, such as AVO and acoustic impedance) by the petrophysical properties of the reservoir unit, such as, porosity, permemability, and fluid saturation, all of which are allowed to vary laterally within each of the sublayers.

1711

is to take a set of stacking or DMO velocity functions picked at analysis locations, edit for any outliers, and interpolate them with some lateral and vertical smoothing. The key to an accurate interval velocity field from Dix conversion, however, is to estimate the rms velocity field, not from unmigrated data, but from prestack time-migrated data. Hence, implicit to the model building strategy advocated in this paper is the requirement for doing prestack time migration for earth modeling prior to prestack depth migration for earth imaging in depth.

ASPECTS OF THE PROPOSED WORKFLOW—A SUMMARY

We have demonstrated a unified workflow that involves processing, inversion, and interpretation of 3-D seismic data. The seismic workflow described in this paper includes earth modeling and imaging in depth with unique aspects. An interval velocity field for depth migration derived from Dix conversion of an rms velocity field often is judged to be inaccurate and, consequently, can require excessive number of iterations for model update to attain a final velocity-depth model with acceptable accuracy. The most likely reason for failing to derive an accurate interval velocity field from Dix conversion to begin with is that the input rms velocity field usually is based on velocity estimation from unmigrated data. For instance, a common approach to derive an rms velocity field

FIG. 14. Image subvolumes extracted from the total image volume derived from 3-D prestack depth migration as in Figures 12c, d: (a) the subvolume associated with the layer bounded by horizon D H 1 on top and D H 2 at the bottom, (b) the subvolume associated with the layer bounded by horizon D H 2 on top and D H 3 at the bottom, and (c) map view of a thin slice extracted from the image subvolume in (a) after the removal of opacity. Refer to Figure 13 for the geometry of the depth horizons themselves.

1712

Yilmaz et al.

A common strategy for estimating rms velocities associated with migrated data is based on DMO correction combined with common-offset migration. Time migration of seismic data requires an rms velocity field associated with events in their migrated positions. This means that you need to migrate the data to estimate the rms velocity field, which is what you need to migrate the data. This paper offers an interpretive approach to circumvent this circuitous argument. Specifically, we have presented a procedure for 3-D migration velocity analysis to

derive the 3-D rms velocity field based on velocity-independent 3-D prestack time migration. We want to obtain the best estimate of a structurally consistent initial earth model in depth where layer velocities are based on rms velocities associated with migrated data so as to minimize the work required to update the model and thus obtain a final earth model in depth. To achieve structural consistency in model building, and thus preserve the vertical and lateral velocity variations in the 3-D interval velocity field derived from Dix conversion, we make use of a set of phantom

FIG. 16. Cross-sections of the reservoir unit represented by a set of thin layers as in Figure 15a along a set of in-line traverses (a) through (i).

FIG. 15. (a) A portion of an in-line section from the image volume derived from 3-D prestack depth migration with the reservoir zone interpreted by a set of thin layers; (b) top and (c) base of the reservoir unit.

FIG. 17. Image subvolume associated with the reservoir unit extracted from the total image volume derived from 3-D prestack depth migration as in Figures 12c, d.

A Unified 3-D Seismic Workflow

depth horizons by subdividing each of the layers bounded by the depth horizons interpreted from the depth-migrated volume of data into a set of sublayers. During model updating, whether by residual moveout analysis or reflection tomography, it suffices to perturb layer velocities only, while retaining the reflector geometries governed by the principal and phantom horizons the same from one iteration to the next. Once the final velocity-depth model is obtained, the final reflector geometries are recovered by performing structural interpretation of the image volume from 3-D prestack depth migration. Albeit that in theory decoupling of common-offset data volumes is not acceptable for prestack depth migration, we have demonstrated that the common-offset strategy for 3-D prestack time migration sometimes can also be applicable to 3-D prestack depth migration. Following NMO and 3-D DMO corrections, each of the common-offset volumes is time migrated using a 3-D rms velocity field, whereas depth migration of the common-offset volumes is done using a 3-D velocitydepth model. ACKNOWLEDGEMENTS

We greatly appreciate the support given to us at various stages in the project by our colleagues from Paradigm Geophysical, especially by Ferudun Kilic and Fugen Zhou. We thank the Shengli Oil Field Division of SINOPEC for granting the permission to work with the 3-D seismic data set used in this case study and for publishing the results of the analysis based on the unified seismic workflow. We also express our sincere gratitudes to Paradigm Geophysical for supporting the project. REFERENCES 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. Beasley, C. J., and Mobley, E., 1998, Spatial dealiasing of 3-D DMO: The Leading Edge, 17, 1590–1594.

1713

Biondi, B., Fomel, S., and Chemingui, N., 1998, Azimuth moveout for 3-D prestack imaging: Geophysics, 63, 574–588. Biondi, B., and Palacharla, G., 1994, 3-D migration using rotated McClellan filters: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1278–1281. Deregowski, S. M., 1990, Common-offset migrations and velocity analysis: First Break, 8, 225–234. Fowler, P., 1984, Velocity-independent imaging of seismic reflectors: 54th Ann. Internat. Mtg., Soc. Explor. Geophys., Expanded Abstracts, 383–385. Hale, D., 1991, 3-D depth migration via McClellan transforms: Geophysics, 56, 1778–1785. Kosloff, D., and Kessler, D., 1987, Accurate depth migration by a generalized phase-shift method: Geophysics, 52, 1074–1084. Lecomte, I., 1999, Local and controlled prestack depth migration in complex areas: Geophys. Prosp., 47, 799–818. Marcoux, M. O., Godfrey, R. J., and Notfors, C. D., 1987, Migration for optimum velocity evaluation and stacking: Presented at the 49th Ann. Mtg., Eur. Assn. Expl. Geophys. Meshbey, V., Kosloff, D., Ragoza, Y., Meshbey, O., Egozy, U., and Cozzens, J., 1993, A method for computing traveltimes for an arbitrary velocity model: Presented at the 45th Ann. Mtg., Eur. Assn. Geosci. Eng. Moser, T. J., 1991, Shortest-path calculation of seismic rays: Geophysics, 56, 59–67. Nichols, D. E., 1996, Maximum-energy traveltimes calculated in the seismic frequency band: Geophysics, 61, 253–263. Podvin, P., and Lecomte, I., 1991, Finite-difference computation of traveltimes for velocity-depth models with strong velocity contrast across layer boundaries—A massively parallel approach: Geophys. J. Internat., 105, 271–284. Reshef, M., 1991, Prestack depth imaging of three-dimensional shot gathers: Geophysics, 56, 1158–1163. Reshef, M., and Kosloff, D., 1986, Migration of common-shot gathers: Geophysics, 51, 324–331. Ronen, S., 1994, Handling irregular geometry: Equalized DMO and beyond: 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1545–1548. Ronen, S., Nichols, D., Bale, R., and Ferber, R., 1995, Dealiasing DMO: Good-pass, bad-pass and unconstrained: 66th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 743–746. Sethian, J. A., and Popovici, A. M., 1999, 3-D traveltime computation using the fast marching method: Geophysics, 64, 516–523. Spitz, S., 1991, Seismic trace interpolation in the f -x domain: Geophysics, 56, 785–794. Vesnaver, A., 1996, Ray tracing based on Fermat’s principle in irregular grids: Geophys. Prosp., 44, 741–760. Vidale, J., 1988, Finite-difference calculation of traveltimes: Bull. Seis. Soc. Am., 78, 2026–2076. Yilmaz, O., 2001, Seismic data analysis—Processing, inversion, and interpretation of seismic data: Soc. Expl. Geophys.