Migration of Seismic Data

IENO GAZDAG AND PIERO SGUAZZERO

Invited Paper

Reflection seismology seeks to determine the structure of the more accurate than finite-difference methods in the
earth from seismic records obtained at the surface. The processing space-time coordinate frame. At the same time, the diffrac-
of these data by digital computers is aimed at rendering them more
tion summation migration was improved and modified on
comprehensiblegeologically. Seismic migration is one of these
processes. Its purpose is to "migrate" the recorded events to their the basis of theKirchhoff integral representation ofthe
correct spatial positions by backward projection or depropagation solutionofthe wave equation. The resultingmigration
based o n wave theoretical considerations. During the last 15 years procedure, known as Kirchhoff migration, compares favor-
several methods have appeared on the scene. The purpose of this ably with other methods. A relatively recent advance is
paper is to provide an overview of the major advances in this field.
Migration methods examined here fall in three major categories: I )
represented by reverse-time migration, which is related to
integralsolutions, 2) depth extrapolation methods, and 3) time wavefront migration.
extrapolation methods. Withinthese categories, the pertinent equa The aim of this paper is to present the fundamental
tions and numerical techniques are discussed in some detail. The concepts of migration. While it represents a review of the
topic of migration before stacking is treated separately with an state of the art, it is intended to be of tutorial nature. The
outline of two different approaches to this important problem.
reader is assumed to have no previous knowledge of seismic
processing.However, some familiarity
with Fourier
I. INTRODUCTION
transforms and their applications is assumed. The organiza-
The purpose of migration is to reconstruct the reflectivity tion of the paper is as follows. First, the basic concepts of
map of the earth from the seismicdata recorded at the seismic data representation are introduced. Next, migration
surface.Theseismicsignal recorded by a receiver (geo- schemes based on integral or summation methods are dis-
phone) is a superposition of seismic waves originating from cussed. Section IV deals with the derivation of the one-way
all possible directions in the subterrain. Thus the recorded wave equations. This is followed by the numerical aspects
events most often are not from reflectors directly below the of migration methods. Migration in object space, known as
receiver but from geological formations far away from the reverse-time migration, is presented in Section VI, and Sec-
point of recording. The term migration refers to the move- tion VI1 is devoted to migration before stacking.
ment of the observed events to their true spatial positions. The effect of migration on a seismic section depends on
Migration is an inverse process, inwhich the recorded the velocity and the reflectivity structure of the subterrain.
waves are propagated back to the corresponding reflector Flat, horizontal reflectors in a medium that has no signifi-
locations. The concept of migration can be summarized in cant velocity variations may not need to be migrated. O n
the following terms. In the process of seismic data acquisi- the other hand, seismic sections from media having steeply
tion the upward traveling waves are recorded at the surface. dipping reflectors and strong lateral velocity variations need
Inthemigration processthese recorded wavesare used to be migrated for correct interpretation. In the latter case,
either as initial conditions or boundary conditions for a migration can have a profound impact on the image, as we
wavefield governed by the wave equation. As a result, these shall see in the example ofmodel 2 (Figs. 11-14) to be
waves are propagated backward and in reverse time, from discussed in Section V.
the surface to the reflector locations.
Until the 196Os, migration was achieved by graphical 11. OF SEISMIC DATA
REPRESENTATION
methods. This was followed by diffraction summation and
wavefront migration based on ray theoretical considera- A.Data Acquisition
tions. In the 197Os, several important developments took
place. Based on thepioneeringworkof Jon Claerbout, Reflection seismology is an echo-ranging technique. An
migration methods based on wave theory were developed. acoustic source (shot) emits a short pulse and a set of
Efficientalgorithms developed fromsimplified finite-dif- recorders (geophones) register the reflected wavesat the
surface. The time series (sampled)dataassociated with a
ference approximations of the wave equations for down-
single shot and receiver is known as a trace. In typical
ward extrapolation. Another important processing advance
marineexploration (Fig. I), a boat tows a source and a
was made by the introduction of Fourier transform meth-
streamer of receivers. As it moves one-half receiver interval
ods. These frequency-domain approaches proved to be
along a seismic line, it fires a shot and records the pressure
at each receiver location. A trace is associated with each
Manuscript received May 3, 1984; revised May 16, 1984. shot and receiver point.
J. Gazdag is with IBMScientific Center, Palo Alto, CA 94304, USA. Let r be the horizontal coordinate of the receiver and s
P. Sguazzero is with IBM Scientific Center, 00149 Rome, Italy. be the horizontal coordinate of the source. Both are mea-

o o 1 8 - 9 2 1 9 / 8 4 / 1 ~ - 1 3 0 2 ~ . M )Ol98-4 IEEE

1302 P R O C E E D I N G S OF THE IEEE, VOL. 72, NO. 10.O C T O B E R 1984

 reconstruction of the wavefield at the time when the reflec-
RECEIVERS SHOT tiontook place. In this process, wave energy is being
COMMON
OFFSET
COMMON
MIDPOINT COMMON OFFSET "migrated" over the ( r , s) plane shown in Fig. 1 .
For practical and economical reasons, conventional
seismic processing techniques have been developed for
groups of traces, called gathers, aligned parallel with one of
the four axes shown in Fig. 1. A set of traces with the same
source ( s = constant) forma common-shot gather.The
traces are alsorearranged into common-midpoint (CMP)
gathers or common-offset gathers as shown in Fig. 1. The
reason is that the two standard operations that are responsi-
ble for the transferring of wave energy to its proper position
are stacking of the CMPgathersand migration ofthe
resulting stacked CMP section. Stacking shifts energy with
I t
RECEIVER COORDINATE r respect to the offset axis, and migration shifts energy with
Fig. 1. Relationship among the horizontal coordinates r, 5 , respect tothemidpoint axis.Thereason for performing
x , and h. All axes represent distances measured along the these operations separately is economic. Ideally, they should
seismic line. Each doton the surfacecorresponds to a be performed together as discussed in Section VII.
seismic trace.
B. Stacking
sured along the seismic line. However, for mathematical Stacking consists of the summation of the traces of each
reasoning, it is helpful to represent the receiver coordinate CMP gather after correcting them to compensate for the
r a n d the source coordinate s on orthogonal axes,as shown offset between source and receiver.This is known as the
in Fig. 1. We also define the midpoint coordinate between normalmoveout (NMO) correction. The NMO process is
+
source and receiver as x = ( r s)/2, and the source/re- illustrated onthe idealized model shown in Fig. 3. The
ceiver half-offset coordinate as h = ( r - s)/2. From these
equations we see that x and h are another set of axes
rotated 45" with respect to the axes r and s.
To understand the reason for migrationlet us consider
the ray paths associated with two point reflectors shown in
Fig. 2.Each trace includes two wavelets at times t measured

I T T T T X

Fig. 3. Illustration of NMO correction and stacking. The

traces of the CMP gather are summed along the NMO
hyperbola as indicated. The normalized stackedtrace is
plotted at h = 0.

wavelets indicate reflected events from the horizontal plane.

Fig. 2. The seismic concept illustrated with two point re- Their two-way travel time is given by
flectors. The recorded traces form a common-shot gather. 4h2
The loci of the recorded wavelets represent two-way travel t2( h ) = t: +-
time as indicated. l
, J

assuming that the velocity of wave propagation v is con-

stant. This equation is the NMO hyperbola. The difference
along the raypaths from source to object and from there to between the reflection time observed for an offset 2h and
the receiver. This example clearly demonstrates that a trace that corresponding to h = 0 is the NMO correction, i.e.,
consists of a superposition of seismic waves propagating
from all possible directions. Consequently, the information AtNMo t( h ) - to
regarding the cross section of reflectivity must be inferred When this amount of time shift is applied to the wavelets
from the interrelationships among an ensemble of traces. shown in Fig. 3, apart from a minor distortion effect, they
To do this, we need efficient and accurate computational appear as if they were recorded with coincident source and
techniques to map the wavelets onto the location of the receiver, i.e., h = 0. The NMO corrected tracesare then
reflecting objects. This mapping canalso be viewed as a summed, or stacked. The ensemble of stacked traces along

