04 Migration

Prévia do material em texto

4 Migration
• Introduction • Exploding Reflectors •Migration Strategies •Migration Algorithms •Migration Parameters • As-
pects of Input data • Migration Velocities •Migration Principles • Kirchhoff Migration • Diffraction Summation
• Amplitude and Phase Factors • Kirchhoff Summation • Finite-Difference Migration • Downward Continuation •
Differencing Schemes • Rational Approximations for Implicit Schemes • Reverse Time Migration • Frequency-Space
Implicit Schemes • Frequency-Space Explicit Schemes • Frequency-Wavenumber Migration • Phase-Shift Migra-
tion • Stolt Migration • Summary of Domains of Migration Algorithms • Kirchhoff Migration in Practice •
Aperture Width • Maximum Dip to Migrate • Velocity Errors • Finite-Difference Migration in Practice •
Depth Step Size • Velocity Errors • Cascaded Migration • Reverse Time Migration • Frequency-Space Migra-
tion in Practice • Steep-Dip Implicit Methods • Depth Step Size • Velocity Errors • Steep-Dip Explicit Methods
• Dip Limits of Extrapolation Filters • Velocity Errors • Frequency-Wavenumber Migration in Practice •
Maximum Dip to Migrate • Depth Step Size • Velocity Errors • Stolt Stretch Factor • Wraparound • Residual Mi-
gration • Further Aspects of Migration in Practice • Migration and Spatial Aliasing • Migration and Random
Noise • Migration and Line Length • Migration from Topography • Exercises • Appendix D: Mathematical
Foundation of Migration • Wavefield Extrapolation and Migration • Stationary Phase Approximations • The
Parabolic Approximation • Frequency-Space Implicit Schemes • Stable Explicit Extrapolation • Optimum Depth
Step • Frequency-Wavenumber Migration • Residual Migration • References
4.0 INTRODUCTION
Migration moves dipping reflections to their true sub-
surface positions and collapses diffractions, thus increas-
ing spatial resolution and yielding a seismic image of
the subsurface. Figure 4.0-1 shows a CMP-stacked sec-
tion before and after migration. The stacked section in-
dicates the presence of a salt dome flanked by gently
dipping strata. Figure 4.0-1 also shows a sketch of two
prominent features — the diffraction hyperbola D that
originates at the tip of the salt dome, and the reflection
B off the flank of the salt dome. After migration, note
that the diffraction has collapsed to its apex P and the
dipping event has moved to a subsurface location A,
which is at or near the salt dome flank. In contrast, re-
flections associated with the gently dipping strata have
moved little after migration.
Figure 4.0-2 is an example with a different type
of structural feature. The stack contains a zone of near-
horizontal reflections down to 1 s. After migration, these
events are virtually unchanged. Note the prominent un-
conformity that represents an ancient erosional surface
just below 1 s. On the stacked section, the unconfor-
mity appears complex, while on the migrated section,
it becomes interpretable. The bowties on the stacked
section are untied and turned into synclines on the mi-
464 Seismic Data Analysis
FIG. 4.0-1. A CMP stack (a) before, (b) after migration; (c) sketch of a prominent diffraction D and a dipping event before
(B) and after (A) migration. Migration moves the dipping event B to its assumed true subsurface position A and collapses
the diffraction D to its apex P. The dotted line indicates the boundary of a salt dome.
FIG. 4.0-2. A CMP stack (a) before, (b) after migration.
Migration unties the bowties on the stacked section and
turns them into synclines (Taner and Koehler, 1977).
grated section. The deeper event in the neighborhood
of 3 s is the multiple associated with the unconformity
above. When treated as a primary and migrated with
the primary velocity, it is overmigrated.
Figure 4.0-3a shows a stacked section that con-
tains fault-plane reflections conflicting with the shal-
low gently-dipping reflections. Note the accurate posi-
tioning of the fault planes and delineation of the fault
blocks on the migrated section in Figure 4.0-3b. From
the three examples shown in Figures 4.0-1, 4.0-2, and
4.0-3, note that migration moves dipping events in the
updip direction and collapses diffractions, thus enabling
us to delineate faults while retaining horizontal events
in their original positions.
The goal of migration is to make the stacked sec-
tion appear similar to the geologic cross-section in depth
along a seismic traverse. The migrated section, however,
commonly is displayed in time. One reason for this is
that velocity estimation based on seismic and other data
always is limited in accuracy. Therefore, depth conver-
sion is not completely accurate. Another reason is that
interpreters prefer to evaluate the validity of migrated
sections by comparing them to the unmigrated data.
Therefore, it is preferable to have both sections dis-
played in time. The migration process that produces
a migrated time section is called time migration. Time
migration, the main theme of Chapter 4, is appropriate
as long as lateral velocity variations are mild to moder-
ate.
When the lateral velocity gradients are significant,
time migration does not produce the true subsurface im-
age. Instead, we need to use depth migration, the output
Migration 465
FIG. 4.0-3. A CMP stack (a) before, (b) after migration. Migration collapses subtle diffractions associated with the growth
faults, moves the fault-plane reflections to the fault positions, and thus makes detailed structural interpretation easier.
466 Seismic Data Analysis
FIG. 4.0-4. A CMP stack (a) before, (b) after time migration. Time migration is adequate for accurate imaging of the
top-salt boundary, whereas depth migration is imperative for accurate imaging of the base-salt boundary (B).
Migration 467
of which is a depth section. Consider the data from an
area with intense salt tectonics in Figure 4.0-4. Time
migration has produced an acceptable image of the re-
gion above the salt. However, note the crossing of events
that is a manifestation of overmigration of the reflection
associated with the base-salt boundary (denoted by B
in Figure 4.0-4b). The improper migration is the result
of inadequate treatment by the time migration of the
effects of severe raypath bendings at the top-salt bound-
ary caused by the strong velocity contrast between the
salt layer and the overlying rocks.
Complex structures associated with salt diapirism,
overthrust tectonics and irregular water-bottom topog-
raphy usually are three dimensional (3-D) in character.
A stacked section really is the seismic response of a 3-
D subsurface on a two-dimensional (2-D) plane of pro-
file. Therefore, 2-D migration is not completely valid
for 3-D data from areas with complex 3-D structures.
Figure 4.0-5a is an inline stacked section from a land
3-D survey. Figure 4.0-5b is a 2-D migration of this
section, while Figure 4.0-5c is the same section after
3-D migration of the entire 3-D survey data. In particu-
lar, note the significant difference in the imaging of the
top of salt T and base of salt B. In 2-D migration, we
assume that the stacked section does not contain any
energy that comes from outside the plane of recording
(sideswipe). Three-dimensional imaging of the subsur-
face is discussed in Section 7.3.
Exploding Reflectors
When a stacked section is migrated, we use the migra-
tion theory applicable to data recorded with a coin-
cident source and receiver (zero-offset). To develop a
conceptual framework for discussing migration of zero-
offset data, we now examine two types of recording
schemes.
A zero-offset section is recorded by moving a sin-
gle source and a single receiver along the line with no
separation between them (Figure 4.0-6). The recorded
energyfollows raypaths that are normal incidence to
reflecting interfaces. This recording geometry obviously
is not realizable in practice.
Now consider an alternative geometry (Figure 4.0-
6) that will produce the same seismic section. Imag-
ine exploding sources that are located along the reflect-
ing interfaces (Loewenthal et al., 1976). Also, consider
one receiver located on the surface at each CMP lo-
cation along the line. The sources explode in unison
and send out waves that propagate upward. The waves
are recorded by the receivers at the surface. The earth
model described by this experiment is referred to as the
exploding reflectors model.
FIG. 4.0-5. A 2-D CMP stack (a) represents a 2-D cross-
section of a 3-D wavefield. Thus, it can contain energy from
outside the plane of the 2-D line traverse. A 2-D migra-
tion (b) is inadequate when this kind of energy is present
on the 2-D CMP-stacked section. (c) Clear imaging of the
salt structure requires both 3-D data collection and 3-D mi-
gration (Section 7.3). (Data courtesy Nederlandse Aardolie
Maatschappij B.V.)
468 Seismic Data Analysis
FIG. 4.0-6. Geometry of zero-offset recording (left), and hypothetical simulation of the zero-offset experiment using exploding
reflectors (right) (Claerbout, 1985).
The seismic section that results from the exploding
reflectors model is largely equivalent to the zero-offset
section, with one important distinction. The zero-offset
section is recorded as two-way traveltime (from source
to reflection point to receiver), while the exploding re-
flectors model is recorded as one-way traveltime (from
the reflection point at which the source is located to
the receiver). To make the sections compatible, we can
imagine that the velocity of propagation is half the true
medium velocity for the exploding reflectors model.
The equivalence between the zero-offset section and
the exploding reflectors model is not quite exact, par-
ticularly in the presence of strong lateral velocity vari-
ations (Kjartansson and Rocca, 1979).
These concepts now are applied to the velocity-
depth model in Figure 4.0-7. Consider source-receiver
pairs placed along the earth’s surface at every tenth
midpoint. In this case, a zero-offset section is mod-
eled. At midpoint 130, five different arrivals are asso-
ciated with rays that are normal incidence to the first
interface. Alternatively, imagine receivers placed along
the earth’s surface at every tenth midpoint and sources
placed along the interface where the rays emerge at the
right angle to the interface (equivalent to the normal-
incidence rays of the zero-offset section). In the latter
case, the velocities indicated in Figure 4.0-7 must be
halved to match the time axis with that associated with
the zero-offset section.
The interface can be sampled more densely by plac-
ing receivers and sources at closer spacing (Figure 4.0-
8a). The next deeper interface can be modeled; that is,
FIG. 4.0-7. A velocity-depth model (top) and the zero-
offset traveltime response (bottom) of the water-bottom re-
flector. Shown also are the normal-incidence rays used to
compute the zero-offset traveltime trajectory. Note the five
arrivals A, B, C, D, and E at CMP 130 — all from the water
bottom.
Migration 469
FIG. 4.0-8. Exploding-reflector modeling of zero-offset traveltimes associated with (a) a water bottom, (b) a flat, and (c)
a dipping reflector. (d) The superposition of the normal-incidence traveltime responses in (a), (b), and (c). Shown on the
velocity-depth models in the left-hand column are the normal-incidence rays used to compute the traveltime trajectories.
The time sections shown on the right-hand column are equivalent to a zero-offset traveltime section with the vertical axis in
two-way time.
470 Seismic Data Analysis
the traveltime trajectory can be computed by placing
sources along this interface and leaving the receivers
where they were on the surface (Figure 4.0-8b). Finally,
the same experiment can be repeated for the third inter-
face (Figure 4.0-8c). To derive the composite response
from the velocity-depth model in Figure 4.0-8d (the left-
hand column), individual responses shown in Figures
4.0-8a, 4.0-8b and 4.0-8c from each interface are su-
perimposed. The result is shown in Figure 4.0-8d (the
right-hand column). We can imagine that sources were
placed at all three interfaces and turned on simultane-
ously. Such an experiment would cause the rays emerg-
ing from the three interfaces to be recorded at receivers
placed on the surface, along the line.
Actually, Figure 4.0-8d (the right-hand column)
represents the modeled zero-offset traveltime section.
Seismic wavefields, however, are represented not only
by wave traveltimes but also by wave amplitudes. Figure
4.0-9a shows the modeled zero-offset wavefield section
based on the same velocity-depth model in Figure 4.0-
8d (the left-hand column). The shallow complex inter-
face (horizon 1 in Figure 4.0-8a) caused the complicated
response of the two simple interfaces (horizons 2 and 3)
in this zero-offset traveltime section.
How valid is the assumption that a stacked sec-
tion is equivalent to a zero-offset section? The conven-
tional CMP recording geometry provides the wavefield
at nonzero offsets. During processing, we collapse the
offset axis by stacking the data onto the midpoint-time
plane at zero offset. For CMP stacking, we normally as-
sume hyperbolic moveout. Figure 4.0-10 shows selected
CMP gathers modeled from the velocity-depth model
in Figure 4.0-8d (the left-hand column). Because of the
presence of strong lateral velocity variations, the hy-
perbolic assumption may not be appropriate for some
reflections on some CMP gathers (Figure 4.0-10a); how-
ever, it may be valid for others (Figure 4.0-10b). We ob-
tain a stacked section (Figure 4.0-9b) that resembles the
zero-offset section (Figure 4.0-9a) to the extent that the
hyperbolic moveout assumption is valid. The assump-
tion that a conventional stacked section is equivalent to
a zero-offset section also is violated to varying degrees
in the presence of strong multiples and conflicting dips
with different stacking velocities (Chapter 5). While mi-
gration of unstacked data is discussed in Chapter 5, our
main focus in this chapter is on migration after stack.
Migration Strategies
In practice, migration of seismic data requires decision
making with regards to:
(a) an appropriate migration strategy,
(b) a migration algorithm compatible with the strat-
egy,
(c) appropriate parameters for the algorithm,
(d) issues concerning the input data, and
(e) migration velocities.
Migration strategies include:
(a) 2-D versus 3-D migration,
(b) post- versus prestack migration, and
(c) time versus depth migration.
The spectrum of migration strategies extend from 2-D
poststack time migration to 3-D prestack depth migra-
tion. Depending on the nature of the subsurface geol-
ogy, any other in-between combination can be selected.
In practice, 2-D/3-D poststack time migration is used
most often for a good reason — it is the least sensitive
to velocity errors, and it often yields results acceptable
for a reliable interpretation. Table 4-1 is an overview of
different migration strategies applied to different types
of seismic data (2-D, 3-D, stacked, and unstacked).
Choice of an appropriate migration strategy re-
quires input from the interpreter as to the structural
geology and stratigraphy in an area. Dipping events on a
stacked section call for time migration. Conflicting dips
with different stacking velocities is one case in which a
conventional stacked section differs from a zero-offset
section. Thus, strictly speaking,poststack migration
Table 4-1. Migration strategies.
Case Migration
dipping events time migration
conflicting dips with prestack migration
different stacking velocities
3-D behavior of 3-D migration
fault planes and salt flanks
Case Migration
strong lateral velocity depth migration
variations associated with
complex overburden structures
complex nonhyperbolic moveout prestack migration
3-D structures 3-D migration
Migration 471
FIG. 4.0-9. (a) The zero-offset wavefield section equivalent
to the zero-offset traveltime section in Figure 4.0-8d (the
right-hand column); (b) the CMP stack generated from the
CMP gathers as in Figure 4.0-10. (Modeling by Deregowski
and Barley, 1981.)
which assumes that the stacked section is equivalent to
a zero-offset section is not valid to handle the case of
conflicting dips. Instead, one needs to do prestack time
migration.
Conflicting dips often are associated with salt
flanks and fault planes, which can have 3-D characteris-
tics. This then requires 3-D prestack time migration. In
Section 5.3, we shall discuss a practical alternative to
2-D prestack time migration strategies. The alternative
sequence includes the application of normal-moveout
(NMO) correction using velocities appropriate for flat
events followed by 2-D dip moveout correction (DMO)
to correct for the dip and source-receiver azimuth ef-
fects on stacking velocities. As a result, conflicting dips
are preserved during stacking, and thus, imaging can be
deferred until after stacking using 2-D poststack time
migration strategies. This series of processing steps is
largely equivalent to 2-D prestack time migration and
results often are comparable. The same workflow also is
applicable to 3-D prestack time migration (Section 7.4).
Accurate imaging of targets beneath complex
structures with strong lateral velocity variations re-
quires depth migration. Aside from the problem of con-
flicting dips with different stacking velocities, strong lat-
eral velocity variations associated with complex over-
burden structures usually cause conventional stacking
based on the hyperbolic moveout assumption to fail.
Therefore, a case of complex overburden structures calls
for depth migration before stacking the data.
Furthermore, complex overburden structures, en-
countered in areas with salt tectonics, overthrust tec-
tonics and irregular water-bottom topographies can
often exhibit 3-D characteristics. Thus, imaging such
structures may require 3-D prestack depth migration.
Field surveys are designed such that line orienta-
tions are, as much as possible, along the dominant strike
and dip directions, so as to minimize 3-D effects. Un-
der these circumstances, the 2-D assumption for migra-
tion can be acceptable. However, if the subsurface has
a truly 3-D geometry, without a dominant dip or strike
direction, then it is imperative to do 3-D migration of
3-D data. In such cases, 2-D migration (whether post-
stack or prestack, time or depth) can lead to potential
problems in interpretation.
A practical alternative to 2-D prestack depth mi-
gration can be a prestack layer replacement to cor-
rect for the complex nonhyperbolic moveout followed
by time migration after stack. This, however, is appli-
cable to situations involving a single overburden layer,
such as irregular water-bottom topography for it to be
reasonably practical.
Migration Algorithms
The one-way-in-depth scalar wave equation is the ba-
sis for common migration algorithms. These algorithms
do not explicitly model multiple reflections, converted
waves, surface waves, or noise. Any such energy present
in data input to migration is treated as primary re-
flections. Migration algorithms can be classified under
three main categories:
(a) those that are based on the integral solution to the
scalar wave equation,
(b) those that are based on the finite-difference solu-
tions, and
(c) those that are based on frequency-wavenumber im-
plementations.
Whatever the algorithm, it should desirably:
(a) handle steep dips with sufficient accuracy,
(b) handle lateral and vertical velocity variations, and
(c) be implemented, efficiently.
Figure 4.0-11 is a migrated CMP stacked section
with a major unconformity. The undermigration — in-
complete imaging of the unconformity, is not because of
472 Seismic Data Analysis
FIG. 4.0-10. Selected CMP gathers modeled from the velocity-depth model in Figure 4.0-8d (the left-hand column). (a)
Gathers from the complex part of the velocity-depth model, (b) gathers from the simpler part of the velocity-depth model.
CMP locations are indicated in Figure 4.0-8d. The CMP stack is shown in Figure 4.0-9b. (Modeling by Deregowski and Barley,
1981.)
erroneously too low velocities. Although we should al-
ways be aware of velocity errors when migrating seismic
data, the undermigration in Figure 4.0-11 is the result
of using a dip-limited algorithm. By using a steep-dip
algorithm, we can achieve a more accurate imaging of
the unconformity (Figure I-9).
The three principle migration techniques are dis-
cussed in this chapter in their historical order of devel-
opment as outlined below. The first migration technique
developed was the semicircle superposition method that
was used before the age of computers. Then came
the diffraction-summation technique, which is based on
summing the seismic amplitudes along a diffraction hy-
perbola whose curvature is governed by the medium ve-
locity. The Kirchhoff summation technique introduced
later (Schneider, 1978), but actually in use earlier, basi-
cally is the same as the diffraction summation technique
with added amplitude and phase corrections applied to
the data before summation. These corrections make the
summation consistent with the wave equation in that
they account for spherical spreading (Section 1.4), the
obliquity factor (angle-dependency of amplitudes), and
the phase shift inherent in Huygens’ secondary sources
(Section 4.1).
Another migration technique (Claerbout and Do-
herty, 1972) is based on the idea that a stacked section
can be modeled as an upcoming zero-offset wavefield
generated by exploding reflectors. Using the explod-
ing reflectors model, migration can be conceptualized
as consisting of wavefield extrapolation in the form of
downward continuation followed by imaging. To under-
stand imaging, consider the shape of a wavefield at ob-
servation time t = 0 generated by an exploding reflec-
tor. Since no time has elapsed and, thus, no propagation
has occurred, the wavefront shape must be the same as
the reflector shape that generated the wavefront. The
fact that the wavefront shape at t = 0 corresponds to
the reflector shape is called the imaging principle. To
define the reflector geometry from a wavefield recorded
at the surface, we only need to extrapolate the wave-
field back in depth then monitor the energy arriving at
t = 0. The reflector shape at any particular extrapola-
Migration 473
FIG. 4.0-11. A CMP stack (a) before, and (b) after migration. Note the undermigration of the unconformity event (U)
caused by the use of a dip-limited algorithm.
tion depth directly corresponds to the wavefront shape
at t = 0.
Downward continuation of wavefields can be im-
plemented conveniently using finite-difference solutions
to the scalar wave equation. Migration methods based
on such implementations are called finite-difference mi-
gration. Many different differencing schemes applied to
the differential operators in the scalar wave equation
exist both in time-space and frequency-space domains.
Claerbout (1985) provides a comprehensive theoretical
foundation of finite-difference migration and itspracti-
cal aspects.
After the developments on Kirchhoff summation
and finite-difference migrations, Stolt (1978) introduced
migration by Fourier transform. This method involves
a coordinate transformation from frequency (the trans-
form variable associated with the input time axis) to
vertical wavenumber axis (the transform variable asso-
ciated with the output depth axis), while keeping the
horizontal wavenumber unchanged. The Stolt method
is based on a constant-velocity assumption. However,
Stolt modified his method by introducing stretching
in the time direction to handle the types of velocity
variations for which time migration is acceptable. Stolt
and Benson (1986) combine theory with practice in the
field of migration with an emphasis on the frequency-
wavenumber methods.
Another frequency-wavenumber migration is the
phase-shift method (Gazdag, 1978). This method is
based on the equivalence of downward continuation to a
phase shift in the frequency-wavenumber domain. The
imaging principle is invoked by summing over the fre-
quency components of the extrapolated wavefield at
each depth step to obtain the image at t = 0.
A reason for the wide range of migration algorithms
used in the industry today is that none of the algorithms
fully meets the important criteria of handling all dips
and velocity variations while still being cost-effective.
Migration algorithms based on the integral solu-
tion to the scalar wave equation, commonly known as
Kirchhoff migration, can handle all dips up to 90 de-
grees, but they can be cumbersome in handling lateral
velocity variations.
Finite-difference algorithms can handle all types of
velocity variations, but they have different degrees of
dip approximations. Furthermore, differencing schemes,
if carelessly designed, can severely degrade the intended
dip approximation.
474 Seismic Data Analysis
FIG. 4.0-12. A CMP stack (a) before, and (b) after migration. Lack of any event to the right of the dotted line on the
migrated section is a result of the finite line length.
Finally, frequency-wavenumber algorithms have
limited ability in handling velocity variations, partic-
ularly in the lateral direction. As a result of limi-
tations of the three main categories of migration al-
gorithms — integral, finite-difference, and frequency-
wavenumber methods, migration software has expanded
further to additional extensions and combinations of the
basic algorithms. Residual migration — phase-shift or
constant-velocity Stolt migration followed by the appli-
cation of a dip-limited algorithm is one example.
Migration Parameters
After deciding on the migration strategy and the ap-
propriate algorithm, the analyst then needs to decide
on the migration parameters. Migration aperture width
is the critical parameter in Kirchhoff migration. A small
aperture causes removal of steep dips; it generates spu-
rious horizontal events and organizes the random noise
uncorrelated from trace to trace.
Depth step size in downward continuation is the
critical parameter in finite-difference methods. An opti-
mum depth step size is the largest depth step with the
minimum tolerable phase errors. It depends on tempo-
ral and spatial samplings, dip, velocity, and frequency.
It also depends on the type of differencing scheme used
in the algorithm.
Finally, the stretch factor is the critical parameter
for Stolt migration. A constant-velocity medium implies
a stretch factor of 1. In general, the larger the vertical
velocity gradient, the smaller the stretch factor needs
to be.
Migration 475
Aspects of Input Data
When migrating seismic data, one needs to be con-
cerned with various aspects of the input data set itself:
(a) line length or areal extent,
(b) signal-to-noise ratio, and
(c) spatial aliasing.
The line length must be sufficient to allow a steeply
dipping event to migrate to its true subsurface location.
It is a fatal error to record short profiles in areas with
complex geology. Also, for 3-D migration, the surface
areal extent of a 3-D survey is almost always larger than
the target subsurface areal extent.
Random noise at late times on a stacked section,
when migrated, potentially can be hazardous for shal-
lower data. One may have to compromise on migration
aperture for deep data to prevent this problem to occur.
Trace spacing must be small enough to prevent
spatial aliasing of steep dips at high frequencies. Al-
though this is not an issue for modern prestack data, a
coarse shot-receiver spacing can degrade the fidelity of
prestack migration severely. Old data and 3-D marine
data in the crossline direction often are trace interpo-
lated prior to poststack migration.
Figure 4.0-12 is a CMP stacked section before and
after migration. From an interpretation viewpoint, the
reliable part of this migrated section is confined to the
upper central part. Lack of any reflection energy to the
right of the dotted line does not mean that there is a
structural discontinuity there. It only means that the
reflections associated with the imbricate structure have
been migrated in the updip direction from right to left.
As a result, a zone with no reflectors to the right of the
dotted line is left behind because of the truncation of
the wavefield represented by the right-hand edge of the
stacked section. The deeper part is useless, because the
noise dominates the section.
Migration Velocities
Horizontal displacement during migration is propor-
tional to migration velocity squared (equation 4-1).
Since velocities generally increase with depth, errors in
migration are usually larger for deep events than shal-
low events. Also, the steeper the dip, the more accurate
the migration velocities need to be, since displacement
is proportional to dip.
Figure 4.0-13 shows a portion of a CMP-stacked
section after 3-D poststack time migration using a per-
cent range of stacking velocities. Note the subtle under-
and overmigration effects on dipping events below the
major unconformity represented by the strong, near-
horizontal reflection. Events A dipping up to the left
and B dipping up to the right cross over one another
on sections that correspond to 90 percent and 95 per-
cent of stacking velocities, indicating undermigration.
The same events are split away from one another in
opposite directions on sections that correspond to 105
percent and 110 percent of stacking velocities. The most
overall acceptable image is seen on the section that cor-
responds to the 100 percent of stacking velocities.
Accuracy in event positioning after migration actu-
ally depends on the combined effects of the performance
of the migration algorithm used and the velocity errors.
For example, the inherently undermigrating character
of a 45-degree finite-difference algorithm can be, for an
event with a specific dip, coincidentally counterbalanced
by the overmigration effect of erroneously too high ve-
locities. In the presence of large vertical velocity gradi-
ents, a two-pass 3-D migration can also cause overmi-
gration of steep dips (in the form of lateral translation)
even with the correct migration velocities.
Figure 4.0-14a shows a portion of a migrated
stacked section. Although this section does not contain
steep dips, accurate imaging of the faults along the low-
relief structures can be important to the interpreter.
Note the slight undermigration, which may be caused
by any of the following: (a) error in migration veloci-
ties, (b) a dip-limited algorithm that failed to focus the
diffraction energy adequately, or (c) a possible 3-D be-
havior of the geometry of the reflector. The section in
Figure 4.0-14a has been migratedwith a dip-limited al-
gorithm. Using a proper algorithm, with the same veloc-
ities, we get the migrated section in Figure 4.0-14b. The
resulting section shows slight overmigration, which can
be attributed to errors in migration velocities. Lowering
the velocities gives the improved, but not completely ac-
curate, image in Figure 4.0-14c. Perhaps, the remaining
issues in imaging may be attributed to 3-D effects.
As demonstrated by the example in Figure 4.0-14,
migration results generally are self-evident — under-
and overmigration often can be recognized on a mi-
grated section. Problems in imaging often are traced to
accuracy in migration velocities. I consider migration
velocities as the weak link between the seismic section
and the geologic cross-section.
In the next section, basic principles of migration
are presented and the Kirchhoff summation, finite-
difference, frequency-space, and frequency-wavenumber
algorithms are reviewed. Practical aspects of the migra-
tion algorithms are expounded in Sections 4.2 through
4.5. Specifically, key parameters for each category of
the migration algorithms are analyzed using appropri-
ate synthetic and field data examples. Further aspects of
476 Seismic Data Analysis
FIG. 4.0-13. A portion of a CMP-stacked section after 3-
D poststack time migration using, from top to bottom, 90,
95, 100, 105 and 110 percent of stacking velocities. Note the
subtle under- and overmigration effects on dipping events
below the major unconformity represented by the strong,
near-horizontal reflection. Note the effect of velocities used
in migration on the positioning of the event A dipping up
to the left and event B dipping up to the right.
migration in practice, including spatial aliasing, migra-
tion response to random noise, line length, and irregular
topography are discussed in Section 4.6. The problem
of conflicting dips with different stacking velocities that
requires dip-moveout (DMO) correction and prestack
time migration, and the accompanying topic on migra-
tion velocity analysis are deferred until Chapter 5. The
problem of imaging beneath complex structures that re-
quires earth imaging and modeling in depth is discussed
in Chapters 8 and 9, respectively.
FIG. 4.0-14. (a) A portion of a migrated CMP stack; note
the subtle undermigration at fault locations A and B caused
by the use of a dip-limited algorithm; (b) same data set but
migrated with an algorithm with no dip limitation; note
the subtle overmigration most likely due to erroneously too
high velocities; (c) same data set migrated with the same
algorithm as in (b) but with velocities adjusted to prevent
overmigration.
4.1 MIGRATION PRINCIPLES
Consider the dipping reflector CD of the simple ge-
ologic section in Figure 4.1-1a. We want to obtain a
zero-offset section along the profile Ox. As we move the
source-receiver pair (s, g) along Ox, the first normal-
incidence arrival from the dipping reflector is recorded
at location A. In this discussion, we assume a normal-
ized constant-velocity medium v = 1 so that time and
depth coordinates become interchangeable. The reflec-
tion arrival at location A is indicated by point C
�
on the
zero-offset time section in Figure 4.1-1b. As we move
from location A to the right, normal-incidence arrivals
are recorded from the dipping reflector CD. The last
arrival is recorded at location B, which is indicated by
point D
�
in Figure 4.1-1b. In this experiment, diffrac-
tions off the edges of reflector CD are excluded to sim-
plify the discussion.
Migration 477
FIG. 4.1-1. Migration principles: The reflection segment
C
�
D
�
in the time section as in (b), when migrated, is moved
updip, steepened, shortened, and mapped onto its true sub-
surface location CD as in (a). (Adapted from Chun and
Jacewitz, 1981.)
Compare the geologic section in Figure 4.1-1a,
which is in depth, with the zero-offset seismic section
in Figure 4.1-1b, which is in time. The true subsurface
position of reflector CD is superimposed onto the time
section for comparison. Clearly the true geologic posi-
tion of reflector CD is not the same as the reflection
event position C
�
D
�
.
From this simple geometric construction, note that
the reflection in the time section C
�
D
�
must be migrated
to its true subsurface position CD in the depth section.
The following observations can be made from the geo-
metric description of migration in Figure 4.1-1:
(a) The dip angle of the reflector in the geologic section
is greater than in the time section; thus, migration
steepens reflectors.
(b) The length of the reflector, as seen in the geologic
section, is shorter than in the time section; thus,
migration shortens reflectors.
(c) Migration moves reflectors in the updip direction.
The example in Figure 4.0-1 demonstrates the
above observations. In particular, the dipping event (B)
has moved in the updip direction, become shorter, and
steepened after migration (A).
As mentioned in the previous section, conventional
migration output is displayed in time, as is the input
stacked section. To distinguish the two time axes, we
will denote the time axis on the stacked section as t
— event time in the unmigrated position, and the time
axis on the migrated section as τ — event time in the
migrated position.
We shall now examine the horizontal and verti-
cal displacements as seen on the migrated time section.
From Figure 4.1-2, consider a reflector segment CD. As-
sume that CD migrates to C
�
D
�
and that point E
�
on
C
�
D
�
migrates to point E on CD. The horizontal and
vertical (time) displacements — d
x
and d
t
, and the dip
∆τ/∆x, all measured on the migrated time section (Fig-
ure 4.1-2), can be expressed in terms of medium velocity
v, traveltime t, and apparent dip ∆t/∆x as measured
on the unmigrated time section (Figure 4.1-2). Chun
and Jacewitz (1981) derived the following formulas:
d
x
=
v
2
t
4
∆t
∆x
, (4− 1)
d
t
= t
�
1−
�
1−
�
v∆t
2∆x
�
2
�
, (4− 2)
∆τ
∆x
=
∆t
∆x
1
�
1−
�
v∆t
2∆x
�
2
. (4− 3)
To gain a quantitative insight into these expres-
sions, we consider a numerical example. For a realis-
tic velocity function that increases with depth, consider
five reflecting segments at various depths. For simplic-
ity, assume that quantity ∆t/∆x is the same for all (10
ms per 25-m trace spacing). From equations (4-1), (4-
2), and (4-3), compute the horizontal and vertical dis-
placements d
x
and d
t
and the dips (in ms/trace) after
migration. The results are summarized in Table 4-2.
Refer to Table 4-2 and equations (4-1), (4-2), and
(4.3) and make the following observations:
(a) The time dip ∆τ/∆x on the migrated section is
always greater than the time dip ∆t/∆x on the
unmigrated section.
(b) The horizontal displacement d
x
increases with
event time t in the unmigrated position. At 4 s,
the horizontal displacement is more than 6 km.
478 Seismic Data Analysis
FIG. 4.1-2. Quantitative analysis of the migration process. Dipping event AB on the unmigrated section (left) is moved to
A
�
B
�
on the migrated section (right). The event after migration also is superimposed on the unmigrated section to compare the
position of the event before and after migration. Point C on dipping reflector AB is moved to C
�
after migration. The amount
of horizontal displacement d
x
, vertical displacement d
t
, and the dip ∆τ/∆x after migration is calculated from equations (4-1),
(4-2) and (4-3).
Table 4-2. Horizontal andvertical displacements as a
result of migration of a series of dipping reflections with
the same apparent dip (10 ms/trace) as measured on
the unmigrated stacked section at various depths, and
changes in dip angle as measured on the migrated sec-
tion in time.
t v d
x
d
t
∆t/∆x ∆τ/∆x
(s) (m/s) (m) (s) (ms/trace) (ms/trace)
1 2500 625 0.134 10 11.5
2 3000 1800 0.400 10 12.5
3 3500 3675 0.858 10 14.0
4 4000 6400 1.600 10 16.7
5 4500 10125 2.820 10 23.0
(c) The horizontal displacement d
x
is a function of the
velocity squared. If there is a 20 percent error in
the velocity used in migration, then the event is
misplaced by an error of 44 percent.
(d) The vertical displacement d
t
also increases with
time and velocity.
(e) The steeper the event dip, the more the horizontal
and vertical displacements after migration.
In Figure 4.1-1a, assume that the zero-offset sec-
tion was recorded only between surface locations A and
B. The time section would include the event C
�
D
�
, but
when migrated, the event would migrate out of the sec-
tion, resulting in a blank migrated section (Figure 4.1-
1b). Therefore, the data on a stacked section are not
necessarily confined to the subsurface below the seismic
line. The converse is even more important; the structure
below the seismic line may not be recorded on the seis-
mic section. Suppose that the data were recorded only
between surface locations O and A. This time, the re-
sulting time section would be blank. So, we should not
record between A and B, and neither should we record
between O and A. Instead, we should record between O
and B in order to record the reflector of interest prop-
erly and also to migrate it properly.
In areas with a structural dip, line length must be
chosen by considering the horizontal displacements of
dipping events from the structures causing the events.
This is an important consideration, especially in 3-D
seismic work. The areal surface coverage of a survey
usually is larger than the areal subsurface coverage of
interest.
To achieve a complete image of a dipping reflec-
tor, also, the recording time must be long enough. For
Migration 479
FIG. 4.1-3. (a) A zero-offset section modeled from the dipping reflectors shown in (b). The medium velocity is constant 3500
m/s and the trace spacing is 25 m. The true dip angles of the reflectors vary from 0 to 45 degrees at 5-degree increment.
Migration of the zero-offset section (a) yields the model of the dipping reflectors (b).
FIG. 4.1-4. A portion of a CMP stack (a) before and (b) after migration. Note the group of events with a range of dips that
fan out from a fault plane. Migration has moved them in the updip direction, made them shorter and steeper.
480 Seismic Data Analysis
FIG. 4.1-5. Curved reflecting interfaces (synclines and an-
ticlines) (a) before and (b) after migration. See text for de-
tails. (Modeling courtesy Union Oil Company.)
FIG. 4.1-6. (a) A velocity-depth model consisting of a syn-
clinal reflector; (b) selected normal-incidence arrivals on the
zero-offset section. Trace the bowtie in the time section.
example, if only OE seconds were recorded (Figure 4.1-
1), then the recorded segment C
�
D
��
would yield only
part of the complete image CD. An excellent example of
recording deeper in time and with longer line length for
steeper dips is shown in Figure 4.0-1. Proper imaging of
the salt dome boundary required that data be recorded
for more than 6 s.
The migration concepts described above are
demonstrated further by the dipping events model in
Figure 4.1-3. The edge diffractions are included here.
The dipping reflectors on the zero-offset section are
steepened, shortened, and moved in the updip direc-
tion as a result of migration. A field data example of a
series of dipping events on a stacked section before and
after migration is shown in Figure 4.1-4. Note that the
steeper the dip, the more the event moves after migra-
tion.
So far, only linear reflectors were considered. We
now consider a more realistic geologic situation that
involves curved reflecting interfaces. Figure 4.1-5 shows
Migration 481
FIG. 4.1-7. A portion of a CMP stack (a) before and (b) after migration. Anticlines appear bigger, while synclines appear
smaller than their actual sizes on the unmigrated section (a).
three synclines and a small anticlinal feature. The syn-
clines appear as bowties on the zero-offset section. By
using the principles deduced from the geometry of Fig-
ure 4.1-1, note that as a result of migration, segment
A of the bow tie moves in the updip direction to the
left. Similarly, segment B moves to the right, while flat-
topped segment C does not move much at all. Conse-
quently, after migration the flanks of bow ties associated
with synclines are opened up. On the other hand, the
small anticline seems to be broader on the zero-offset
section than it is on the migrated section. Again note
that segmentD moves updip to the right, while segment
E moves updip to the left as a result of migration. Thus,
synclines broaden and anticlines compress as a result of
migration. Migration velocities also affect the apparent
size of the structure; higher velocities mean more mi-
gration and, hence, smaller anticlinal structure.
Why does a syncline look like a bowtie on the
stacked section? The answer is in Figure 4.1-6, where
a symmetric syncline is modeled. Given the subsurface
model in Figure 4.1-6a, the normal-incidence rays can
be computed to derive the zero-offset traveltime section
in Figure 4.1-6b. Only five CMP locations are shown
for clarity. At locations 2 and 4, there are two distinct
arrivals, while at location 3, there are three distinct ar-
rivals. By filling in the intermediate raypaths, the bow
tie character of the syncline can be constructed on the
time section. Complete the procedure by tracing the
traveltime trajectory in Figure 4.1-6b.
Two field data examples containing synclinal and
anticlinal structures are shown in Figures 4.1-7 and
4.1-8. In Figure 4.1-7, note that the synclinal feature
broadens and the anticlinal feature narrows as a result
of migration. In Figure 4.1-8, the bow ties associated
with two small synclinal basins A and B grow larger in
depth. After migration, the bowties are untied and the
synclines are delineated.
Kirchhoff Migration
Claerbout (1985) uses the harbor example in Figure 4.1-
9 to describe the physical principles of migration. As-
sume that a storm barrier exists at some distance z
3
from the beach and that there is a gap in the barrier.
Imagine a calm afternoon breeze that comes from the
ocean as a plane incident wave. The wavefront is par-
allel to the storm barrier. As we walk along the beach
line, we see a wavefront different from a plane wave.
The gap on the storm barrier has acted as a secondary
source and generated the semicircular wavefront that is
propagating toward the beach.
If we did not know about the storm barrier and the
gap, we might want to lay out a receiver cable along the
beach to record the approaching waves. This experiment
482 Seismic Data Analysis
FIG. 4.1-8. A portion of a CMP stack (a) before and (b) after migration. Migration unties the bowties and turns them into
synclines below A and B.
is illustrated in Figure 4.1-10 with the recorded time
section. Physicists call the gap on the barrier a point
aperture. It is somewhat similar to a point source, since
both generate circular wavefronts. However, the ampli-
tudes on the wavefront that propagate outward from
a point source are isotropic, while those from a point
aperture are angle-dependent. The point apertureon
the barrier acts as a Huygens’ secondary source.
From the beach experiment, we find that Huygens’
secondary source responds to a plane incident wave and
generates a semicircular wavefront in the x − z plane.
The response in the x− t plane is the diffraction hyper-
bola shown in Figure 4.1-11.
Imagine that the subsurface consists of points along
each reflecting horizon that behave much as the gap
on the storm barrier. From Figure 4.1-12, these points
act as Huygens’ secondary sources and produce hyper-
bolic traveltime trajectories. Moreover, as the sources
(the points on the reflecting interface) get closer to each
other, superposition of the hyperbolas produces the re-
sponse of the actual reflecting interface (Figure 4.1-13).
In terms of the harbor example, this is like assuming
that the barrier is wiped out by a storm so that the
primary incident plane wave reaches the beach with-
out modification. The diffraction hyperbolas, which are
caused by sharp discontinuities at both ends of the re-
flector in Figure 4.1-13, remain. These hyperbolas are
equivalent to diffractions seen at fault boundaries on
stacked sections.
In summary, we find that reflectors in the subsur-
face can be visualized as being made up of many points
that act as Huygens’ secondary sources. We also find
that the zero-offset section consists of a superposition
of the many hyperbolic traveltime responses. Moreover,
when there are discontinuities (faults) along the reflec-
tor, diffraction hyperbolas often stand out.
Migration 483
FIG. 4.1-9. The gap in the barrier acts as Huygens’ sec-
ondary source, causing the circular wavefronts that approach
the beach line. (Adapted from Claerbout, 1985.)
FIG. 4.1-10.Waves recorded along the beach generated by
Huygens’ secondary source (the gap in the barrier in Figure
4.1-9) have a hyperbolic traveltime trajectory.
FIG. 4.1-11. A point that represents a Huygens’ secondary
source (a) produces a diffraction hyperbola on the zero-offset
time section (b). The vertical axis in this section is two-way
time, while the vertical axis in the time section in Figure
4.1-10 is one-way time.
FIG. 4.1-12. Superposition of the zero-offset responses (b)
of a discrete number of Huygens’ secondary sources as in
(a).
FIG. 4.1-13. Superposition of the zero-offset responses (b)
of a continuum of Huygens’ secondary sources as in (a).
484 Seismic Data Analysis
Diffraction Summation
Huygens’ secondary source signature is a semicircle in
the x − z plane and a hyperbola in the x − t plane.
This characterization of point sources in the subsurface
leads to two practical migration schemes. Figure 4.1-
14a shows a zero-offset section that consists of a single
arrival at a single trace. This event migrates to a semi-
circle (Figure 4.1-14b). From Figure 4.1-14, note that
the zero-offset section recorded over a constant-velocity
earth model consisting of a semicircular reflecting inter-
face contains a single blip of energy at a single trace as
in Figure 4.1-14a. Since this recorded section consists of
an impulse, the migrated section in Figure 4.1-14b can
be called the migration impulse response. An alternate
scheme for migration results from the observation that
a zero-offset section consisting of a single diffraction hy-
perbola migrates to a single point (Figure 4.1-15b).
The first method of migration is based on the su-
perposition of semicircles, while the second method is
based on the summation of amplitudes along hyperbolic
paths. The first method was used before the age of dig-
ital computers. The second method, which is known as
the diffraction summation method, was the first com-
puter implementation of migration.
The migration scheme based on the semicircle su-
perposition consists of mapping the amplitude at a sam-
ple in the input x− t plane of the unmigrated time sec-
tion onto a semicircle in the output x − z plane. The
migrated section is formed as a result of the superposi-
tion of the many semicircles.
The migration scheme based on diffraction sum-
mation consists of searching the input data in the x− t
plane for energy that would have resulted if a diffract-
ing source (Huygens’ secondary source) were located at
a particular point in the output x−z plane. This search
is carried out by summing the amplitudes in the x − t
plane along the diffraction curve that corresponds to
Huygens’ secondary source at each point in the x − z
plane. The result of this summation then is mapped
onto the corresponding point in the x − z plane. As
noted early in this section, within the context of time
migration, however, the summation result actually is
mapped onto the x− τ plane, where τ is the event time
in the migrated position.
The curvature of the hyperbolic trajectory for am-
plitude summation is governed by the velocity function.
The equation for this trajectory can be derived from the
geometry of Figure 4.1-15. A formal derivation also is
provided in Section D.2. Assuming a horizontally lay-
ered velocity-depth model, the velocity function used
to compute the traveltime trajectory is the rms veloc-
ity at the apex of the hyperbola at time τ (Section 3.1).
FIG. 4.1-14. Principles of migration based on semicircle
superposition. (a) Zero-offset section (trace interval, 25 m;
constant velocity, 2500 m/s), (b) migration.
FIG. 4.1-15. Principles of migration based on diffraction
summation. (a) Zero-offset section (trace interval, 25 m; con-
stant velocity, 2500 m/s), (b) migration. The amplitude at
input trace location B along the flank of the traveltime hy-
perbola is mapped onto output trace location A at the apex
of the hyperbola by equation (4-4).
Migration 485
From the triangle COA in Figure 4.1-15a, we note that
t
2
= τ
2
+
4x
2
v
2
rms
. (4− 4)
Having computed the input time t, the amplitude
at input location B is placed on the output section at
location A, corresponding to the output time τ at the
apex of the hyperbola.
From Section 3.1, reflection traveltimes in a layered
earth approximate small-spread hyperbolas. This may
seem to impose a serious restriction on the aperture
width — the lateral extent of the diffraction hyperbola,
in the summation process. However, the small-spread
approximation is valid even at large distances from the
apex, and the errors associated with it are insignificant
at late times. In practice, this approximation is not usu-
ally an issue.
Amplitude and Phase Factors
Now consider several factors associated with the am-
plitude and phase behavior of the waveform along the
diffraction hyperbola. From Figure 4.1-9, given the al-
ternative of standing at location A or B, we intuitively
think that it is safer to stand at location B. This is
because the wave amplitude at location A, which is on
the z-axis, is stronger than the wave amplitude at loca-
tion B, which is at an oblique angle from the z-axis. As
mentioned earlier, this is one difference between a point
source with uniform amplitude response at all angles
and the point aperture that produces a wavefront with
angle-dependent amplitudes. This angle dependence of
amplitudes, which is described by the obliquity factor,
should be considered before summation. To correct for
the obliquity factor, the amplitude at location B in Fig-
ure 4.1-15a is scaled by the cosine of the angle between
BC and CA before it is placed at output location A.
Another factor is the spherical divergence of wave
amplitudes. Again, from Figure 4.1-9, given the alterna-
tive of standing at location B or C, we prefer to stand at
location C. The reason for this is that the wave ampli-
tude along the wavefront at location C, whichis farther
from the point aperture source, is weaker than the wave
amplitude at location B. Wave energy decays as (1/r
2
),
where r is the distance from the source to the wavefront,
while amplitudes decay as (1/r). Thus, amplitudes must
be scaled by factor (1/r) before summation for wave
propagation in three dimensions.
Finally, there is a third factor that involves the
inherent property of Huygens’ secondary source wave-
form. This factor is difficult to explain from a physical
viewpoint. Nevertheless, it is obvious from Figure 4.1-
13 that Huygens’ secondary sources must respond as a
wavelet along the hyperbolic paths with a unique phase
and frequency characteristic. Otherwise, there would be
no amplitude cancelation when they are close to one an-
other. The waveform that results from the summation
must be restored in both phase and amplitude.
In summary, we must consider the following three
factors before diffraction summation:
(a) The obliquity factor or the directivity factor, which
describes the angle dependence of amplitudes and
is given by the cosine of the angle between the di-
rection of propagation and the vertical axis z (Fig-
ure 4.1-15).
(b) The spherical spreading factor, which is propor-
tional to
�
1/vr for 2-D wave propagation, and
(1/vr) for 3-D wave propagation.
(c) The wavelet shaping factor, which is designed with
a 45-degree constant phase spectrum and an am-
plitude spectrum proportional to the square root of
the frequency for 2-D migration. For 3-D migration,
the phase shift is 90 degrees and the amplitude is
proportional to frequency.
Kirchhoff Summation
The diffraction summation that incorporates the obliq-
uity, spherical spreading and wavelet shaping factors
is called the Kirchhoff summation, and the migration
method based on this summation is called the Kirchhoff
migration. To perform this method, multiply the input
data by the obliquity and spherical spreading factors.
Then apply the filter with the above specifications and
sum along the hyperbolic path that is defined by equa-
tion (4-4). Place the result on the migrated section at
time τ corresponding to the apex of the hyperbola. In
practice, the order of the filter application, specified by
factor (c), and summation can be interchanged without
sacrificing accuracy because the summation is a linear
process and the filter is independent of time and space.
The velocity used in equation (4-4) is taken as the
rms velocity, which can be allowed to vary laterally.
However, lateral variation in velocity distorts the hy-
perbolic nature of the diffraction pattern and somehow
must be considered. The value for the rms velocity typi-
cally is that of the output time sample; that is, the apex
time τ of the hyperbola.
What was determined from a physical point of view
in the preceding discussion can be rigorously described
by the integral solution to the scalar wave equation.
Schneider (1978), Berryhill (1979) and Berkhout (1980)
are excellent references for the mathematical treatment
of the Kirchhoff migration method. The integral solu-
tion of the scalar wave equation yields three terms; the
far-field term which is proportional to (1/r), and two
other terms which are proportional to (1/r
2
). Hence,
486 Seismic Data Analysis
it is the far-field term that makes the most contribu-
tion to the summation that is used in practical imple-
mentation of Kirchhoff migration. The output image
P
out
(x
0
, z = vτ/2, t = 0) at a subsurface location (x
0
, z)
using only the far-field term is computed from the 2-D
zero-offset wavefield P
in
(x, z = 0, t), which is measured
at the surface (z = 0), by the following summation over
a spatial aperture
P
out
=
∆x
2π
�
x
�
cos θ
√
v
rms
r
ρ(t) ∗ P
in
�
, (4− 5)
where v
rms
is the rms velocity at the output point
(x
0
, z) and r =
�
(x− x
0
)
2
+ z
2
, which is the distance
between the input (x, z = 0) and the output (x
0
, z)
points. The asterisk denotes convolution of the rho fil-
ter ρ(t) with the input wavefield P
in
.
The rho filter ρ(t) corresponds to the time deriva-
tive of the measured wavefield, which yields the 90-
degree phase shift and adjustment of the amplitude
spectrum by the ramp function ω of frequency (Ta-
ble A-1 of Appendix A). For 2-D migration, the half-
derivative of the wavefield is used. This is equivalent
to the 45-degree phase shift and adjustment of the am-
plitude spectrum by a function of frequency defined as
√
ω. Since the rho filter is independent of the spatial
variables, it actually can be applied to the output of the
summation in equation (4-5). Finally, the far-field term
in equation (4-5) is proportional to the cosine of the an-
gle of propagation (the directivity term or the obliquity
factor) and is inversely proportional to vr (the spherical
spreading term) in three dimensions. In two dimensions,
the spherical spreading term is
√
vr.
Equation (4-5) can be used to compute the wave-
field at any depth z. The ouput image P
out
is computed
at (x
0
, z = vτ/2, t = 0) using the input wavefield P
in
at (x, z = 0, t− r/v). To obtain the migrated section at
an output time τ , equation (4-5) must be evaluated at
z = vτ/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 com-
plete migrated section is obtained by performing the
summation in equation (4-5) and setting t = 0 for each
output location. The range of the summation is called
the migration aperture.
Finite-Difference Migration
To describe the physical basis of finite-difference migra-
tion, recall the harbor example of Figure 4.1-9. Instead
of taking the section recorded along the beach, which
contains the diffraction hyperbola, then collapsing it to
get the migrated section in Figure 4.1-15, consider the
following alternative procedure. Again, start with the
wavefield recorded along the beach (Figure 4.1-16a). As-
sume that the barrier is 1250 m from the beach. Now
move the recording cable into the water, 250 m from the
beach. Start recording at the instant the plane wave hits
the barrier. The recorded section is shown in Figure 4.1-
16b. Move the cable 500 m from the beach and record
the section in Figure 4.1-16c, followed by a recording
750 m from the beach to obtain the section in Figure
4.1-16d. Finally, 1000 m from the beach, record the sec-
tion shown in Figure 4.1-16e.
Note that each recording yields a hyperbola in
which the apex moves closer to zero time. The actual
extent of the recording cable is denoted by the solid line
on top of each frame. Had we recorded at the barrier
(1250 m from the beach), the apex of the hyperbola
would be positioned at t = 0.
In Kirchhoff migration, the diffraction hyperbola
is collapsed by summing the amplitudes, then placing
them at the apex. The alternative approach implied
by the result of the experiment shown in Figure 4.1-
16 is to use the hyperbola recorded a distance away
from the beach to construct the hyperbola that would
be recorded at another distance closer to the source of
the diffraction hyperbola. The process is stopped when
the hyperbola collapses to its apex. In the harbor ex-
periment, this collapse occurs when the receiver cable
coincides with the barrier, or, equivalently, when t = 0.
As stated in the introductory section, this is called the
imaging principle.
Downward Continuation
The harbor experiment described above can be simu-
lated in the computer. Pretend that moving the receiver
cable from thebeach into the water closer to the barrier
is like moving the receiver cable from the surface down
into the earth closer to the reflectors. Think of the gap
on the barrier as equivalent to a point diffractor on a re-
flecting interface causing the diffraction hyperbola (Fig-
ure 4.1-17a). Start with the wavefield recorded at the
surface and move the receivers down to depth levels at
finite intervals. Downward continuation of the upcom-
ing wavefield at the surface, therefore, can be considered
equivalent to lowering the receivers into the earth.
The computer-simulated wavefields at these differ-
ent depths are shown in Figure 4.1-17. By applying the
imaging principle at each depth, the entire wavefield is
imaged. The final output from this process is the mi-
grated section. The last section (panel f) at 1250 m
has only one arrival at t = 0. The recording cable is
on the storm barrier and the arrival from the gap oc-
curs at t = 0. As the cable moved into the ocean and
Migration 487
FIG. 4.1-16. Moving the receiver cable in the harbor experiment (Figure 4.1-9) from the beach into the water at discrete
intervals parallel to the beach line. Numbers on top indicate the distance of the receiver cable from the beach line.
FIG. 4.1-17. Computer simulation of the experiment illustrated in Figure 4.1-16. Here, we downward continue the receivers
at discrete depth intervals. The numbers on top indicate the distance of the receiver cable from the surface, z = 0.
recorded closer to the barrier, the recorded diffraction
hyperbola arrived earlier, and became shorter and more
compressed. It collapsed to a point when the receivers
coincided with the storm barrier over which the source
point forms a gap.
There is one important difference between the
physical experiment in Figure 4.1-16 and the computer-
simulated downward-continuation experiment in Figure
4.1-17. The receiver cable is the same at each step in
Figure 4.1-16, whereas the effective cable length gets
shorter and shorter toward the source (the gap in the
barrier) in Figure 4.1-17. This is because we started by
recording the wavefield at the surface (Figure 4.1-16a)
with a finite cable length. The recorded information is
confined to within the two raypaths depicted on the
section in Figure 4.1-17a. As the cable moves closer to
the source, the effective receiver cable containing the
information is confined to smaller and smaller lengths.
Although receivers are lowered vertically, energy moves
down along raypaths it originally took on the way up.
To relate these recordings at different depths (Fig-
ure 4.1-17), we superimpose them as shown in Figure
4.1-18a. Moreover, the recordings can be shifted so that
the apexes of the hyperbolas coincide and are positioned
at a time that is equivalent to the distance from the sur-
face to the diffractor as shown in Figure 4.1-18b. This
is called time retardation.
Reconsider the results from the computer simula-
tion of the harbor experiment in Figure 4.1-17. Sup-
pose we stopped recording at a depth of 1000 m before
FIG. 4.1-18. (a) Superposition of the time sections in Fig-
ure 4.1-17; (b) removing the translational effect by retar-
dation to place the energy at the apex of the hyperbola
obtained initially along the beach line.
488 Seismic Data Analysis
the barrier. The original hyperbola in Figure 4.1-17a
was partially collapsed at this depth (Figure 4.1-17e).
Therefore, downward continuing to a depth short of the
true depth of the source causes undermigration. Diffrac-
tions and dipping events also are undermigrated if in-
correctly low velocities are used for migration.
Assume that the recording continues and passes
beyond barrier position z
3
(Figure 4.1-9). We infer that
the focused energy on the section at this depth (Figure
4.1-17f) would propagate through the focal point and
turn into hyperbolas that are the mirror images of those
in Figures 4.1-17a through 4.1-17e. We have downward
continued more than necessary. This yields overmigra-
tion, which also is caused by incorrectly high velocities.
From these observations, note that downward continu-
ing to a wrong depth is like downward continuing with
the wrong velocity (Doherty and Claerbout, 1974).
Another important issue to consider is how often
the extrapolated wavefield should be computed. When
going from one frame to another in Figure 4.1-17, what
should the depth step size be? This is discussed in detail
later in Section 4.3.
Differencing Schemes
Finite-difference migration algorithms are based on dif-
ferential solutions to the scalar wave equation that are
used to downward continue the input wavefield recorded
at the surface. A simple numerical example illustrates
the finite-difference method of solving differential equa-
tions (Claerbout, 1985). Assume that you have $100
today. Given an annual inflation rate of 10 percent, for
the same buying power next year, you need $110. A
computer algorithm can determine the face value of the
present $100 in future years. Table 4-3 shows the results
of extrapolation from one year to the next.
Given the present value, 100, find future values in
the data column. The following equation solves for the
unknown x:
(1.0)× x+ (−1.1)× (100) = 0, (4− 6a)
which yields x = 110.We used a two-point operator and
aligned it with the data column as indicated in Table
4-3. Similarly, we have
(1.0)× x+ (−1.1)× (110) = 0, (4− 6b)
which yields x = 121. By using the new value for x, we
obtain
(1.0)× x+ (−1.1)× (121) = 0, (4− 6c)
which yields x = 133, and so on. By moving the oper-
ator down in the time direction as shown in Table 4-3,
we extrapolate the data column into the future.
Table 4-3. Simple extrapolation of a data vector in
time.
Operator Data Time Step
-1.1 100 0
1.0 x 1
100 0
-1.1 110 1
1.0 x 2
100 0
110 1
-1.1 121 2
1.0 x 3
Equation (4-6a) is generalized as
(1.0)× P (t+ 1) + (−1.1)× P (t) = 0, (4− 7a)
which is rewritten in the form
P (t+ 1)− P (t) = (0.1)× P (t), (4− 7b)
where t is the time variable and P is the quantity being
extrapolated. Instead of defining the time interval as
one unit, we can define it as an arbitrary increment
of time ∆t. Also, assume that the inflation rate is a.
Equation (4-7b) then takes the more general form
P (t+∆t)− P (t) = aP (t). (4− 8a)
Alternatively, we could use the average of the
present and future values on the right side of this equa-
tion:
P (t+∆t)− P (t) =
a
2
�
P (t+∆t) + P (t)
�
. (4− 8b)
Equations (4-8a) and (4-8b) now can be put into
the form of equation (4-6a) as
P (t+∆t) +
−1− a
�
P (t) = 0, (4− 9a)
and
�
1−
a
2
P (t+∆t) +
�
−1−
a
2
P (t) = 0. (4− 9b)
By using either equation (4-9a) or equation (4-9b), we
compute the future values of P (t) from a given initial
value as shown in Table 4-4.
The operator in which the coefficient of the future
value P (t +∆t) is unity is called the explicit operator.
Stability of the finite-difference solution — the problem
Migration 489
Table 4-4. Application of two-point implicit and ex-
plicit operators to extrapolate data P from t to t+∆t.
Explicit Implicit Data
Operator Operator Column
−1− a −1− a/2 P (t)
1 1− a/2 P (t+∆t)
of wave amplitudes growing from one extrapolation step
to another, can be an issue with this type of operator
(Section D.6) An implicit operator produces stable re-
sults because of averaging on the right side of equation
(4-9b), known as the Crank-Nicolson scheme. For the
differential equations used in finite-difference migration
algorithms, such as the parabolic equation described in
Section D.3, scalar a becomes a matrixcoefficient. Im-
plicit schemes require inversion of this matrix. However,
no inversion is needed with explicit schemes, since fu-
ture values can be written explicitly in terms of only
past values.
Equation (4-9a) is rewritten by redefining scalar a
as a∆t to obtain
P (t+∆t)− P (t)
∆t
= aP (t). (4− 10)
The left side of equation (4-10) is the discrete repre-
sentation of the continuous derivative of P with respect
to time, dP/dt. Therefore, equation (4-10) is the finite-
difference equation that corresponds to the differential
equation
dP
dt
= aP (t). (4− 11)
We have derived the differential equation that de-
scribes the inflation of money (equation 4-11). Now con-
sider the analysis in reverse order. We start with the dif-
ferential equation (4-11), and write the corresponding
difference equation (4-10), which is the equation that
is solved in the computer. This equation is written in
either the explicit (equation 4-9a) or implicit (equation
4-9b) form to extrapolate the present value of P to the
future.
This example illustrates how finite-difference
schemes can solve differential equations in the com-
puter. The scalar wave equation can be treated in a
similar, but more complicated manner. Complications
arise because it is a partial differential equation that
contains the second derivatives of the wavefield with re-
spect to depth, time, and spatial axes. Setting up the
computer algorithm is more involved and is not dis-
cussed here. Claerbout (1976, 1985) provides details of
various aspects of the finite-difference migration meth-
ods.
Rational Approximations for Implicit Schemes
The scalar wave equation is a two-way wave equation in
depth that describes propogation of both upcoming and
downgoing waves. If we consider the resulting wavefield
from the exploding reflectors model as the upcoming
waves, then we are really interested in a one-way wave
equation to downward continue the upcoming waves. In
fact, we normally use some rational approximation to
the one-way wave equation in finite-difference imple-
mentations.
To get the actual differential equation to be used
in downward extrapolation of the upcoming waves, and
therefore to perform a finite-difference migration, the
general strategy is as follows:
(a) Start with the two-way scalar wave equation:
∂
2
P
∂x
2
+
∂
2
P
∂z
2
−
1
v
2
(x, z)
∂
2
P
∂t
2
= 0, (4− 12)
where x and z are the space variables, t is the time
variable, v is the velocity of wave propagation, and
P (x, z, t) is the pressure wavefield.
(b) Assume constant velocity and perform 3-D Fourier
transform of the pressure wavefield. This is equiv-
alent to substituting the plane-wave solution
exp(ik
x
x+ ik
z
z− iωt) to equation (4-12). The sub-
stitution yields the dispersion relation between the
transform variables
k
z
= ∓
�
ω
2
v
2
− k
2
x
, (4− 13a)
where k
x
and k
z
are the wavenumbers in the x
and z directions, and ω is the angular temporal
frequency.
(c) We are interested in upcoming waves, hence we
only need one of the two solutions. We also want to
invoke the exploding reflector model by replacing v
by v/2 to obtain the following paraxial dispersion
relation
k
z
=
2ω
v
�
1−
�
vk
x
2ω
2
, (4− 13b)
where the horizontal wavenumber k
x
has been nor-
malized with respect to 2ω/v.
(d) Make a rational approximation to the square-root
expression in equation (4-13b) so as to derive a
differential equation (Sections D.3 and D.4). This
490 Seismic Data Analysis
approximation imposes a dip limitation to the dif-
ferential equation. One approximation to the dis-
persion relation given by equation (4-13b) is ob-
tained by Taylor expansion of the square root and
retaining the first two terms in the series (Section
D.3)
k
z
=
2ω
v
�
1−
1
2
�
vk
x
2ω
2
�
. (4− 14a)
This dispersion equation is known as the 15-degree
approximation and is the basis for the first finite-
difference time migration algorithm developed by
Claerbout and Doherty (1972). Albeit no longer in
use, we shall review the 15-degree finite-difference
algorithm for its historical significance.
(e) Operate on the pressure wavefield P (k
x
, k
z
, ω) with
the approximate form of the dispersion relation
given by equation (4-14a), and inverse Fourier
transform in the z direction to get the differential
form of the approximate one-way wave equation.
∂
∂z
P (k
x
, z, ω) = −i
2ω
v
�
1−
1
2
�
vk
x
2ω
2
�
P (k
x
, z, ω).
(4− 14b)
(f) Recall from Figure 4.1-17 that, after each
downward-continuation step, we retard the wave-
field by translating it in time so that after migra-
tion, events appear in their correct depth locations.
The time retardation is done by applying a linear
phase shift to the pressure wavefield P
Q = P exp(−iωτ), (4− 15a)
where the retarded time is
τ =
	