C A Z D A C A N D S C U A Z Z E R O : M l C R A T l O NO F SEISMIC DATA 1303

the midpoint axis is referred to as the CMP stacked section. method, which is based on an integral representation of the
A very important benefit of this operation is the significant solution of the wave equation. First, however, a geometric
improvement in the signal-to-noise ratio of the CMP sec- view of wavefield reconstruction is studied in some detail,
tion in comparison with the unstacked data. though at an introductory level. Our guiding principle is to
conveythe concepts ofmigration in a simple descriptive
C. The Exploding Reflector Model manner. The specific techniques used to implement these
ideas in the form of practical and efficient computer pro-
A CMP section mayberegarded as data obtained from grams are not considered in this paper.
coincident sources and receivers ( h = 0), as shown in Fig.
qa). In this zero-offset model, the energy travel path from
A. Geometric View of Wave Reconstruction
source to reflector is identical to that from reflector to
Let us consider the space-time distribution of energy in a
simple seismic experiment shown in Fig. 5. The ( x , z ) object

lX plane represents the cross section of the earth, and the

( x , t ) image plane is the record section. In the example
shown, there is one point diffractor in theobject plane.
When such a point diffractor is excited by a surface source,
it sets off outgoing waves in all directions. Using the ex-
ploding reflector model, such a diffractor canbe approxi-
mated by a buried point source. Using v/2 as the velocity

(4
Fig. 4. (a)The CMP section may be regarded as data ob-
tained with coincident (zero-offset) sources and receivers.
(b) The exploding reflector model of the CMP data for
migration.

receiver. The assumption is that all sourcesare activated

simultaneously, but each receiver records signalsoriginating Fig. 5. The space-time distribution of the energy of a point
source activated at t = 0 forms a cone. The intersection of
from the same source-receiver point. this cone with the image plane z = 0 defines the diffraction
Such zero-offset data do not correspond to any wavefield hyperbola.
resulting from a single experiment. As a result, it is helpful
to create a hypothetical physical experiment to provide an
intuitive picture of zero-offset migration. Suchan experi- of wave propagation, the waves observed and recorded at
ment is known as the “exploding reflector model” [I]. In z = 0 result in a zero-offset image of the diffractor.
this model, shown in Fig. qb), the energy sources are not at Assuming that v = constant, after the point source is set
the surface, but they are distributed along thereflecting off at t = 0, the resulting wavefronts are concentric circles,
surfaces. In other words, the reflectors are represented by which expand with increasing time, as shown in Fig. 5.
buried sources, which are activated at the same time t = 0. From this diagram it is seen that the space-time ( x , z , t )
Therefore, one needs to be concerned only with upward distribution of the wave energy lies on a cone. The intersec-
traveling waves. Since the record section involves two-way tionof this cone with the z = 0 plane is a hyperbola
travel time, it needs to be converted to oneway travel expressed by
time. In practice, thetime scale of CMPsections is kept
unchanged. Instead, the velocity of wave propagation is
scaled down by a factor of two.
The upward propagating part of the energy from the point
source lies along a hyperbola in the record section. The
FOR MGRATIONO F CMP DATA
I\\. INTEGRAL METHODS
objective of migration is to recreate from the observed data
With the help of the exploding reflector model outlined in the image plane the situation that existed in the object
above, migration of CMP data can be defined as the map- plane at t = 0. In other words, migration may be regarded
ping of the wavefields recorded at the surface back to their as aprocedure of mappingthe waveenergy distributed
origin at t = 0. The basic theory of wavefield reconstruction along diffraction hyperbolas into the loci of the correspond-
is due to Hagedoorn [2]. Based on hisideas, graphical ing point sources.
methods were developed, which were later implemented Another way to establish the relationship betweenthe
ondigital computers. Our objective is to describe two recorded image and the object i t represents is to examine
simple classical methods followed by that of the Kirchhoff the significance of a point, say a delta function,inthe

1304 PROCEEDINGS OF THE IEEE, VOL. 7 2 , NO. 10, O C T O B E R 1%

 inout t r a c e s

Fig. 7. Diffraction summation migration.The input traces

are summedalongthe diffraction hyperbolacorresponding
to a single scatterer at location x to produce an output trace
Fig. 6. The space-time distribution of the energy of an
at each location x .
inward propagatingcircular wavefront (centered at z = 0)
lies on a cone. A t time Tall the energy is focused in the
neighborhood of a point. The reverse of thisprocessforms
the basis of wavefront migration. A counterpart of diffraction migration is wavefront inter-
ference migration, which can be understood from the geo-
metrical view shown in Fig. 6. This approach to migration is
image space. Consider thehalf-circle centered at z = 0 based onthe observation that each pointin the image
shown in Fig. 6. The initial conditions at t = 0 are arranged plane must be mapped into a semi-circle in the object
so that the acoustic energy distributed uniformly along the plane. Thus the wavefront method [3] is a “one to many”
half-circle propagates inward, forming smaller concentric mapping. For the constant velocity case, the procedure is
semicircles, until at t = T the energy is focused in one the following. We take a sample of a seismic trace in the
point. Thus the energyobserved atsome t = T on the record section and distribute it over a semi-circle (equal
image plane may be thought of as being originated from a travel times) in the object plane as illustrated in Fig. 8.
wavefront lyingon a half-circle centered at z = 0 and
having a radius r = vT/2, as shown in Fig. 6. Consequently,
in the process of reconstructing the wavefront of t = 0, a i n p u t traces
point of the image plane must be mapped onto a half-circle
in the objectplane. This means that a point diffractor in the
subsurface would cause a hyperbolic seismic event, but a
point on a seismic section must have been causedby a
half-circle reflector in the subsurface.

B. Classical Methods of Migration

We haveseen that for a constant-velocity medium the