z
0
dz
v¯(z)
, (4− 15b)
andQ is the retarded wavefield. The velocity v¯(z) is
the horizontally averaged v(x, z). Substitute equa-
tion (4-15a) into (4-14b) to obtain the differen-
tial equation associated with the 15-degree finite-
difference algorithm in two parts
∂
2
Q
∂z∂t
=
v
4
∂
2
Q
∂x
2
, (4− 16a)
and
∂Q
∂z
= 2
�
1
v¯(z)
−
1
v(x, z)
�
∂Q
∂t
, (4− 16b)
where Q is the retarded wavefield. Derivation
of equations (4-16a,b) is based on the assump-
tion that velocity varies vertically. Nevertheless,
in practice, the velocity function in equations (4-
16a,b) can be varied laterally, provided the varia-
tion is smooth. Equation (4-16a) accounts for col-
lapsing diffraction energy to the apex of the trav-
eltime curve only. Hence, it is referred to as the
diffraction term. When lateral velocity variations
are significant, the diffraction curve is somewhat
like a skewed hyperbola with its apex shifted later-
ally away from the diffraction source. This lateral
shift is accounted for by the thin-lens term given
by equation (4-16b) (Section D.3). If the lateral
velocity variations are significant, then the thin-
lens term is not negligible. Migration algorithms
that implement both the diffraction and thin-lens
terms represented by equations (4-16a,b) generally
are two-step schemes that alternately solve these
two terms. To propagate one depth step, first apply
the diffraction term on wavefield Q. The thin-lens
term then is applied to the output from the diffrac-
tion calculation. A migration method that includes
the effects of the thin-lens term is called depth mi-
gration, since the output section is in depth. Depth
migration is warranted if there are strong lateral
variations in velocity; in this case, the coefficient
of the thin-lens term cannot be negligible. If we
assume that velocity varies only in the vertical di-
rection, then v¯(z) = v(x, z). This makes the thin-
lens term of equation (4-16b) vanish, and we are
left with the diffraction term of equation (4-16a).
A migration method that implements the diffrac-
tion term (equation 4-16a), only, is known as time
migration, the output of which is in time τ of equa-
tion (4-15b). When recast in terms of the τ vari-
able, equation (4-16a) takes the form
∂
2
Q
∂τ∂t
=
v
2
8
∂
2
Q
∂x
2
. (4− 17)
This is the parabolic equation for time migration.
(g) Finally, write down the difference forms of the dif-
ferential operators either in implicit form to be used
in finite-difference solution of the parabolic equa-