zero-offsetrecord section of a point reflector is char-
acterized by a diffraction hyperbola. Thus to account for all
the energy due to the point source in estimating its in-
tensity, weperformthe following summation. Foreach
( x , z ) point of the object plane, we construct a diffraction outputtraces
hyperbola in the image plane. We determine where the Fig. 8. Wavefrontmigration. Each input trace is mapped
hyperbola intersects each trace. Then we take the value of onto a semi-circle corresponding to theproper wavefront.
each trace at the point of intersection and sum all these
values together. The result of this summation is taken as the The diffraction summation and the wavefront inter-
value ofthe migrated section at (x,z), and this value is ference methods had notable success in the late 1960s and
placed in the object plane at that point. This method is early 1970s [4].They alsohave a number of undesirable
called diffraction summation (or diffractionstack) migra- characteristics [3, p. 1241. The reason for their shortcomings
tion and represents a straightforward, albeit not theoreti- is rooted in the fact that while these migration procedures
cally sound, approach. Fig. 7 shows how, in diffraction make good sense and are intuitively obvious, they are not
summation migration, an output trace is generated from the based on a completely sound theory.
input traces. The input traces represent the stacked section,
whereas the output traces define the depth section.
C. The Kirchhoff Method
It should be noted that, in practice, depth is most often
measured in unitsof the two-way vertical travel time, The Kirchhoff integral theorem expresses the value of the
T = 2z/v. Under the constant velocity assumption, T corre- wavefield at an arbitrary point in terms of the value of the
sponds to the apex of the diffraction hyperbola. Thus dif- wavefield and its normal derivative at all points on an
fraction summationmigration amounts to summing wave arbitrary closed surface surrounding the point [5, p.3771. In
amplitudesalongthediffraction hyperbola and mapping practice, measurements are made only at the surface, there-
the result at the apex of the same hyperbola. fore, instead of a closed surface, the integration must be

CAZDAC A N D SCUAZZERO: MIGRATION OF SEISMIC DATA 1305

limited to the surface of the earth. furthermore, in seismic arranging these traces into some kind of gathersor con-
practice, only the wavefield is recorded, and its normal stant-offset sections. The computer generation of seismic
derivative is not axtailable.The Kirchhoff integral canbe data calls for the solution of the equations governing the
expressed more easily for three-dimensional data, i.e., where wave propagation in the mediumdefined in themodel
the pressure wavefield is defined over some area of the under consideration. It is important, however, to ensure
earth’s surface rather than along a line as we assumed in the that the modeling is not simply an inverse of the migration
two-dimensional case. Let p(x, y,z = 0, t ) be the recorded from the numerical and algorithmic point of view. Other-
data, where x and yare surface coordinates, z is depth, and wise, errors in the forward process may be partially canceled
t is time. Then the migrated wavefield atsome point in the inverse process, which can lead to incorrect error
(xl, y,, z,, t) is obtained from the following expression [6], estimates. To avoid these pitfalls, the synthetic seismic
V I
,[8, P. 1251, PI,as: records presented in this paper are obtained by means of an
algorithm that permits the computation of two-way reflec-
P(xly,,zl,t= 0) tion times with precision limited only by roundoff error.
=jj=?g (x,y,z= O , t = r/c)dxdy (4) In our approach to forward modeling, the reflecting ob-
jects are represented by scattering surfaces in three dimen-
where sions, or lines in two-dimensions. The utility of approximat-
ing acontinuousreflecting sheetby discrete scattering
r = [(x, - x)’ +( y, - y)2 + z:] 1 /2
points.in a scaled seismic test tank was shown by Gardner
et a/. [ I l l . The computer implementation of this method is
c =. v/2 is the half velocity (assuming that t is the two-way based upon the ray theory of diffraction and the superposi-
travel time), and 9 is the angle between the z axis and the tion principle. The modeling of a complex subsurface re-
line joining (xl, y,,zl) and (x,y,z= 0). flectivity structure is done by representing it as a set of
To obtain the two-dimensional version of (4), one must
independent point scatterers. To obtain the zero-offset sec-
assume that p is independent of y. After integrating with tion of an elementary diffractor, .following the exploding
respect to y, one obtains the following asymptotic ap- reflector principle, rays are traced in every direction through
proximation [8, p. 1261: the medium up to the surface. A wavelet whose amplitude
is proportional to the strength of the diffractor is assigned
~ a 9: / ~ p ( x , z = O , t = ~ / c ) d x
p ( x , , ~ ~ , t = O ) - / -cos
G (after taking into account geometrical spreading) to the grid
point of the time section corresponding to the arrival of the
(5) ray at the surface. The contributions of all point diffractors
where a:/2
denotes the half derivative with respect to t. A are then summed to obtain the zero-offset section of the
simple expression for the operator a:/’ can be given in the structure.
frequency domain, where the transfer function of is An example of this modeling approach is Model 1 repre-
given by senting a syncline shown in fig. ?a). The synthetic zero-off-
set section of Model 1 is shown in Fig. yb), assuming a
uniform velocity of v = 3 km/s. Fig. 10 illustrates the same
It should be observed that, apart from the half-derivative
operator and a weighting factor, ( 5 ) represents an integra-
tion alongthe diffraction hyperbola. In this respect it is x (km)
0 2 4 6 8 1 0 1 2
quite similar to thediffraction summation method. It is r
important to note, however, that the Kirchhoffmethod
yields significantly higher quality migration than the diffrac-
tion summation approach discussed earlier.
When the propagation velocity c varies with the space
coordinates x, y, and z, (4) and (5) must be suitably gener-
:j
3 v
alized. for instance, in two dimensions we obtain [IO]

P(Xl,Zl,t = 0)

=jWa;/*p(x,z=O,t= T(x,,z,,x,z=~))~x (7)

where W is a weighting factor, and T( x,, z,, x , z) represents

thetwo-way travel time measured along a raypath that
connects the points( xl,zl) and (x,z) obeying Fermat’s least
time principle [ 5 , p. 1281.

0.Synthetic Seismic Records

As we havediscussedearlier, migration is an inverse
process. To test a migration scheme and evaluate its perfor-
mance, we need a set of seismic data,e.g., a zero-offset (4
section obtained from an idealized model with known Fig. 9. (a) Schematic of representinga sync,ine in
reflectivity and velocities. This is usually done by simulating a constant velocity medium, (b) Synthetic zero-offset time
the forward process, collectingthe results into traces, and section of Model 1.

1306 PROCEEDINGS O F THE IEEE, VOL. 72, NO. 10, OCTOBER 1964
 xlkml canexpress the behavior of the waves in terms of the
3.0 4.0 5.0 6.0 7.0 8.0 9.0
source coordinates as

Assuming that the velocity function has no lateral Lfaria-

tions, that is, it does not depend on r and s, we may
Fourier-transform (8) and (9) with respect to r , s, and t. This
operation gives
Fig. 10. Migrated section obtained from the zero-offset
section of Model 1 (Fig. 9(b)) bymeansof the Kirchhoff
method.

and
section after migration using the Kirchhoff summation
method. We note that all the waveenergy is distributed
along the syncline, which demonstrates the high accuracy
that can be expected from the Kirchhoffmigration ex-
where k , and k, are the wavenumbers (or spatial frequen-
pressed by (4) and (5).
cies) with respect to r and s, respectively, and w is the
temporal frequency. P ( k , , k,, w , z r ,z,) is the Fourier trans-
Iv. WAVE-EQUATION
MIGRATION
IN IMAGE SPACE
form of p( r , s,t, z,, z,).
Since we are dealing with wave phenomena, as one can Equations (IO) and (11) have two independent solutions
reasonably expect, migration canalso be described as a that propagate in opposite senses in ( r , s, t) space. One
method of obtaining the numerical solution of some par- moves in the positive t-direction, the other in the negative
tial-differential equation. These partial-differential equa- t-direction. These pairs of solutions are generally referred to
tions can be solved numerically either in the imagespace as upgoing and downgoing waves [12, p. 1691, [13, p. 4281.
or in the object space. For historical and practical reasons, We reconstruct the wavefield by lowering the source and
migration schemes formulated in image space are consider- recorder coordinates to thelocationofa reflector under
ably more popular than those defined in object space. A consideration [14]. Coincident source and recorder posi-
variety ofmigration techniques based on solving a wave tions imply t = 0. Therefore, the only meaningful solution
equation over a computationalgrid in image spacehas is the one that returns to zero time, the one that represents
evolved in the last decade or so. These methods are gener- regressing waves whenmovingdownward. The one-way
ally referred to as wave-equation migration. This name can wave equations that govern these regressive waves are
be misleading since it seems to suggest that the Kirchhoff
method is not based on wave theory. To be sure, all
migration methods with solid theoretical foundation must
be, by definition, wave-equation methods. and
In what follows we shall derive the equations for wave
extrapolationunder the simplifying assumption that the
velocityofthe medium has no lateral dependence. The
treatment of laterally varying media will be presented in the where we have introduced
latter part of the paper. The equations under consideration
are for extrapolation of seismic data obtained from multiple
sources and multiple recorders as illustrated in Fig. 1.
to simplify the notation.
A. Equations for Wavefield Extrapolation When recorders and sources are lowered simultaneously,
their position may be denoted by a common depth coordi-
We shall follow the notation used in describing the data nate z, i.e., one may let
acquisition illustrated in Fig. 1. We let r and s be the
horizontal coordinates of the receiver and the source, re- z, = z and z, = z. (15)
spectively. We also definethe midpoint coordinate x =
This permits us to write
(r + s)/2 and the half-offset h = ( r - s)/2. The variable z
represents depth into the ground. We let p ( r , s,t, z,,z,)
represent the wavefield, where z, and z, are depth coordi-
nates for recorder and sources, respectively. For fixed
sources, the wavesseen by the recorders are governed by Substituting (12) and (13) into (16) we obtain
the scalar wave equation

which expresses the simultaneous downward extrapolation

in which v stands for the velocity of wave propagation. The of recorders and sources. This expression is known as the
reciprocity principle permits the interchange of sources and doublesquareroot equation [15]. The solution of (17) can
receivers. Using this principle and fixing the receivers we be expressed as

C A Z D A C A N D SCUAZZERO: MIGRATION OF SEISMIC D A T A 1307

 v.NUMERICAL FOR ZERO-OFFSET MIGRATION
METHODS IN
SPACE
IMAGE

+(I - K:)~’~] Az) Migration of zero-offset seismic data in image space can
be summarized as follows. The zero-offset data areex-
If thedata are defined in the (x, h ) domain, it is desirable to trapolated from the surface downward to some depth z =
express K, and K, as functions of k, and k, which are the nAz in n computational steps. The wave extrapolation calls
wavenumbers with respect to x and h, respectively. Using for the numerical solution of (26) or some approximation
the relationship among receiver and source coordinates r thereof. The computations are performed in the image
and s, and midpoint and half-offset coordinates x and h, space, using the zero-offset section as initial condition. In
that is this process, as z increases, recorded events are shifted
laterally toward thier correct position as they move in the
x=(r+s)/2 h=(r-s)/2 (19) negative t direction.Downward extrapolation to depthz
we find results in awavefield that would have been recorded, i f
both sources and recordershad been located at depth z.
Thus events appearing at t = 0 areat their correct lateral
position. Therefore, the extrapolated zero-offset dataat
Once the surface data have been extrapolated by means t = 0 are taken as being the correctly migrated data at the
of (18) to the entire half-space z > 0, migration can be current depth. These data ( t = 0) are then mapped onto the
accomplished by imaging the wavefield at t = 0, r = s. depth section atz, the depth of extrapolation. This map-
Thus the correctly migrated zero-offset data are given by ping process is also referred to as imaging. Since imaging is
a standard procedure, migration schemes differ from each
p(x,h=O,t=O,z)=~~~P(k,,kh,w,z)exp(ik,x).
other in theirapproach to wave extrapolation, the complex-
kz kh *
ityofwhich depends largely onthe migrationvelocity
(21) function.
If the migration velocity has no horizontal dependence,
B. Wave Equation for Zero-Offset Data
the extrapolation of zero-offset seismicdatacanbe ex-
In current seismic practice, chiefly for economic reasons, pressed by the exact wave-extrapolation equation (26) in
migration is performed after NMO and stacking on CMP the wavenumber-frequency domain. This equation has a
sections. The N M O correction transforms data into zero- simple analytic solution, whose implementation calls for a
offset under the assumption that phase shift applied to the Fourier coefficients of the zero-
offset section. In the presence of lateral velocity variations,
p(x,h,t’,z=O)=p(x,h=O,t,z=O) (22)
the exact wave-extrapolation equation (26) is no longer
inwhich t’ differsfrom t by the amount ofthe NMO valid. To circumvent this problem, the exact expression can
correction, i.e., be approximated by truncated series expansions [16],[18],
t’ = t + AtNMo (23) which can accommodate horizontalvelocity variations.
These equations are then solved numerically, either in the
where AtNMo is defined by (2) and illustrated in Fig.3. If space-time domain or in the space-frequency domain.
we let p denote the NMO corrected version of p, then in
view of (22), p can be regarded as being independent of h A. Migration in the (k,, 0) Domain
P(x,h,t,z=O)=p(x,h=O,t,z=O). (24) Letting p(x, t,z = 0) represent the zero-offset section and
Therefore, its Fourier transform is zero for all nonzero P(k,,u,z = 0) its double Fourier transform, the wave-
wavenumbers kh, i.e., extrapolation equation (26) becomes
~ ( k , , k h , w , Z =f O
Or)kI O
h #, o . (25)
Note that (17) is still the correct extrapolation equation for
F. Substituting (20) into (17) and using (25) we obtain The problem is to reconstruct the wavefield p ( x , t = O,z),
which existed at t = 0 from the wavefield p ( x , t , z = 0)
observed at z = 0. We shall assume that within one ex-
+
trapolation step, say from depth z to z Az, the velocity is
The one-way wave equation (26) is the fundamental equa- constant, i.e.,
tion for downward extrapolation of zero-offset data. It is
expressed in the wavenumber-frequency domain (k,,o),
v ( S ) = constant, zg S<z + Az. (28)
and does not havean explicit representation in the mid- Then the solution of (27) can be expressed as
point-time domain ( x , t ) . We have obtained it as a special P(k,,o,z+ Az) = P(k,,o,z)
case of a general extrapolation equation for multioffset
data. It can alsobeeasily derived by making use ofthe
exploding reflector model of zero-offset data as shown
elsewhere [16], [19]. It should be observed that the velocity
variable that figures in (26) is the half velocity, as one would (29)
expect, from the exploding reflector model for zero-offset in which P( k, a,z) is the zero-offset section extrapolated
data. to depth z. This analytic solution states that P is extrapo-

1308 PROCEEDINGS O F THEIEEE, VOL. 72, NO. 10, OCTOBER 1 9 8 4

 +
lated from z to z Az by simply rotating its phase by a 8. Migration in the (x, w ) Domain
specified amount. Therefore, migration schemesbased on
As we have seen, in the absence of horizontal velocity
this principle are referred to as phase-shift methods [19],
dependence, the set of independent ordinary differential
[20]. From P( k,, o,z) we obtain p ( x , t = 0 , z ) by inverse equations (27) governs the extrapolation of the zero-offset
Fourier transform with respect to k, and a summation with
seismic data in the (k,,w) domain. The simple analytic
respect to o,i.e.,
solution expressed by (29) is not valid for velocity fields
p(x,t = 0,z) =~~P(k,,w,z)exp(ik,x). (30) with lateral variations. In this case, the square-root expres-
k 0 sion in (27) must be approximated in some form, for in-
Although (29) is valid only under the assumption expressed stance, by a quadratic polynomial (Claerbout’s [22] para-
in (28), the velocity can vary from one Az step to another. bolic approximation). The expansion can then be formally
Thus the phase-shift method can accommodate media with converted into a numerical method in the ( x , o ) domain.
vertical velocity variations, i.e., in which v = v(z). After the approximation, the velocity v is allowed to vary
Undercertain circumstances, it may be acceptable to with x as well as with z, and a wave-extrapolation equation
migrate with a velocity that is constant over theentire for general media is obtained.
domain (object plane) of interest. This uniform migration Another alternative is rational approximation by trun-
velocity assumption, if properly exploited, enables one to cated continued fractions [23]. In this latter approach, an
develop a very efficient migration algorithm. This migration approximation of (27) is
scheme, often referred to as the F-K method, was first
reported by Stolt [21]. Essentially, it is a fast “direct” method
for the evaluation of (29) and (30). Let us rewrite (29) using
+
z and z = 0 in place of z Az and z, respectively, i.e.,
P(k,,o,z) = P(k,,o,z=O)exp(ik,z) (31)
where (36)
A practical approach to solving extrapolation equations like
(36) is by Marchuk splitting. This is done by decomposing
(36) into two extrapolators
which represents the dispersion relations of the one-direc-
tional wave equation (27). The migrated section p ( x , t = (37)
0,z) can be obtained from P(k,, o , z = 0) as follows:
which is known as the thin lens term, and
p(x,t=O,~)=//dodk,P(k,,w,z=O)

eexp { i [ k,( w ) z + k,x]} (33)

which corresponds to (30), except that the summations are which is the Fresnel diffraction term. Advancing to greater
replaced by integrations. In order to take advantage of the depths i s done by applying (37) and (38) alternately in small
fast Fourier transform (FFT) algorithm, we want to integrate Az steps. Equation (38) canbeexpressed in the ( x , ~ )
(33) with respect to k, instead of a. Using (32) we can domain first eliminating fractions and then Fourier trans-
rewrite (33) as forming it with respect to k,. The result is

p( x , t = 0 , z ) = // dk, dk, i(k,k,) exp { i [ k,z + k,x]}

(34)
where At each step Az, (37) advances P through the multiplication
by a phasor corresponding to a shift backward in time of
Vk, magnitude AT = 2Az/v. To advance P with (39), for
4k,,k,) = P(kx,w,z= 0). (35)
2[ k t + k:]’” numericalstability considerations, an implicit Crank-
Nicolson difference scheme is most often used.
In summary, the F-K method consists of three steps: The approximationof the square-root operator in (27)
1) the Fourier transformation of p ( x , t , z = 0) to obtain withthe rational operator appearing in (36) and the
finite-difference implementation of (39) alter the dispersion
P(k,,w,z = O),
2) theinterpolation and scaling of the Fourier coeffi- relations ofthe wave equation andgenerate numerical
cients as expressed in (35), errors increasing with the expansion parameter sin9 =
3) the inverse FFT of the interpolated results as expressed vkX/2w, where 9 represents the angle with respect to the
vertical axis of the wavefront under consideration. When
in (34).
the dip of the imaging wavefield exceeds 40“-45’, disper-
Migration with this fast method is limited to homogeneous sion effects separate the low frequencies fromthehigh
media with constant velocity. To overcome this limitation, frequencies during the wavefront depropagation. As a re-
Stolt [21] suggested coordinate transformations to cast the sult, the image of steeply dipping reflectors may be accom-
wave equation in a form that is approximately velocity paniedby “ghosts.” The phenomenon is demonstrated
invariant. through a numerical example, themodel ofwhich is

C A Z D A C A N D S C U A Z Z E R O : M I G R A T I O N O F SEISMIC DATA 1309

depicted in Fig. 11. The synthetic zero-offset section of the x(krnl
0.0 20 4.0 6.0 12.0
8.0 10.0
model is shown in Fig. 12. The migrated depth section 0.0
obtained by a finite-difference implementation of (37) and
(39) is illustrated in Fig. 13. Note that the non-steeply
dipping events are clear, but the steeply dipping event has 2.0
“ghosts.“
Better dispersion properties can beobtained including
higher order terms in the formal expansion [16]. An alterna-
4.0
tive approach to wavefield extrapolation in general media is
described in [24]. A t each Az step, the wave extrapolation is
accomplished in two stages. In the first stage the wavefield
P( x, w , z) given at depth z is extrapolated to z Az by the + 6.0

phase-shift method using G reference velocities v,,y;.., v, Fig. 14. Depthmigration section obtained from the time
spanning the velocity range at depth z. This stage generates section shown in Fig. 12 using the phase shift plus interpola-
tion method. All reflectors are migrated correctly.

9
0 L 4 6 8 10 12
x (km)
G reference extrapolated wavefields at z A z , namely, +
v = 2.0 4, e,---, P(. in the second stage, the definitive wavefield
P ( x , w , z + A z ) is constructed byinterpolationof the G
2 r reference wavefields. Thisphase shift plus interpolation
v >
= 3.0
(PSPI) method is unconditionally stable and has good dis-
41 v - 4.0 persion relation properties worth its relatively high compu-
tational cost (for each frequency w , and for eachstep in
depth, G+ 1 Fourier transforms in the x direction are re-
i quired). The synthetic zero-offset section of Fig. 12 migrated
6 %-( km) with the PSPI method is shown in Fig.14. One can notice
Fig. 11. Schematic of model 2, representing a dipping mul- significantimprovement in themigration of the steeply
tilayer example. The thick-line segments indicate where the dipping reflector in comparison with the one shown in
reflector segments have been ”turned on.” Fig. 13.

xlkrn) C. Migration in the (x, t ) Domain

0.0 2.0 4.0 6.0 8.0 10.0 12.0
Migration as a wave depropagation process is best de-
scribed in the “transformed” domains ( k , , w ) and ( x , @ ) . It
was introduced, however, in the early 1970s using paraxial
wave equations in the “physical” space-time ( x , t ) [22],
[14]. We will reproduce here those equations, deriving them
from their frequency-domain counterparts.
An inverse Fourier transform of the wavefield P ( x , o , z )
in (37) and (39) results in the following two equations in
terms of p ( x , t , z ) :

Fig. 12. Zero-offset time section ofmodel 2 a+--=--.

v2 v a3p a3p
(41)
at2az 4 ax2& 16 a x 2 a z
x(krnl
0.0 2.0 4.0 6.0 8.0 10.0 12.0 For gentle dips, theright-hand side of (41)can be ne-
glected, and it becomes
aZp azP
wz+--=0.
4 a2
Equations (40) and (41) or (40) and (42) are specialized wave
equationsthat depropagate waveenergy only within a
small angle about the vertical axis, hence are called paraxial
approximations: details on their numerical integration with
finite-difference methods can be found in [12, pp. 211-2121.
Seismic interpreters oftenfindit useful toworkwith
subsurface migrated sections with depth measured in units
Fig. 13. Depth migration obtainedfrom the zero-offset time of time,for ease of comparison with unmigrated CMP
section of model 2 shown in Fig. 12. The ( x , o ) domain
finite-difference migration methodusedinthisexample sections. Such maps are obtained by introducing a pseudo-
corresponds to (37) and (39). Note that the steeply dipping depth r (vertical travel time) which is related to true depth
reflector is not imaged correctly. z by

1310 PROCEEDINGS O F THE IEEE, VOL. 72, NO. 10, OCTOBER 1984
 t = T , rather than at t = 0. Time migration methods based
(43) on thisimagingprinciple are often described as being
solved in a “downward moving” coordinate frame.
where V(z) is an approximation, independent of x , of the
velocity v ( x , z ) . With the introduction of the vertical coor-
VI. WAVE-EQUATION
MIGRATION IN OBJECT SPACE
dinate T and the help of an additional approximation, (40)
and (42) become Migrationof zero-offset data p ( x , t , z = 0) is accom-
plished by the methods described in Sections IV and V
through a downward extrapolation of surface data. Advanc-
ing in depth, image plane sections p ( x , t , z= const) are
computed and the final migrated section is given by the
extrapolated wavefield at time zero p ( x , t = 0,z). An alter-
Wavefieldextrapolation with (44)and (45) lead to time native approach calls for a reversetime extrapolation of the
migration as opposed to depth migration obtained with wavefield using the CMP section as a boundary condition.
(40) and (42). Time migration schemes depropagate the Marching backward in time, object plane sections (time-
wavefield as if the medium were horizontally layered and slices) p ( x , t = const, z ) are computed. Calculations begin
do not correctly refract rays at the velocity interfaces, but at time T, which defines the lastsample ofthe record
constitute a viable and economic tool in most practical section, and continue in the negative t direction until time
situations. An excellent example of finite-difference (time) zero. At that time the amplitudes in object space are con-
migration of a CMPstack section is illustrated in Figs. 1 5 sidered as the final migrated section.
and 16. Notice the decrease in width of the diffractions on In this section, dealing with reverse-time migration, for
the sake of clarity, we shall deviate slightly from the nota-
tions used throughout this paper. We shall let 9 ( x , t ) repre-
8 8 sent the stacked CMPdata,and p ( x , t , z ) the wavefield.
Consequently, p ( x , t = 0, z ) will represent the migrated
CMP section.
In object space, just as in image space, migration schemes
differ from each other only in the type of wave equations
being used for wave extrapolation. For example, Hemon
[28] and McMechan [29] suggested the “full” wave equation

stack TraceSpacing=33m which corresponds to our (8)or(9)expressed in the mid-

point variable x . This equation is then integrated numeri-
Fig. 15. CMP stacked section from data recorded i n the
cally with the following initial conditions:
Santa Barbara channel. Steepest dips are approximately 25”.
(From Hatton et a/. [17]; courtesy of Western Geophysical.) p(x,tT
=,z) =0, forzto (47)
and

q x , t = T,z) = 0.
at (4.8)
The boundary conditions at the surface are given by
p(x,t,z=O) = q(x,t). (49)
The problem with (46) is that it allows for reflections from
interfaces wherethe velocity changes rapidly. We recall
that the exploding reflector model shown in Fig. 4 assumes
that the surface record includes only upward propagating
M(
g- TraceSpacing=33rn waves. Therefore, in the depropagation phase all wave
Fig. 16. Finite-differencetimemigration of the CMP sec-
components should move downward in reverse time, and
tion illustrated i n Fig. 15. (From Hatton et a/. [17]; courtesy of noupward reflected waves are admissible. This problem
Western Geophysical.) can be easily avoided if one uses, instead of (46), a one-
directional wave equation. Baysal et a / . [30]and Loewenthal
et a/. [31] used an equation for downward propagating
therighthand of the migrated section and the develop- waves
ment of synclines at the left. When a better positioning is
required, after time migration, a correction is made using
Hubral’s image ray time-to-depth conversion, [25, p. 1031,
*
a t = ‘2( ~ ’ i k , [ l +( k x / k , ) 2 ] 1 ’ 2 . F ) p (50)

[26], [27]. The solution of (44) can be written as p( x , t , T + where 9 and 9-’ are operators representing the double
A T ) = p ( x , t + A7,7) which represents a uniform translation Fourier transform from the ( x , z ) domain to the ( k , , k , )
of the wavefield p along the negative t axis as T increases. domain and its inverse. Unfortunately, (50) cannot be writ-
There is no need to perform this time shift explicitly. ten without the transform operators 9 and s-’,since the
Equation (44) can be accounted for by imaging the data at multiplication in the Fourier space by ik,[l (kx/k,)2]1/2 +

C A Z D A C A N D SCUAZZERO: MIGRATION O F SEISMIC DATA 1311

cannotbe expressed in the ( x , z ) space. Further details The theoretical basis for migration of unstacked (multi-
regarding (50) and its numerical solution are found in [32] offset) data was discussed in Section IV. in this process, the
and [33]. The initial and boundary conditions for reverse- multioffset data are converted (see, for example, (17) and
time migration with (50) are given by (47) and (49), respec- (21)) into a migrated zero-offset section. This is a relatively
tively. Using a third-order Runge-Kutta algorithm, for ex- costly process and does not provide any intermediate re-
ample, to solve (50) on the zero-offset section shown in Fig. sults. To overcome these shortcomings of the “one-step”
12 generates the migrated depth section illustrated in Fig. migration of unstacked data, a number of special-purpose
17. before-stack migration techniques were developed. Forex-
ample, the prestack partial migration [IS] and the offset
dkml
continuation algorithm [37],[38]are applied to common-
0.0 2.0 4.0 6.0 8.0 10.0 12.0 offset sections after they have been NMO corrected. These
0.0
stack-enhancement [39] techniques share the advantage of
yielding an unmigrated CMP stack section that helps the
interpreter in resolving spuriousevents generated on the
20
final migrated section by inaccurate velocity estimates. They
are computationally more expensive than the conventional
processing since they introduce an additional process be-
4.0 tween N M O and stacking but leave themigration cost
unchanged.
The methods described in this section belongto the
6.0 former category in which the separate components of the
conventional processing (NMO, stacking, zero-offset migra-
Fig. 17. Depthmigration section obtainedfrom the time tion) are merged into a unified process. They represent the
section shown in Fig. 1 2 using reverse-time migration based
best methods available for imaging the subsurface in a
on (50).
variable velocity medium. They are, unfortunately, also the
most expensive in terms of computation and data handling.

Another approach to reducing reflections was reported

A. Migration of Multioffset Data
by Baysal etal. [35], who propose the use ofa two-way
nonreflecting wave equation. Their approach begins consid- The theoretical considerations for migration of multioff-
ering the general wave equation [34, p. 1711 set data were presented in Section IV. The wavefield ex-
trapolation is governed by (13, and the imaging is accom-
plished by implementing (21). The solution by phase shift
of (13, expressed by (18), can accommodate only vertical
where p is the density of the medium. Since p does not velocity variations from one Az step to another. In the
play any role in migration, itsvaluemaybe chosen to presence of lateral velocity variations, (17) needs to be
minimize reflection. This is done by setting p so that the generalized as follows [15]:
+ - ( I1
‘acoustic impedance
K( x , z ) = p ( x , z ) v ( x , z )= constant (52)
-ap
az
= -( 1I
v(r,z)
-I C y 2
v( st.)
- IC:)
l’*I
P

over the entire object plane. (53)

The basic conceptof using CMP data as a boundary
where IC, and IC,are given by (14). This impliesthatwe
condition and computing its “response” in the object plane
distinguish between the velocity alongthe shot axis and the
is not new. The heuristic method called wavefront inter-
velocity along the shot axis and thevelocity along the
ference migration discussed in Section Ill is based on the
receiver axis. We hasten to point out, however, that in (53)
same principle. The only important difference is that the
the algebraic expressions in the operators IC,and I C ~are not
reverse-time migration methods are based on a wave-equa-
immediately defined when v ( r , z ) and v( s,z) have lateral
tion solution rather than the intuitive approach depicted in
variations. This is the same limitation as the one we have
Fig. 8.
encountered with the single-square-root equation (27). For-
tunately, most ofthesolution methods for zero-offset
VII. MIGRATION BEFORESTACKING
migration are applicable to (53). In other words, multioffset
The task of imaging the subsurface in conventional data can be extrapolated by the PSPl method or by finite-
seismic data processing is partitioned intwo steps: the difference methods developed in the ( r , s , ~ )domain or in
stacking of the CMP gathers, and the migration of the CMP the ( r , s, t ) domain.
stacked section. This procedure, however, is correct only for Having decided on a suitable extrapolation scheme, be-
a simplified subsurface model consisting of horizontal layers fore-stack migration can be summarized as follows [ a ] . The
of laterally uniform velocity and ofhorizontal reflectors. algorithm operates on the entire data volume p(r, s, t,z = 0)
Levin [36] has shown that the stacking velocity depends not processing it in the alternating directions r a n d s. For each
only on the propagation velocity of the overburden but also depth level z, in increments of Az, beginning with z = 0,
on the dip angle of the reflectors. in structurally complex the following steps are executed:
formations, the correction for all dip anglesbecomes im- 1) Order the data into common-shot gathers p(r, so, t , z )
possible, and stacking loses its effectiveness. Under these and extrapolate downward from z to z Az each com- +
circumstances, it may be advantageous to consider migra- mon-source gather s = so using (12).
tion before stacking. 2) Reorder the data obtained from step 1 into common-

1312 PROCEEDINGS OF THE IEEE. VOL 72, NO. 10,OCTOBER 1 9 8 4

 xtkm)
receiver gathers r = ro and extrapolate downward by Az 0.0 1.0 2.0 0
.
3
each common-receiver gather r = r, using (13). 0.0 ~

3) Perform the imaging m( x, z) = p( r = x, s = x, t = 0, z),

where m(x,z) is the migrated CMP section.
0.5

6. Migration of Common-Shot Gathers

Migration before stacking canalso be accomplished by
migrating the individualcommon-shot gathers and sum-
ming these migrated gathers. The downward extrapolation
of the recorders is governed by (12) or some approximation 1.5
thereof. The extrapolated wavefield is focused at the reflec-
tor location, not at time zero, but at the propagation time
between the shot and the reflector. To simplify matters, let 20
us consider a single point reflector, for example Pl in Fig. 2,
Fig. 19. A syntheticcommon-shotgather obtained from
whose position (r',z') is specified in the receiver-depth model 3 with source location s = 1.5 km.
coordinate system.The event corresponding to the point
reflector 6 is observed at t = t, + t 2 , where tl(s; f , f ) is
the travel time between the source and the reflector, and x(km)
0.0 1.0 2.0 30
t 2 ( c t',z') is the travel time between the reflector and the
receiver. Under these assumptions, the downward extrapo-
lation focuses the recorded wavefield in the neighborhood
of the point ( r ' , ~ ' ) at time r,(s;r',z'), which is the correct
time for imaging.
We are now in a position to set forth an algorithm for
migration before stack. For each common-shot gather
p(r, F, t,z = 0) we execute the following steps:
1) Computethe arrival time of the primary wavefield
tl = tl($ I ' , z') at all points ( r ' , ~ ' )of the subsurface.
2) Extrapolate downward (for instance with (12)) in depth
the common-shot data p(r, F, t,z = 0), thus regenerating for
each point of thesubsurface the scattered field p ( / , 5, t,z').
3) Perform the imaging Fig. 20. Depth section obtained bymigration before stack.
First, the common-source gathers of model 3 with source
m ( f , f ; < )=p(t,5,t1(5;t,f),f). (54) locations ranging from 1 to 2 km are migrated. The section
shown is the result of the superpositionof 40 migrated
Then sum over a set of migrated common-source gathers, gathers.
that is,
rn(f,f)=xm(r',z';<). (55) shows the superposition of 40 migrated common-shot
E gathers spanning the interval between s = 1 km and 5 =
The procedure is illustrated in the following example. Fig. 2 km.
18 shows a model consisting offour layersseparated by
reflecting interfaces: the velocity of the medium ranges C Relative Merits of the Two Approaches
from 2.0 to 3.5 km/s. Fig. 19 shows the common-shot
The treatment of the multioffset seismicdata in the
gather that is obtainedsetting s = 1.5 km. Thereare 96
source-receiver coordinate system, as described above, has
traces with a spacing of Ax = 25 m, the sampling time is
a rather good theoretical foundation. Unfortunately, it in-
A t = 10 ms, and 256 samples define a trace. Finally, Fig. 20
volves very large amounts of data, the repeated transfer of
which from primary storage to secondary storage, together
x ( km)
0.0 0 . 5 1.0 1.5 22. 0
.5 3.0 with the necessary data transpositions, can become rather
0.0 I
v - 2.0
burdensome and taxing on the computing system.The
migration of individual common-shot gathers, on the other
0.5 i
1
v - 2.5
hand, requires relatively little storage space and the compu-
tations are simpler in comparison with the former approach,
l.O which call for the simultaneous migration of a large set of
-
~

1.5 v 3.0 common-shot gathers. The only serious problem for which
one must be prepared is related to imaging. The imaging of
2.0 i v - 3.5 the extrapolated wavefield, expressed by (54), requires the
knowledge of the arrival time of the primary wavefield. So
2.5 1 long as this arrival time rl( s; r, z) is a single-valued function
of r and z, there is no problem. If, however, there are
3.0-
z(km) multiple raypaths between the source s and some points of
the ( r , z ) domain, t, may assume multiple values, and there
Fig. 18. Schematicof model 3, representinga multilayer
example with a gentle dip in the second reflector. The is no longer a uniquely defined imaging time. In this case,
thick-line segments denote reflectors. following [14], the primary wavefield must alsobe com-

GAZDAC OF SElSMlC D A T A 1313

puted and the imaging mustbe performeddeconvolving [3] J. D. Johnson and W. S. French, “Migration-The inverse
the scattered wavefield by the primary wavefield. method,” in Concepts and Techniques in Oil and Cas Ex-
ploration, K. C. lain and R. J. P. de Figueiredo, Eds.Tulsa,
OK: SOC. Exploration Geophysicists, 1982, pp. 115-157.
VIII. CONCLUDINGREMARKS [4] W. A. Schneider, “Developments in seismicdata processing
and analysis: 1968-1970,” Ceophys., vol. 36, pp. 1043-1073,
We have presented an overview of the major advances in 1971.
seismic migration. Our aim was to provide a logical se- [5] M. Born and E. Wolf, Principles of Optics, 3rd ed. New York:
quence of development rather than a historical one. We Pergamon, 1%5.
have put emphasis on the fundamental concepts of migra- [6] W. S. French, ”Computer migration of oblique seismic reflec-
tion profiles,” Ceophys., voi. 40, pp. 961-980, 1975.
tion and have carefully avoided passing judgement or ap-
[q W. A.Schneider, “Integral formulation for migration in two
praising value based upon their current populatiry. The and three dimensions,” Ceophys., vol. 43, pp. 49-76,1978.
merit of a migration method is largely decided by economic [8] A. J. Berkhout, Seismic Migration. Amsterdam, The Nether-
factors. These factors change rapidly with the advancement lands: Elsevier 1980.
of recording technology, computer architecture, and cost- [9] A. J. Berkhout,”Wave field extrapolation techniques i n seismic
migration, a tutorial,” Ceophys., vol. 46, pp. 1638-1656, 1981.
performance indices of data processing systems. [lo] J. A. Carter and L. N. Frazer, “Accomodating lateral velocity
Mainly for reasons of cost, migration in the 1970s was changes i n Kirchhoff migration by means of Fermat’s princi-
appliedto stacked CMP sections, and lateral velocity ple,” Ceophys., vol. 49, pp. 46-53, 1984.
changes were not correctly taken into account in routine [ I l l C. H. F. Cardner, W. S. French, and T. Matzuk, “Elements of
migration and velocity analysis,” Ceophys., vol. 39, pp,
migration programs. With the continuing trend toward lower 811-825, 1974.
cost-performance ratio in advanced computer architec- J. F. Claerbout, Fundamentals of Geophysical Data Processing.
tures, high-accuracy migration methods can be expected to New York: McGraw-Hill, 1976.
gain acceptance in the near future. Another area that will E. A. Robinson and M. T. Silvia, Digital Foundations of Time
receive increased attention is migration before stacking. In Series Analysis: Vol2, Wave Equation Space Time Processing.
San Francisco, CA: Holden-Day, 1981.
addition to better migration, a major benefit resulting from J. F. Claerbout, “Toward a unified theory of reflector map-
migrationbefore stacking can be derived fromimproved ping,” Ceophys., vol. 36, pp. 467-481, 1971.
velocity information. Ideally, we would like to have a pro- 0. Yilmaz and J. F. Claerbout, “Prestack partial migration,”
cess that simultaneously images velocities and migrates Geophys., vol. 45, pp. 1753-1 779, 1980.
J. Cazdag, “Wave equation migration with the accurate space
data to correct subsurface locations. This calls for imaging derivative method,” Ceophys. Prospect., vol. 28, pp. 60-70,
part of the subsurface recursively step by step, using migra- 1980.
tion and velocity analysis procedures [41], [42). [I 71 L. Hatton, K. L. Larner, and B. S. Gibson, “Migration of seismic
Seismic exploration methods and processing techniques data from inhomogeneous media,” Ceophys., vol. 46, pp.
751-767,1981.
are growing rapidly. Innovations can be observedover a
A. J. Berkhout and D. W. Van W. Palthe, “Migration in terms
broad range from recording techniques to displays. Record- of spatial deconvolution,” Ceophys. Prospect., vol. 27, pp.
ing techniques have been developed to provide areal cover- 261 -2W, 1979.
age. Acquisition and migrationof seismicdata in three J. Cazdag, “Wave equationmigrationwiththe phase-shift
dimensions are becoming widespread. These methods are method,” Geophys., vol. 43, pp. 1342-1351, 1978.
C. Bolondi, F. Rocca, and A. Savelli, “A frequency domain
particularlyuseful in the presence of complex structures. approach to two-dimensional migration,” Ceophys. Prospect.,
Migration algorithms for three-dimensional data are natural VOI.26, pp. 750-772, 1978.
extensions of those used for two-dimensional data [43], [MI. R. H. Stolt, “Migration by Fourier transform,” Ceophys., vol.
Increasing demand for oil and decreasing availability of 43, pp. 23-48, 1978.
J. F. Claerbout, ”Coarse grid calculation of waves i n inhomo-
hydrocarbon deposits are the motivating forces behind most geneous media with application to the delineation of com-
innovative techniques which have one thing in common: plicated seismic structures,” Ceophys., vol. 35, pp. 407-418,
they are aimed at improving resolution. A decade ago, 1970.
geophysicists might have been satisfied with mapping the R. Clayton and B. Engquist, “Absorbing boundary conditions
gross structure of the subterrain. Today, they hope to de- for acoustic and elastic wave equations,” Bull. Seis. Soc.
Amer., vol. 67, pp. 1529-1540,1977,
termine rock parameters, hydrocarbon content, as well as J. Cazdag and P. Sguazzero, “Migration of seismicdata by
entrapment structure. These demands have stimulated re- phase shift plus interpolation,” Ceophys., vol. 49, pp. 124-131,
search in more general and theoretically appealing inverse 1984.
methods [45]-[52]. These works, when developed into prac- P. Hubral and T. Krey, Interval Velocities from Seismic Reflec-
tion Time Measurements. Tulsa, OK: Soc. Exploration Ceo-
tical processing tools, have the potential of providing in- physicists, 1980.
sight into long-standing unsolved geophysical problems. P. Hubral, “Time migration-Some ray theoretical aspects,”
Ceophys. Prospect., vol. 25, pp. 738-745, 1977.
K. L. Larner, L. Hatton, B. S. Gibson, and I-C. Hsu, “Depth
ACKNOWLEDGMENT
migration of imaged time sections,” Ceophys., vol. 46, pp,
734-750, 1981,
The authors wish to thank B. Gibson of Western Geo- [28] C. Hemon, “Equations d’onde et modeles,” Ceophys. Pros-
physical Company for providing Figs. 15 and 16. pect., vol. 26, pp. 790-821, 1978.
[29] C. A. McMechan, “Migration by extrapolation of time-depen-
dentboundary values,” Ceophys. Prospect., vol. 31, pp.
REFERENCES 413-420, 1983.
[30] E. Baysal, D.D. Kosloff, and J. Sherwood, ”Reverse time
[I] D. Loewenthal, L. Lou, R. Robertson, and J. W. Sherwood, migration,” Ceophys., vol. 48, pp. 1514-1524, 1983.
“The wave equation applied to migration,” Ceophys. Pros- [31] D. Loewenthal and I. R. Mufti, “Reversed timemigration i n
pect., vol. 24, pp. 380-399, 1976. spatial frequency domain,” Ceophys., vol. 48, pp. 627-635,
[2] J. C. Hagedoorn, ”A process of seismic reflection interpreta- 1983.
tion,” Geophys. Prospect., vol. 2, pp. 85-127, 1954. [32] J. Cazdag, “Extrapolation of seismic waveforms by Fourier

1314 PROCEEDINGS OF THE IEEE, VOL. 72, NO. 10,OCTOBER 1984

 methods,“ ISM 1. Res. Develop., vol. 22, pp. 481-486, 1978. [441 H. Jakubowicz and S. Levin, “A simple exact method of 3 - 0
-, “Modeling of the acoustic wave equation with migration-Theory,” Geophys.Prospect., vol. 31, pp. 34-56,
transform methods,” Ceophys., vol. 46, pp. 854459,1981. 1983.
L. M. Brekhovskikh, Waves in Layered Media. New York: J. K. Cohen and N. Bleistein, “An inverse methodfor de-
Academic Press, 1 W . termining small variations in propagation speed,” SIAM ].
E. Baysal, D. D. Kosloff, and J. Sherwood, “A two-way nonre- Appl. Math., vol. 32, pp. 784-799, 1977.
flecting wave equation,” Ceophys., vol. 49, pp. 132-141,1984. -, “Velocity inversion procedure for acoustic waves,”
F. K. Levin, “Apparent velocity from dipping interface reflec- Ceophys., vol. 44, pp. 1077-1087,1979.
tions,” Ceophys., vol. 36, pp. 510-516, 1971. S. Raz, “Direct reconstruction of velocity and density from
G. Bolondi, E. Loinger, and F. Rocca, “Offset continuation of scattered field data,” Ceophys., vol. 46, pp. 832-836,1981.
seismic sections,”Ceophys.Prospect., vol. 30, pp. 813-828, R. W. Clayton and R. H. Stolt, “A Born-WKBJinversion method
1982. for acoustic reflection data,” Ceophys., vol. 46, pp. 1559-1567,
S. M. Deregowski and F. Rocca, “Geometrical optics and 1981.
wave theory of constant offset sections in layered media,” M. Lahlou, 1. K. Cohen, and N. Bleistein, “Highly accurate
Ceophys. Prospect., vol. 29, pp. 374-406, 1981. inversion methods for three-dimensional stratified media,”
P. Hood, “Migration,” i n Developments in Geophysical Ex- SIAM J. Appl. Math., vol. 43, pp. 726-758, 1983.
ploration Methods-2, A. A. Fitch, Ed. London, England: F. C. Hagin and J. K. Cohen, “Refinements to the linear
Applied Science Publ. 1981, pp. 151-230. velocityinversion theory,”Ceophys., vol. 49, pp. 112-118,
P. S. Schultz and J. W. C. Shenvood, “Depth migration before 1984.
stack,” Ceophys., vol. 45, pp. 376-393, 1980. A. Bamberger, C. Chavent, and P. Lailly, “About the stability
P. S. Schultz, “A method for direct estimation of interval of the inverse problem in 1-D wave equations-Application
velocities,” Ceophys., vol. 47, pp. 1657-1671, 1982. to the interpretation of seismic profiles,” Appl. Math. Optim.,
J. Gazdag and P. Sguazzero, “Interval analysis by wave ex- VOI.5, pp. 1-47, 1979.
trapolation,” Ceophys. Prospect., vol. 32, pp. 454479,1984, P. Lailly, “The seismicinverse problem as a sequence of
6. Gibson, K. Larner, and S. Levin, “Efficient 3-D migration in before stack migrations,” presented at the Conf. on Inverse
two steps,” Ceophy. Prospect., vol. 31, pp. 1-33, 1983. Scattering, Univ. Tulsa, Tulsa, OK, M a y 1983.

G A Z D A C A N D S C U A Z Z E R O : M I G R A T I O N OF SEISMIC DATA 1315