Applications of Velocity-Stack

Methods to Seismic Data Processing

U.S. GEOLOGICAL SURVEY BULLETIN 2102

 AVAILABILITY OF BOOKS AND MAPS OF THE U.S. GEOLOGICAL SURVEY
Instructions on ordering publications of the U.S. Geological Survey, along with prices of the last offerings, are given in the current-year
issues of the monthly catalog "New Publications of the U.S. Geological Survey." Prices of available U.S. Geological Survey publications re-
leased prior to the current year are listed in the most recent annual "Price and Availability List." Publications that may be listed in various
U.S. Geological Survey catalogs (see back inside cover) but not listed in the most recent annual "Price and Availability List" may no longer
be available.
Reports released through the NTIS may be obtained by writing to the National Technical Information Service, U.S. Department of
Commerce, Springfield, VA 22161 ; please include NTIS report number with inquiry.
Order U.S. Geological Survey publications by mail or over the counter from the offices listed below.

BY MAIL OVER THE COUNTER

Books Books and Maps

Professional Papers, Bulletins, Water-Supply Papers, Tech- Books and maps of the U.S. Geological Survey are available
niques of Water-Resources Investigations, Circulars, publications over the counter at the following U.S. Geological Survey offices, all
of general interest (such as leaflets, pamphlets, booklets), single of which are authorized agents of the Superintendent of Docu-
copies of Earthquakes & Volcanoes, Preliminary Determination of ments.
Epicenters, and some miscellaneous reports , including some of thP.
foregoing series that have gone out of print at the Superintendent of
Documents, are obtainable by mail from • ANCHORAGE, Alaska-Rm. 101, 4230 University Dr.

U.S. Geological Survey, Map Distribution • LAKEWOOD, Colorado-Federal Center, Bldg. 810
Box 25286, MS 306, Federal Center
Denver, CO 80225 • MENLO PARK, California-Bldg. 3, Rm. 3128,
345 Middlefield Rd.
Subscriptions to periodicals (Earthquakes & Volcanoes and
Preliminary Determination of Epicenters) can be obtained ONLY • RESTON, Virginia-USGS National Center, Rm. 1C402,
from the 12201 Sunrise Valley Dr.

Superintendent of Documents • SALT LAKE CITY, Utah-Federal Bldg., Rm. 8105,

Government Printing Office 125 South State St.
Washington, DC 20402
• SPOKANE, Washington-U.S. Post Office Bldg., Rm. 135,
(Check or money order must be payable to Superintendent of West 904 Riverside Ave.
Documents.)
• WASIDNGTON, D.C.-Main Interior Bldg., Rm. 2650,
Maps 18th and C Sts., NW.

For maps, address mail orders to Maps Only

U.S. Geological Survey, Map Distribution Maps may be purchased over the counter at the following
Box 25286, Bldg. 810, Federal Center U.S. Geological Survey offices:
Denver, CO 80225

Residents of Alaska may order maps from • FAIRBANKS, Alaska-New Federal Bldg, 101 Twelfth Ave.

U.S. Geological Survey, Earth Science Information Center • ROLLA, Missouri-1400 Independence Rd.
101 Twelfth Ave., Box 12
Fairbanks, AK 99701 • STENNIS SPACE CENTER, Mississippi-Bldg. 3101
Applications of Velocity-Stack
Methods to Seismic Data Processing

By Myung W. Lee and Warren F. A gena

U.S. GEOLOGICAL SURVEY BULLETIN 2102

UNITED STATES GOVERNMENT PRINTING OFFICE, WASHINGTON : 1994

 U.S. DEPARTMENT OF THE INTERIOR
BRUCE BABBITT, Secretary

U.S. GEOLOGICAL SURVEY

Gordon P. Eaton, Director

For sale by U.S. Geological Survey, Map Distribution

Box 25286, MS 306, Federal Center
Denver, CO 80225

Any use of trade, product, or firm names in this publication is for descriptive purposes only and
does not imply endorsement by the U.S. Government

Library of Congress Cataloging-in-Publication Data

Lee, Myung W.
Applications of velocity-stack methods to seismic data processing
by Myung W. Lee and Warren F. Agena.
p. cm.~(U.S. Geological Survey bulletin : 2 102)
Includes bibliographical references.
Supt. of Docs. no. : I 19.3:2102
I. Seismology-Mathematical models-Data processing. I. Agena, W. F.
II. Title. Ill. Title: Velocity-stack methods. IV. Series.
QE75.B9 no. 2102
[QE539.2.D36]
557.3 s-dc20
[551.2'2'0 15118] 94-3663
CIP
 CONTENTS

Abstract ................................................................................................................................ .
Introduction .... ............... .. ...... ..... ...... ......... .. ........ .. ..... .. ......... .. ...... ...... ...... ..... ...... ......... .. ..... 1
Acknowledgments ........... ..................... ........... ........................................................ ............. 2
Theory................................................................................................................................... 2
Inverse and Transpose Operators ... ...... ..... ..... ............... ...... .. ......... ............ .. ..... ......... .. 4
Applications.......................................................................................................................... 6
Multiple Suppression.................................................................................................... 6
Model Description .......... .. ... ... ..... ......... ...... .. ...... ..... ...... ............... ........ ... .. ......... .. 6
Noise-Free Model................................................................................................. 7
Noise Model ....................... .......... .. ................... ........... ......................................... 8
Real Data Example ..... ...... ...... ... .. ...... ......... .. ...... .. ... .......... .. ....... .. .. ...... ......... .. ..... 10
Coherency Filtering .......... .. ....... ..... ..... ...... ... .. .......... .. ..... .. ............. ... ......... .. ..... ......... .. 11
Trace Interpolation .. ...... ..... ........ .. ........... .. ..... ............... ...... .. ......... .. .......... .. ..... ...... ..... 14
Velocity Analysis......................................................................................................... 15
Discussion............................................................................................................................. 18
Conclusions .... ........... ...... ...... ..... ..... ....... ..... ...... ............... ... .. .. ...... ................ .......... .. ........... 24
References Cited................................................................................................................... 24

FIGURES
1. Examples of inverse operator and transpose operator........................................................................................... 5
2. Response of inverse operator and transpose operator using parabolic moveout................................ ................. .. 6
3. Example of performance of inverse and transpose operator................................................................................. 7
4. Multiple and primary models of synthetic model filtered with a symmetrical wavelet........................................ 8
5. Models for multiple suppression........................................................................................................................... 9
6. Multiple suppression for three multiple models with and without static shifts..................................................... 10
7. Amplitude analysis of multiple suppression for model 2 ...................................................... ................. ............... 11
8. Multiple suppression for model 3 using an inverse operator ................................................................................ 12
9-10. Multiple suppression for model1 with random-amplitude noise and static shifts using:
9. An inverse operator...................................................................................................................................... 12
10. A transpose operator..................................................................................................................................... 13
11. Straight stacked section of a profile acquired at Lake Baikal ........................... ......... ................... ............... .. ....... 13
12. Multiple suppression after applying static correction .......................................... ~................................................ 14
13. Expanded plot of figure 12 .................................................................................................................................... 15
14. Example of coherency filtering............................................................................................................................. 16
15. Detailed plot for coherency filtering..................................................................................................................... 17
16. Two-dimensional frequency-wavenumber analysis.............................................................................................. 18
17. Example of trace interpolation ... ......... .. ...... ............... ........... .. ............... ...... .. ...... ..... ......... ...... .. ..... ....... ..... ... .. ..... 19
18-19. Examples of velqcity analysis using the:
18. Noise-free model shown in figure 5B .......................................................................................................... 20
19. Noise model with random-amplitude noise and static shifts........................................................................ 21
20. Example of velocity analysis for real data ........... ........... ......................... ............................ ................................. 22
21. Example of multiple suppression with and without static correction.................................................................... 23
22. Effect of random noise and static correction in velocity analysis......................................................................... 24

III
 APPLICATIONS OF VELOCITY-STACK METHODS
TO SEISMIC DATA PROCESSING
By Myung W. Lee and Warren F. Agena

ABSTRACT a family of hyperbolic curves and is described as a velocity

stack using a transpose operator (Thorson, 1984; Thorson
A velocity stack is originally defined as the summation and Claerbout, 1985). Conventional constant velocity analy-
of a seismic event along a hyperbolic moveout trajectory, but sis is one form of velocity stack using a transpose operator.
this definition is extended to include both linear and para- It is a linear transformation and maps the data into velocity
bolic moveouts in this report. Thus, our velocity stack can be space. One problem in using a transpose operator is the
used for modeling events with a variety of moveouts. The smearing effect it has on the data. To better discriminate sig-
ability to model via a velocity stack provides us a number of nal and noise in velocity space, Thorson ( 1984) proposed an
applications in the area of seismic data processing, such as inverse op~rator, which essentially deconvolves the spatial
multiple suppression, coherency filtering, trace interpola- mixing effect of the transpose operator.
tion, and velocity analysis. To take advantage of the compu- Thorson ( 1984) described two approaches in estimating
tational efficiency in the frequency domain, only linear and an inverse operator: the generalized inverse and stochastic
parabolic moveouts are considered here. inverse. The original formulation of Thorson ( 1984) requires
Velocity stacking can be accomplished through the use inversion of a large matrix. Hampson ( 1986) showed that,
of either a transpose operator, such as in conventional veloc- for a typical 96-trace seismic shot of 3-s duration, the matrix
ity analysis, or by an inverse operator. These operators order would be 150,000 by 150,000. Thorson and Claerbout
depend only on the source-receiver offsets of seismic traces (1985) suggested an iterative solution using a conjugate-gra-
and the range of velocities. The choice between the transpose dient method, which is still time consuming but more effi-
and the inverse operator for computational purposes is imma- cient. In order to reduce computational time in estimating an
terial because either operator need only be computed once. inverse operator further, Hampson (1986) proposed a para-
Some advantages in using the inverse operator over the trans- bolic trajectory rather than a hyperbolic trajectory. The par-
pose operator include preserving the input waveform (the abolic trajectory is not a physically observable traveltime
transpose operator tends to bias the frequency content toward curve in the offset domain, but it can be approximated by
lower frequencies) and minimizing any mixing effects. These using NMO corrections with appropriate velocities. This
favorable properties of the inverse operator also allow for approximation enables the use of the Fourier transform, and
better coherency filtering and trace interpolation. a computationally efficient method can be derived.
Because velocity stacking works well at all source- Computing an inverse operator is more time consuming
receiver offsets, its performance in suppressing multiples is than computing a transpose operator. However, because both
better than two-dimensional filtering in the frequency-wave- computations depend only on offset and number of velocities
number domain. When working with good-quality seismic (or slownesses ), these operators need only be computed
data, using an inverse operator for velocity analysis has sig- once, and, hence, the choice between the two methods based
nificant advantages over a transpose operator. However, as on computing time alone is immaterial.
the quality of the data decreases with the addition of more To a first-order approximation, seismic data in the off-
noise, the use of the inverse operator is only marginally better. set domain can be considered as a linear superposition of
many events whose traveltimes describe certain trajecto-
ries. The collection of events are defined as a model in the
INTRODUCTION velocity domain. The inverse operator allows us to accu-
rately estimate such a model from seismic data contami-
Normal moveout (NMO) correction and stacking are nated with random noise, and it provides us with a number
standard, industry, seismic-data-processing sequences. of processing techniques, such as multiple removal, coher-
Combining these two processes is identical to stacking along ency filtering, trace interpolation, and velocity analysis.
2 APPLICATIONS OF VELOCITY-STACK METHODS TO SEISMIC DATA PROCESSING

The performances of both the inverse and transpose opera- a shot can be modeled as an event whose traveltime is repre-
tors for each application are compared, and the advantages sented by a hyperbolic curve. The only difference between
of the inverse operator are discussed. the above two examples is the functional relationship of
Thorson (1984) described m!lltiple suppression by arrival times with respect to the appropriate coordinate sys-
stacking along the primary velocity in velocity space. He uti- tem. In the data space or offset space, we use offset, h, and
lized the property of better velocity discrimination of an two-way time, t, as variables, and in the model space or
inverse operator in his multiple suppression. Another velocity space we use slowness, p, and zero-offset time, 't.
approach to multiple suppression is to model the multiple Therefore, the seismic data d(h,t) can be represented by:
energy by an inverse operator and subtract this multiple
energy from the original data. This approach of multiple sup- d (h, t) = JJu (p, 1:) () (f(t, h, 1:, p)) dpdt + n (h, t)
(1a)
pression by velocity stack is described well by Hampson
(1986), and the actual performance for a field data set is doc- where
umented by Geist and others (1990). This paper employs d(h,t) is the measured seismogram at offset, h, and two-
Hampson's approach and uses models to quantify the perfor- way time, t,
mance of multiple suppression-we emphasize the use of u(p,'t) are seismic models with slowness, p, and zero-off-
static corrections for better multiple elimination. set time, 1:,
The random noise in seismic data cannot be modeled as () is the Dirac delta function,
coherent events using a velocity stack. By transforming the flt,h,'t,p) is the traveltime relationship between data and
data into model space and back to the offset domain, random model spaces, and
noise will be reduced. This property can be used as a coher- n(h,t) is the measurement noise at offset, h, and two-way
ency-filtering or a signal-enhancement technique. Kong and ·time, t.
others (1985) and Monk (1992) described nonlinear, coher- Three functional relationships, flt,h,'t,p), in equation
ency-filtering techniques using the invertible Radon trans- 1a are considered by numerous authors. Stoffa and others
form or 't-p transform. Seismic events can also be predicted (1981) usedflt,h,'t,p) = 1:-(t-ph) as their functional relation-
along various moveouts in model space, and this property ship (this is the well-known slant stack or 1:-p operation),
can be used in trace interpolation. An inverse operator used whereas Thorson (1984) used a hyperbolic relationship,
for coherency filtering and trace interpolation produces data 2 2
flt,h,'t,p) = 1:-J?- p h , for his inverse velocity stack.
with reduced trace mixing (smearing) and preserves the
high-frequency components. Hampson (1986) used a parabolic approximation to a
The compact response, or reduced side lobes, of an hyperbolic traveltime curve (which isflt,h,'t,p) = 't-(t-sh 2 ),
inverse operator in the slowness domain is ideally suited for with a new parameter, s) in his study of inverse velocity
high-resolution velocity analysis (Thorson, 1984 ). However, stack for multiple removal. The meaning of s will be dis-
noise introduced into real data sets degrades the performance cussed later. In this section, it is treated as a parameter
of an inverse operator. This limitation in the presence of similar to slowness, p.
noisy data will be discussed later. Thorson (1984) defined the integral relationship in
equation 1a without a noise term as an operator, L, such that
d(h,t) = Lu(p,'t) and defined a transpose operator, LT, as
u(p,'t) = LTd(h,t). Or:
ACKNOWLEDGMENTS T
u (p, 1:) = L d (h, t) =
We thank John J. Miller and William J. Stephenson
for constructive reviews. All data were processed using
CogniSeis DISCO software with additional programs
ffd (h, t) () (f(t, h, 't, p)) dhdt (lb)

developed by U.S. Geological Survey personnel. In the operator L, the variable 1: is expressed as a function of
t, and in the operator L T, the variable tis expressed as a func-
tion of 1:. Equation 1b indicates that performing normal
THEORY moveout correction and stacking, or stacking along the
moveout defined by flt,h,'t,p), are methods of transforming
Detailed theoretical development of inverse velocity data into velocity space or model space. Under what condi-
stacking can be found in Thorson (1984) and Hampson tions is this operation invertible? In other words, what oper-
(1986). To a first-order approximation, considering only ator will correctly invert the velocity space into the data
arrival times, seismic data can be modeled as a superposition space? As pointed out by Thorson (1984), the transpose oper-
of laterally coherent events. For example, a plane-wave ator defined in equation 1b does not work well for a truncated
response of dipping layers consists of laterally coherent data set, such as any real data set. The purpose of this section
events whose stepout times are linear with distance. Also, a is to derive an efficient inverse operator that performs better
seismic event recorded on a horizontal-receiver spread from than the transpose operator defined in equation 1b.
 THEORY 3

A computationally efficient inverse operator is diffi- Let us take a Fourier transform of equation 2, and let
cult to derive if a hyperbolic moveout is assumed in equa- upper-case variables represent the Fourier-transformed
tion lb. But, for a linear moveout or parabolic moveout, an variables. Then:
efficient method can be derived in the frequency domain
(Hampson, 1986). If we use a parabolic relationship for fin
equation 1a, then:

2
d(h,t) = IIu(s,t)O(t- (t-sh ))dsdt+n(h,t)
I d(h, t) e
-irot
dt =

2
= Iu(s,t-sh )ds+n(h,t) (2)

Equation 2 indicates that the seismic data can be modeled

2
by summing the model along a parabolic trajectory. Equa-
D (h, ro) = I I u (s, t) e-irote -irosh dsdt + N (h, ro)
tion 2 also indicates that the data space can be modeled by
a convolution of the model u(p,t) with the parabolic trajec-
2
tory in the time domain. This convolutional property using -irosh
= JU (s, ro) e ds + N (h, ro) (3)
the parabolic or linear traveltime curve in equation 1 is use-
ful in the frequency domain because convolution in the
time domain can be represented by multiplication in the
frequency domain. Equation (3) is valid for a continuous variable in h and s.
Let us assume that the offset, h, has only discrete values
ht.h2 , .... ,hN, and s has discrete values St.S2, .... ,sM. Then
equation 3 can be written, using a matrix notation (a matrix
is denoted as a bold-faced, upper-case letter) as:

D = LU+N (4)
where
D is a column vector of (D([Link]),D(h2,ro), .... ,D(hN,ro))
with dimension N,
U is a column vector of (U([Link]),U(s2,ro), .... ,U(sM,ro))
with dimension M,
N is a column vectorof(N([Link]),N(h2,ro), .... ,N(hN,ro)), and
2
L is an NxM matrix with element of Lnm =e-irosmh n.

Explicitly, equation 4 can be written:

-iros 1h: -iros 2 h: -irosMh:

D (hp ro) e e ·e U(spro) N (hp ro)
-iros 1 h~ -iros 2 h~ -irosMh~
D (h 2, ro) e e ·e U (s 2, ro) N(h 2, ro)

+ (5)

-iros 1 h~ -iros 2 h~ -irosMh~

D (hN, ro) e e ·e U (sM, ro) N(hN' ro)
4 APPLICATIONS OF VELOCITY-STACK METHODS TO SEISMIC DATA PROCESSING

Now the question is how to compute model space, U, U(sp ro) G 11 G 12 · G 1M D (hp ro)
from the noisy discrete data, D, while minimizing the misfit
U (s2, ro) G2I G22 · G 2M D (h 2 , ro)
function or noise, N. Thorson and Claerbout (1985) demon-
strated that this problem can be solved as an overdeter-
= (9)
mined linear system. The standard least-squares solution to
equation 4 is:

(6) GNIGN2

Equation 6 shows the desirable inverse operator, G = For a given frequency, each element of G is a complicated
(LTL )-1 L~ in the velocity stack. function of h and s, which is easily seen from equation 5,
The meaning of the s in the parabolic moveout is not as but does not involve actual data. In other words, for the
straightforward as the p in the linear or hyperbolic moveout same field geometry (same offset) and range of slownesses,
because this is not a physically observable moveout. The lin- the same inverse operator can be applied to the data set to
ear relationship between the time variables ('t and t) makes compute U or D. Therefore, we can pre-compute an inverse
possible the frequency domain operation in equation 1. operator or a transpose operator before applying it to the
Under which circumstances can a hyperbolic moveout be data. Because we only have to compute the operator once,
approximated by a parabolic moveout that is linear in 't and the fact that computing an inverse operator is more expen-
t? The hyperbolic-moveout equation can be expanded using sive that computing the transpose operator is rather trivial.
a Taylor series: Figure 1A shows the elements of an inverse operator,
J2't + p h Gu, convolved with a zero-phase band-pass wavelet for a
22
= 't J1 + (p 22
h I 't 2 ) parabolic moveout, and figure 1B shows the transpose oper-
ator. The operators were computed with 12 offsets and nine
slownesses. The near-offset and far-offset distances shown
in figure 1 are 0 m and 275m, respectively, and the moveout
time at the far offset ranges from -60 ms to 60 ms. As can be
If (p 2h2!t2 ) << 1, dropping higher order terms can be justi- seen from figure 1, the elements of the transpose operator is
fied. The only adjustable variable that will satisfy this con- essentially a band-pass wavelet located at a trajectory of the
dition is the slowness variable, p, and this adjustment can be parabolic moveout, whereas those of the inverse operator are
accomplished by applying a normal moveout (NMO) cor- a very complicated function of offset and slowness along the
rection with a slowness, Pc• as shown in Hampson (1986). parabolic moveout trajectory. Notice the waveform changes,
Therefore the time (t') after the NMO correction with Pc can including polarity reversal, for the inverse operator. The
be written using equation 7: waveforms of the transpose operator are identical for all off-
sets and slowness values, so the application of a transpose
operator is the same as straight stacking along a parabolic
moveout. The waveforms of the inverse operator are all dif-
ferent from each other, so the application of an inverse oper-
ator is similar to a weighted stack along a parabolic moveout.
If the numerical values of D in equation 9 are identical
Therefore, if Pr is very small or p is close to Pc• a parabolic for all offsets, as for a horizontal event with the same wave-
lets, the model or velocity space can be derived by sum-
approximation can be justified. The parameter s defined in
ming the rows of the inverse operator or summing along the
equation 2 is (p,2/2't), and, strictly speaking, it is not a slow-
offsets. Summing along the offsets is the same operation as
ness parameter. Because the coefficient h2 in the parabolic summing each slowness panel in figure 1. The output of
approximation is a function of zero-offset time, the shape of model space for the operators shown in figure 1 is shown in
a parabola for a given s is a function of zero-offset time, dif- figure 2A for the inverse operator and is shown in figure 2C
ferent from the slant stack. Also, smearing can be expected for the transpose operator. The output of the inverse opera-
as events deviate from the ideal parabolic form. tor is symmetrical in time, but this is not true for the trans-
pose operator. The amplitude of the inverse operator decays
quickly with increasing moveout time, whereas the ampli-
INVERSE AND TRANSPOSE OPERATORS tude of the transpose operator decays slowly with increas-
ing moveout time. The increased lateral spread of the
Denoting the element of inverse operator G as Gij• transpose operator in the slowness domain causes lateral
equation 6 can be written as: smearing in the velocity stack.
 THEORY 5

MOVEOUT TIME AT FAR OFFSET. IN MILLISECONDS

-60 -45 -30 -15 0 15 30 45 60
I I _I .I
I !
I
r
en ~
)
)

~(
0
u
( l(
~ • (

UJ
en
_..
::::::i t •
~
o·
C)
N

l
A

I
I I
I

r •
en
c
z
c
,
u
UJ
en
_..
::::::i • •
~
C)
C)
N

l
B

Figure 1. A, Inverse operator; B, Transpose operator convolved with a zero-phase band-pass wavelet of 4/6-851100 Hz for a para-
bolic moveout with moveout time of -60 ms to 60 ms at the far offset. For each moveout time at the far offset, 12 offsets are shown;
these consist of, from left to right, the near offset (equal to 0 m) to the far offset (equal to 275m).

The difference in performance between the inverse and shows a bias toward lower frequency. This is one reason
the transpose operator can be easily examined by summing why the conventional velocity-analysis method, summing
U over all slownesses. The result is shown in figure 2B for along the assumed velocity trajectory, shows a lack of reso-
the inverse operator and in figure 2D for the transpose oper- lution, as mentioned in Thorson (1984). For horizontal
ator. Equation 3 indicates that integration with respect to events, this result can be represented by:
slowness in the model space with h = 0 will result in the M M N
zero-offset data in the offset domain. Therefore, if accurate
U is modeled from D (data in the offset domain), a summa-
D (h = 0, ro) = L u (sm, ro) = L L Gnm (10)
m=l m=ln=l
tion of the model U along the slowness will reproduce the
input data. The result of this operation for the inverse opera- This equation can be used to examine whether the computed
tor, shown in figure 2B, clearly shows the waveform as iden- inverse operator, G, is correct or not.
tical to the input, which is a zero-phase band-pass wavelet An application of inverse and transpose operators in
with 4/6-851100 Hz, but the result of the transpose operator creating a model space is shown in figure 3. Figure 3A shows
6 APPLICATIONS OF VELOCITY-STACK METHODS TO SEISMIC DATA PROCESSING

MOVEOUT TIME. IN MILLISECONDS · MOVEOUT TIME. IN MllliSECONDS

-60 0 60 -60 0 60
I I !

i
U)

E
0
0

1
A B c D
Figure 2. Response of inverse and transpose operator using parabolic moveout trajectory for a horizontal event. A, Response in slowness
domain of inverse operator; B, Summation of response of inverse operator with respect to slowness; C, Response in slowness domain of
transpose operator; D, Summation of response of transpose operator with respect to slowness. Note difference of the waveform in Band D.

the noise-free input data; figure 3B shows the modeled data MODEL DESCRIPTION
using the inverse operator; and figure 3C shows the result of
applying the transpose operator. As indicated in figure 2, the In order to analyze the performance of the multiple sup-
result of the transpose operator fails to model the higher fre- pression technique by a velocity stack, a 24-channel com-
quency part of the input data. Loss of the high-frequency mon midpoint (CMP) gather, with near offset of 400 m and
component is evident at the far offsets. Results from using an far offset of 2, 700 m, was generated. Figure 4 shows the syn-
inverse operator, however, show an almost perfect reproduc- thetic model filtered with a 4/6-42/54-Hz symmetrical
tion of the input data.
wavelet. Figure 4A shows the multiple model, and figure 4B
shows a primary model. In figure 4A, the event near 2 s at the
near offset represents a reference trace to measure the perfor-
APPLICATIONS mance of the multiple suppression technique, and events
The ability to extract model space from a noisy data set below 2.5 s represent multiple energy. The average root-
can be used for a variety of seismic-data-processing tech- mean-square (RMS) amplitude of multiple energy is very
niques, and a number of applications of inverse velocity close to the RMS amplitude of the reference trace. The RMS
stack, such as velocity analysis, is discussed in Thorson velocity of multiples is modeled by a linearly increasing
(1984). Hampson (1986) presented an application to remove function with depth: 1,400 m/s at 2,500 ms and 1,600 m/s at
multiples using a parabolic approximation. In this paper, 4,500 ms. Four primary events with identical amplitude and
applications of the inverse and transpose operators will be polarity are chosen and shown in figure 4B. The amplitude
discussed in the area of multiple suppression, coherency fil- of the primaries in figure 4B is 20 percent of the reference
tering, velocity analysis, and spatial interpolation. amplitude in figure 4A, and the locations and RMS velocities
of four primaries in figure 4B are: 3,000 ms with 1,600 m/s;
3,420 ms with 1,800 m/s; 3,800 ms with 2,000 m/s; and
MULTIPLE SUPPRESSION
3,950 ms with 2,200 m/s.
The methodology of multiple suppression using the In our model study, figure 4A was mixed with figure 4B
velocity stack can be summarized in the following five steps: for a variety of scale factors. Three synthetic models were
1. Apply NMO correction using stacking velocities of generated for the analysis.
multiples. Modell: Multiple strength 14 dB higher than primaries.
2. Model the multiple event using the inverse or trans-
pose operator. Model 2: Multiple strength almost equal to primaries.
3. Subtract the modeled multiple events in step 2 from Model3: Multiple strength 20 dB higher than primaries.
the original data. Figure 5A shows model 1, and figure 5B shows model
4. Remove the NMO correction applied in step 1. 2. Notice that, in modell, primary events cannot be observed
5. Stack with the primary velocity. by the eye. The amplitude relationship in figure 5A is similar
The performance of multiple suppression outlined above is to the field data example presented later. Figure 5C is model
analyzed using models and a real data set. 1 with random noise added and will be used later.
 APPLICATIONS 7

OFFSET. IN METERS
400 1,600 2,600 400 1,600 2.600 400 1,600 2,600

5.00

A B c
Figure 3. Example of performance of inverse and transpose operator. A, Original input; B, Modeling by an inverse operator; C, Modeling
by a transpose operator.

NOISE-FREE MODEL primaries. Also, the result shown for model2, "a" part, indi-
cates that straight stacking will be sufficient to suppress mul-
The results of straight stacking for noise-free models tiples when the multiple strength is almost the same as that
are shown in the "a" parts of figure 6, and results of multiple of the primaries (model 2). However, it demonstrates that, if
suppression by a velocity stack using the inverse operator the multiple strength is more than 14 dB higher than the pri-
(inverse velocity stack) are shown in the "c" parts of figure maries, straight stacking is not enough to enhance the prima-
6. The results for the noise-free model (the "c" parts of fig. ries (fig. 6, models 1 and 3, "a" parts), and some kind of pre-
6) indicate that multiple suppression by the inverse velocity stack multiple-suppression technique, such as an inverse
stack suppresses the multiple energy enough to recognize the velocity stack, is required.
8 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING

OFFSET. IN METERS Gaussian random noise were introduced into the CMP
400 1,000 1,600 2,200 400 1,000 1,600 2,200 gather. The static shift was chosen for our noise model for
the following reasons:
1. The velocity stack or other multiple-suppression tech-
nique keys on the alignment of signal and coherent
noise, so the performance depends on the degree of
misalignment.
2. Other kinds of misalignment, such as waveform
changes due to different ghosting, can be modeled
approximately by a static shift.
3. The effect of inaccurate multiple velocities can be
(/)
treated as a static problem.
0
:z 4. Static shifts account for a variation of streamer depth
0
(..')
during data acquisition .
.~ 3.00
~ The performance of velocity stack with static shifts is
u.i
~ also shown in figure 6. The "b" parts of figure 6 represent
i=
the results of straight stacking, and the "d" parts of figure 6
represent the results of multiple suppression by the inverse
velocity stack. Notice that static shifts on the order of 4 ms
does not significantly affect the performance of straight
4.00
stacking (compare the "a" parts with the "b" parts of fig. 6).
4
However, the performance of the inverse velocity stack is
highly dependent on the static shift. In the "d" parts of fig-
ure 6, the primaries are still recognizable for model 1, but
not for model 3. The performance of the velocity stack
under the influence of static shifts is shown in figure 7B.
5.00 ---t+t++t+tt+H-H-++t++t-+t+t---i-++t++++++-H++-++1-++-<i-++-1--l--
The strength of primary events in figure 7A (noise free) and
A B figure 7B (static noise) are similar, but the amount of multi-
ple suppression under the static shifts is on the order of 20
Figure 4. Synthetic model filtered with a 4/6--42/54-Hz sym- dB (it is about 40 dB in the noise-free case). In other words,
metrical wavelet. A, Multiple model except the event at 2 s, which
if the multiple energy contaminated by static shifts is more
serves as a reference amplitude; B, Primary model with four
than approximately 20 dB higher than that of primaries, it
events. See text for additional explanation.
is difficult to enhance primaries over multiples using the
inverse velocity stack.
Notice that the strength of multiple energy is almost the
In order to see the details of multiple suppression, an
amplitude analysis of the result of velocity stack for model 2 same as that of the primaries in figure 6, model 3, "d" part,
(fig. 6, model 2, "a" part) was performed, and the result is after multiple suppression. Recalling that the multiple
shown in figure 7A. This figure indicates that the velocity strength before the multiple suppression for model 3 is 20 dB
stack suppressed the multiple energy on the order of 40 dB. higher than the primary, the signal-to-noise ratio for model3
Notice the five peaks near the -10-dB line in figure 7A, indi- can be considered as a threshold for recognizing the prima-
cating almost perfect primary recovery. On the other hand, a ries. However, in order to differentiate primaries from mul-
straight stack only suppress the multiple energy on the order tiples, the amplitudes of primary events after multiple
of 12 dB. Approximately a 28-dB enhancement is achieved suppression should be a couple of decibels higher than those
by the velocity stack. of the multiples. Figure 6, modell, "d" part, whose multiple
The above analysis is for the ideal case in which we strength is about 14 dB higher than the primaries, clearly
know exact multiple and primary velocities. For this ideal shows the primaries after multiple suppression even though
case, there is neither random nor other coherent noise, and there are 4-ms static shifts. Therefore, in processing real data
each channel has an identical waveform. These ideal condi- sets, we could not enhance the primaries if the strength of
tions are never met with real data. multiples is about 17 dB greater than that of the primaries
(see fig. 6, model1, "d" part and fig. 6, model3, "d"·part).
NOISE MODEL Approximately a 20-dB loss of performance of the
velocity stack for the noise model is due to small static shifts.
In order to evaluate the performance of the velocity In theory, static shifts can be computed and static corrections
stack in a "real" data situation, noise in the form of random can be applied before the velocity stack. One of the results of
static shifts, in the range of 4 ms, and random amplitude with applying static shifts computed by non-surface-consistent
 APPLICATIONS 9

OFFSET. IN METERS
400 1.600 2.600 400 [Link] 2.600 400 1,600 2,600

5.00

A B c
Figure 5. Models for multiple suppression. A, Multiple modell: The multiple strength is about 14 dB higher than primaries; B, Multiple
model2: The multiple strength is almost the same as primaries; C, Additive random noise to model!. The ratio of primary signal to random
noise (not coherent multiple noise) is about 1.

static, or "trim static" in industry terminology, is shown in result of straight stacking with random noise, and figure 9B
figure 8 for model3. As expected, the result (fig. 8D) is com- shows the result of straight stacking with both random
parable to the result for the noise-free model (fig. 8B). noise and static shifts in the range of 4 ms. There is not
The effect of random-amplitude noise in multiple sup- much difference between figures 9A and 9B. Both are too
pression is analyzed using model 1 by adding Gaussian ran- noisy to recognize primary events. The result of applying
dom noise with maximum amplitudes to those of primaries an inverse operator to the random noise section without
(fig. 5C). The result of multiple suppression using an static shifts is shown in figure 9C. All primary events are
inverse operator is shown in figure 9. Figure 9A shows the clearly observable. The effect of both random-amplitude
10 APPLICATIONS OF VELOCITY-STACK METHODS TO SEISMIC DATA PROCESSING

1111111 IIIII

.DE.
1111111 IIIII 1111111 111111
MODEL2
MODE~

2.00 --:11iiilililiP-~iiliiliiP----~-P-~iiiiiiP-

(( ((
. ;..

a b c
Ill
d a I~, c
Ill 1~, ~: 9,

Figure 6. Multiple suppression for three multiple models with and without static shifts in the range of 4 ms using an inverse operator in
the velocity stack: a, Straight stack without static shifts. b, Straight stack with static shifts. c, Multiple suppression without static shift. d,
Multiple suppression with static shifts.

noise and static shifts is shown in figure 9D, and the result source. Unstacked data consists of 24-fold CMP gathers
of application of static correction is shown in figure 9E. with a near offset distance of 400 m and a far offset dis-
Comparing with figure 6 for model 1, it is obvious that tance of 2, 775 m.
noise in the form of 4-ms static shifts has a more pro- Our model study indicates that correcting for static
nounced effect than an equal amount of random noise. shifts is a very important step in suppressing multiple
The result of multiple suppression using the transpose energy. Thus, one additional step was added to the steps out-
operator is shown in figure 10. In comparing figures 9C and lined at the outset of the multiple suppression section: static
9E with figures IOC and IOE, the advantage of the inverse computation and application after step 1. The static compu-
operator over the transpose operator is evident. But when tation was performed in the CMP domain.
there is too much noise, the results shown in figure 9D and Figure 12 shows the result of multiple suppression. Fig-
figure 1OD imply that there is not much difference between ure 12A is the result of applying two-dimensional velocity
the inverse and transpose operators. Hence, for high-quality filtering, or F-K filtering, a technique similar to the one
input data, the inverse operator works better than the trans- described by Ryu (1982), and figure 12B shows the result of
pose operator. But, as the amount of noise (particularly in the the application of the inverse velocity stack. For the inverse
form of static time shifts) increases, the advantage of the velocity stack, the multiples are modeled with seven slow-
inverse operator over th~ transpose operator decreases. nesses, and the range of parabolic moveouts is -20 ms to 40
ms at the far offset. Both techniques performed adequately
REAL DATA EXAMPLE for the multiple suppression below about 4,600 ms. Overall
signal-to-noise ratio, however, is better in the section pro-
Figure 11 shows an original stacked section generated duced by the application of inverse velocity stack. The
using a standard processing sequence. This data set was details of figure 12 near CMP 2,500 is shown in figure 13,
acquired in Lake Baikal using a tuned air-gun array as a and this figure clearly demonstrates that the inverse velocity
 APPLICATIONS 11

0 0~----------------------------------~

-10 - -10 f.-

T
-20 - Tz
-20 f.-

.
~
[Link]
C/) C/)
-LU
Ocr::
Zo...
0...
~
C/)

~
C/)
.....J 0
.~
LU
CD Ci5
u

-~~~
LU
0
C/)
LU
0::
a:
... ..
-
~··
~ -30 -30 -

·.·~
~
uJ C/)
LU
~
0 C/)
~
1- (5
:::::i z
0...
::2 .]
<(

·.~ .. ~ ·.

~
-40 - -40-

~·
-50 f.-
~: ~ u... · ~-
.v .
-50 -
~ ~
~· ~ ·. A
-60 I

2
J I •

3
I

4 5
-60--------~
I

2
..
L-----L--------L-----L~
3 4
I

5
I ~
TIME. IN SECONDS
A B

Figure 7. Amplitude analysis of multiple suppression for model2 using an inverse operator in the velocity stack. Each dot represents a
root-mean-square (RMS) amplitude within 16ms, and solid lines are running medians ofthree amplitudes. A, Without static shift; B, With
static shift in the range of 4 ms.

stack outperforms the F-K operation in multiple suppression. COHERENCY FILTERING

Notice much stronger and clearer dipping events are imaged
in figure 13B. It is well known that stacking suppresses random noise
We tested many seismic profiles using the inverse proportional to the root of the number of traces in the stack.
velocity stack versus the F-K domain approach. In almost If all reflections are horizontal, this is a good method for
every case, the inverse velocity stack worked better than enhancing signal-to-noise ratio (this is the fundamental con-
the F-K domain approach, particularly for data sets with cept applied in CMP stacking). In the stacked section, this is
high noise and low frequency. A similar conclusion was not an appropriate method because there are numerous non-
reached by Geist and others (1990) in their study of long- horizontal events. In order to enhance dipping events, stack-
period (> 1 s), first-order, multiple removal over a conti- ing should be done along different moveouts. Stacking along
nental slope. Better performance of the inverse velocity the linear moveouts can be effectively implemented using a
stack for the lower frequency range has an important impli- slant stack in the 't-p domain. This kind of coherency filter-
cation for deep crustal studies. As mentioned in Hampson ing is necessary for data with a low signal-to-noise ratio.
(1986), the better performance of the inverse velocity stack In this section we will discuss coherency filtering in the
is due to the fact that inverse operators work well for all context of model extraction from noisy data. As indicated in
offsets, whereas the F-K approach leaves significant multi- equation 6, the model space can be estimated by applying the
ple energy on the near-offset traces. inverse operator to the offset-domain data set. This estimated
12 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING

1.00 1.00

~ ~
2.00 "i'((lft'ft
~

J
2.00
\' ( ~
~ '(

it
(

en
0
~
~
z ~ <
0
u
~ 3.00
~
w
,, ~
~
C/) { { (
Cl
:z I ;.
0
.~
(..)

~ 3.00
~ I
u.i
~
i=
~\·
·r; . .
P'
!'I
4.00 "'
~

'1 ~~~ ~

) •)•)

t
4.00

fr
5.00
!
H
'i
A B c D E
) ) ) )

Figure 9. Multiple suppression for model 1 with random-

amplitude noise and static shifts using an inverse operator in the
velocity stack. Primary signal-to-random-noise ratio is about 1,
and the range of static shifts is ±4 ms. A, Straight stack with ran-
dom noise only; B, Straight stack with random noise and static
shift; C, Multiple suppression with random noise only; D, Multi-
5.00
ple suppression with random noise and static shift; E, Multiple
suppression with random noise and static shift after static correc-
tion with "trim static" (see text) in the CMP domain.
A B c D

Figure 8. Multiple suppression for model 3 using an inverse op-

erator in the velocity stack. A, Straight stack without static shift;
B, Multiple suppression without static shift; C, Multiple suppres-
sion with static shift in the range of 4 ms; D, Multiple suppression
with static shift after correction with "trim static" (see text) in the
CMPdomain.
 APPLICATIONS 13

1.00 COMMON MIDPOINT

1,375 1,975 2,575
2.00

2.00 ~
~
~

----- en
0
z
4.00

~~~ ~~» 0
(._)
5.00
~ t UJ
~ en
en
~ ~
a ~ ~/
z
0
f < u.i
~
(._)
i=
~ 3.00
~ 6.00
u.i
~
i=
{ ~· <'(, :t (
.
('\W~i· •}'
~
4.00 "'

tUi ~~~1~} Figure 11. Straight stacked section of a profile acquired at Lake
Baikal in the fall of 1992. Notice strong, first-order, water-bottom
multiple and pegleg multiples after 4,600 ms. See text for acquisi-
tion parameters.

5.00

A B c D E

Figure 10. Multiple suppression for model 1 with random-

amplitude noise and static shifts using a transpose operator in the
velocity stack. Primary signal-to-random-noise ratio is about 1,
and the range of static shifts is ±4 ms. A, Straight stack with ran-
dom noise only; B, Straight stack with random noise and static
shift; C, Multiple suppression with random noise only; D, Multi-
ple suppression with random noise and static shift; E, Multiple
suppression with random noise and static shift after static correc-
tion with "trim static" (see text) in the CMP domain.
14 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING

COMMON MIDPOINT
1,375 1,975 2,575 1,375 1,975 2,575
2.00

4.00

Cl)
c
:z
0
u
U.J
Cl) 5.00 ·.··
~
u.i
~
i=

6.00

A B

Figure 12. Multiple suppression after applying static correction to data shown in figure 11. A, Multiple suppression using two-
dimensional frequency-wavenumber domain filtering; B, Multiple suppression using an inverse operator in the velocity stack.

model or signal can be inverse-transformed back to the offset offset of the center trace is 0 m). The output is estimated at
domain data set using equation 4 without the noise term. The the center trace. Both figures 14A and 14B show a reduction
presence of the noise in equation 4 indicates that the trans- of background random noise or the enhancement of local
form generated by equation 6 is not a complete space. Some coherency estimated within five traces. Details of coherency
of the noise, n(h,t), that has not been included in U is lost filtering in the time domain are shown in figure 15, and fig-
when transforming back to the offset domain. This operation ure 16 shows the frequency-wavenumber domain represen-
can be written as: tation of the data shown in figure 15. The low-frequency
enhancement and mixing effects produced by the transpose
(11) operator are apparent in figures 15 and 16.

where
fi is the coherency-filtered output. TRACE INTERPOLATION
Figure 14A shows the result of applying an inverse
operator as a coherency filter to the section shown figure The concept of coherency filtering by inverse velocity
12B, and figure 14B shows the result of applying a transpose stack can be applied to trace interpolation. The only differ-
operator. The inverse or transpose operator was designed ence is the application of the L operator in equation 11. In
using linear moveout with 11 slowness steps and five offsets. equation 11, the element of L applying to model space U to
The range of linear moveout at the far offset is ±24 ms (the transform back to the offset domain is the same element used
 APPLICATIONS 15

COMMON MIDPOINT
2,275 2.475 2,675 2.275 2.475 2,675
5.00

en
Cl
z
0
(..)

~ 6.00 -----i~'\9-J
~
uJ
~
i=

7.00

A B

Figure 13. Expanded plot of figure 12 near CMP 2,500. A, Multiple suppression using two-dimensional frequency-wavenumber
domain filtering; B, Multiple suppression using an inverse operator in the velocity stack.

for computing G. But in trace interpolation, a different L can VELOCITY ANALYSIS

be used for the output. In other words, a different offset
value, not included in the input data, can be used for the Moveouts in the CMP domain depend on the interval
interpolation. Assume that a trace is interpolated using two velocities of the layered Earth, and the model space is esti-
input traces, and the interpolation is done at the middle of mated in the slowness domain. Thus, the inverse velocity
two traces using linear moveout. Then the L used for com- stack is a logical choice for velocity analysis. When using
puting G can be written as: hyperbolic moveouts, the slowness can be directly related to
the Earth's velocity. But, in using parabolic moveouts, it
only provides a pseudo or apparent velocity after NMO cor-
(12)
rection with velocity, Vc (or slowness, Pc-see eq. 8).
Parameter "sh 2" in parabolic moveouts is presented as
The element L* for the interpolation can be written as: the moveout time at the far-offset distance; thus, a relation-
*)
L* =
l e
-iws 1h 1* -iws 2 h 1
e (13)
ship is required to derive an actual velocity from the velocity
stack using parabolic moveouts. The relationship between
where apparent and true velocity can be derived from equation 8:
h* =(hi+ h2)/2 1 1
-+-
An example of trace interpolation using an inverse 2 2 2
operator is given in figure 17. The trace interval of figure 17
v vr vc
is 25 m, and five traces were used for the interpolation. An Using the parameters:
inverse operator was designed to pass seismic events whose
linear moveouts are within ±4 ms per trace. Highly dipping 2
v (14)
events at the left side and random noises are reduced in the
interpolated section (fig. 17 B). Also, trace mixing, which is
obvious when using the transpose operator, is not apparent in Because the parameter is given by the moveout time at the
the interpolated section. far offset, equation 14 can be written as:
16 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING

COMMON MIDPOINT
1,375 1,975 2,575 1,375 1,975 2,575
2.00

4.00

en
Cl
z
0
u
UJ
en 5.00
~
u.i
~
t=

6.00

A B
Figure 14. Example of coherency filtering of the section shown in figure 12B. Coherency filter was designed using five offsets with
the moveout time at the far offset of ±24 ms. The offset at the center trace is 0 m, and the trace interval is 50 m. A, Inverse operator; B,
Transpose operator.

shows linearly increasing velocity for the multiples and four

primary events. The moveout time of P4 is -168 ms at two-
v = (15)
way time of 3,950 ms. Putting this number into equation 15,
the correct velocity by a parabolic moveout is 2,197 m/s,
which is very close to the input velocity of 2,200 m/s. How-
ever, the output from the transpose operator or conventional
where
velocity analysis, shown in figure 18B, show a smearing
Mt is moveout time at the far offset, and effect of the adjacent moveouts. The above analysis was
hmax is the far-offset distance. applied to the noisy model, random noise (primary signal-to-
As indicated in equation 15, the true velocity depends on the random-noise ratio of about 2), and static shifts (in the range
zero-offset time as well as the correction velocity for NMO of 4 ms), and the result is shown in figure 19. Although the
and moveout time. result from the inverse operator is better than the result from
A velocity analysis for the noise-free model shown in the transpose operator, the difference is not as dramatic as
figure 5B is given in figure 18 with linearly increasing veloc- that shown in figure 18.
ity ofvc = 1,500 m/s at 2,500 ms and Vc = 1,700 m/s at 4,500 An application to a real data set is shown in figure 20.
ms. The range of moveout times at the far offset is -240 ms The result using the application of an inverse operator (fig.
to 120 ms, with a 12-ms interval. The output of the applica- 20A) produces better resolution but is noisier than the result
tion of the inverse operator, shown in figure 18A, clearly produced by the transpose operator (fig. 20B)
 COMMON MIDPOINT
2.275 2.475 2.675 2.275 2.475 2,675 2.275 2.475 2.675
5.00

6.00
en
0
z ;J>
0
u
UJ
en
::g
t""'
-zL..iJ .......
n
;J>
::E >-3
i= .......
0
z
en

7.00

A B c
Figure 15. Detailed plot for coherency filtering near CMP 2,500 shown in figure 14. A, Original data without coherency filtering; B, Coherency filtering by an inverse operator;
C, Coherency filtering by a transpose operator.

"""'
......,J
18 APPLICATIONS OF VELOCITY-STACK METHODS TO SEISMIC DATA PROCESSING

20.00

15.00
q

N
~·
li:
U.J
:I:
z
~~I
;:. 10.00
(_)
z I ~~ .
~~ I

~~ ~
U.J

~·
:::::>
0
U.J
a:
u..
I
,, IIII,
5.00
Itt,

0.00

-0.5 0.0 0.5 -0.5 0.0 0.5 -0.5 0.0 0.5

WAVENUMBER
A B c
Figure 16. Two-dimensional frequency-wavenumber analysis for the data shown in figure 15. A, Original data without coherency filtering;
B, Coherency filtering by an inverse operator; C, Coherency filtering by a transpose operator.

DISCUSSION a computer program for equation 6, (LTL+A)- 1 is used for

the computation of the inverse operator. All examples in this
Many conventional seismic processing methods are paper are generated by adding 1 percent of each diagonal ele-
based on the transpose operator or conjugate operator ment to the diagonal terms. The amount of noise can be
(Claerbout, 1992). Although the transpose operator applied adjusted to tune the performance of the inverse operator.
to perfect data generally gives an imperfect result, some- Thorson and Claerbout (1985) defined this form of inverse
times the transpose operator performs better than the as "stochastic inverse" (Aki and Richards, 1980) and suggest
inverse operator because the conjugate operator tolerates using the ratio of noise variance to variance of the points in
imperfections of input data owing to the limited aperture, the model space for A. The method proposed by Thorson
truncation, nonstationary noise, missing data, etc. On the and Claerbout (1985) may provide an optimum inverse oper-
other hand, an inverse operator gives a perfect result for a ator for each data set, but the operator cannot be pre-com-
perfect data set but has the possibility of noise amplifica- puted because A cannot be determined without data.
tion when using actual data and requires more information Two basic assumptions in inverse velocity stacking are
to do a better job. As can be seen from equation 6, the that, along various moveouts, events are aligned and have
transpose operator is the first step of inversion, and, if uniform amplitudes. Multiple suppression using real data
(LTL)- 1 is close to an identity matrix, the inverse operator sets shows improved results after static time corrections are
can be approximated by a transpose operator. applied to improve alignment. However, the improvements
To reduce noise amplification of the inverse operator were marginal in most cases. One example is shown in figure
applied to an imperfect data set or to stabilize the inverse 21. This suggests that simple static corrections alone will not
operator, a random noise factor, A, such as the one used for suffice in certain cases. This is because our waveform
spiking deconvolution, can be applied. When implementing changes in a CMP gather. This waveform change violates the
 DISCUSSION 19

COMMON MIDPOINT
1.250 1,300 1,350 1.400 1,250 1.300 1,350 1.400

A B

Figure 17. Example of trace interpolation using an inverse operator. Trace interval is 25 m, and five traces are used for the interpolation.
A, Input data; B, Interpolated data.

basic assumptions of inverse velocity stacking. Thus, pre- others (1985) proposed a nonlinear signal estimator by mul-
processing, such as amplitude balancing and waveform tiplying semblance (Taner and Koehler, 1969) to the slant-
deconvolution, will be beneficial to the performance of the stacked data and inverse-transforming the weighted slant-
inverse operator. stack data. Monk (1992) presented a diversity slant stack,
When inverse velocity stacking is applied for irregular which is a weighted slant stack based on the power in the
offset gathers similar to land data or in a marine data set, slant-stack domain. Both Kong and others (1985) and Monk
whose CMP fold is less than the number of channels, it is (1992) applied a transpose operator in equation 11. In order
impractical to compute different inverse operators for each to preserve the original frequency content, which is altered
offset gather. For example, 24-fold marine data acquired by applying a transpose operator, Monk (1992) and Kong
with 96 channels have 4 sets of offsets in the CMP domain, and others (1985) applied an appropriate Rho filter (Claer-
and 4 sets of inverse operators are required. To handle irreg- bout, 1984) to the output. The main purpose of weighting,
ular geometries, one effective method is to apply differential proposed by Monk (1992), is to reduce the smearing effect
NMO correction to a reference offset gather. of the transpose operator and to modify the transpose opera-
The estimation of model space from the noisy input tor, more or less, to look similar to an inverse operator. A
data by means of the velocity stack is defined as a process of weighted slant stack using an inverse operator provides an
coherency filtering. Coherency filtering improves the signal- alternative to the conventional approach to signal enhance-
to-noise ratio, as shown in figure 15, and helps subsequent ment. The price we have to pay for a weighted slant stack for
processing steps, such as migration. Most signal-enhance- coherency filtering in this formulation is that U should be
ment techniques use weighted velocity stacks. Kong and converted in the 't-p domain to compute and apply weight.
20 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING

MOVEOUT TIME, IN MILLISECONDS

-240 120 -240 120

.. /11
2.00 ~--+-+-1--+-4-~~-f--4-.i~~~~~,.......,,,lliiiillll,.r-- --+-+++-+~1-+++-+-++H~~r'w'irrn•-'E-f'r'!Jlnr~r'rlr,~
~cP> )-, ~f

en
Cl
z
0
u
UJ
en
~
ui
~
i=

MULTIPLE

A B
Figure 18. Example of velocity analysis using the noise-free model shown in figure 5B. The range of moveout time at the far offset is
between -240 ms and 120 ms, with a 12-ms interval. The correction velocities for NMO to approximate parabolic moveouts are linearly
increasing velocities with 1,500 rnls at 2,500 ms to 1,700 rnls at 4,500 ms. A, Inverse operator; B, Transpose operator.

The interpolation method presented here by means of method has the ability to keep high-frequency components
the velocity stack implicitly uses a local coherency and dip (as opposed to other interpolation methods based on the
estimation (or slowness in this paper) and is subject to the transpose operator).
aliasing problem when interpolation distance is too large. One problem for the velocity analysis of real data using
The formulation of interpolation is cast as a prediction in an inverse operator is low signal-to-noise ratio compared to
the offset domain in this paper, but the essential operation the output by a transpose operator. One way to overcome this
is similar to interpolating the amplitude along the direction difficulty is to adjust the random noise component when
of lateral coherence of each event (Barden, 1987; Lamer deriving the inverse operator. More random noise in the
and others, 1981). This approach differs from Spitz (1991), design of the inverse operator increases the signal-to-noise
who derived a method based on the fact that linear events ratio but also increases lateral smearing. An optimum trade-
may be interpolated exactly, regardless of spatial interval, off between resolution and smearing can be found by trial
without any information regarding their true dip. He inter- and error. Figure 22 shows a velocity analysis with an
polated traces using forward-backward, one-step prediction inverse operator using 2 percent noise for A, instead of 1 per-
filters in the frequency domain. However, the proposed cent used for figure 20A with static corrections. The output
 DISCUSSION 21

MOVEOUT TIME. IN MILLISECONDS

-240 120 -240 120

< i
I ) } ~ ~
< < L< ~ ~ r1 J 1

''~(
I ••

:
rr'''''''f
((~~~~({ f
~

(I)
c
z
8
LU
Cf)

~
u.i
~
i=

A B

Figure 19. Example of velocity analysis for a noise model. The noise model was generated by adding random-amplitude noise and
static shifts to the noise-free model shown in figure 5B. Primary signal-to-random-noise ratio is about 2, and the range of static shifts is
±4 ms. The range of moveout time at the far offset is between -240 ms and 120 ms, with a 12-ms interval. The correction velocities for
NMO to approximate parabolic moveouts are linearly increasing velocities with 1,500 rn!s at 2,500 ms and 1,700 rn!s at 4,500 ms. A,
Inverse operator; B, Transpose operator.
22 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING

MOVEOUT TIME. IN MILLISECONDS

-100 -50 0 50 100
2.00

2.50

en
c
z
0
(_)

~ 3.00
~
w
~
i=

3.50

4.00
A
2.00

2.50

en
c
z
0
(_)

~ 3.00
~
w
~
i=

3.50

4.00
B

Figure 20. Example of velocity analysis for real data having 24 folds and near-offset distance of 400 m and far-offset distance of 2,777 m.
The constant correction velocity of 1,500 m/s is used for NMO to approximate parabolic moveouts, and the range of moveout time at the far
offset is between -100 ms and 100 ms, with a 10-ms interval. A, Inverse operator; B, Transpose operator.

Figure 21 (facing page). Example of multiple suppression with

and without static correction using an inverse operator. A, Straight
stack; B, Multiple suppression without static correction; C, Multiple
suppression with static correction.
 DISCUSSION 23

g
M

1-
z
00...
0
~g
zM
0
~
~
0
(..)

0
0
M

0 0
0 0
N Lri
SONOJ3S Nl '3V\IIl
24 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING

MOVEOUT TIME. IN MILLISECONDS

-100 -50 0 50 100
2.00

2.50

en
c
z
0
<...>
~ 3.00
u..i
~
i=

3.50

4.00

Figure 22. Plot showing the effect of random noise and static correction in a velocity analysis using the inverse operator. This figure is the
same as figure 20A except that ( 1) the amount of random noise in the inverse operator is 2 percent instead of 1 percent, as used in figure 20A,
and (2) static correction was applied to the input data.

of the transpose operator with static corrections is not much 2. Because velocity stack relies on the alignment of seis-
different from the result shown in figure 20B. However, the mic events along the moveouts, static correction of
result shown in figure 22 demonstrates that the improved the CMP gathers is important in multiple suppression.
performance of the inverse operator can be obtained with Model studies indicate that about 4-ms static shifts in
static corrections and by tuning the inverse operator with a gather degrades multiple suppression about 20 dB.
varied random noise. 3. Coherency filtering or trace interpolation by an
inverse operator can be considered as a high-resolu-
tion method and minimizes trace mixing inherent to
CONCLUSIONS output by the transpose operator.
4. Velocity analysis by an inverse operator shows high
A number of applications of velocity stack, such as discrimination ability of close velocities, but the
multiple suppression, coherency filtering, trace interpola- amount of noise, such as static shifts, controls the per-
tion, and velocity analysis, are described in this paper. The formance rather dramatically. To obtain high-resolu-
advantage of using an inverse operator over a transpose tion velocity analysis for field data, pre-processing,
operator is in its ability to accurately estimate the model such as static correction and amplitude adjustment,
space from noisy seismic data. A transpose operator has a should be applied.
bias toward lower frequencies and more spread in the slow-
ness domain. These properties result in trace mixing or
smearing in the velocity stack. However, the advantage of REFERENCES CITED
using an inverse operator over the transpose operator
decreases as noise increases. Based on this study, the follow- Aki, K., and Richards, P.G., 1980, Quantitative Seismology, Theory
ing conclusions can be derived. and Methods: San Francisco, W.H. Freeman Company, 911 p.
1. Because of the effectiveness of velocity stack at all Barden, V., 1987, Trace interpolation in seismic data processing:
offsets using an inverse operator, multiple suppres- Geophysical Prospecting, v. 35, p. 343-358.
sion by velocity stack is better than the frequency- Claerbout, J.F., 1984, Imaging the Earth's Interior: Blackwell Sci-
wavenumber domain approach. entific Publications, 398 p.
 REFERENCES CITED 25

---1992, Earth Sounding Analysis, Processing Versus Inver- Ryu, J.V., 1982, Decomposition (DECOM) approach applied to
sion: Blackwell Scientific Publications, 304 p. wave field analysis with seismic reflection records: Geophys-
Geist, E.L., Hart, P.E., and Fisher, M.A., 1990, A comparison of ics, v. 47, p. 869-883.
multiple suppression methods across a continental slope: Eos,
Spitz, S., 1991, Seismic trace interpolation in the F-X domain: Geo-
Transactions, American Geophysical Union, v. 71, p. 1589.
physics, v. 56, p. 785-794.
Hampson, Dan, 1986, Inverse velocity stacking for multiple elimi-
nation: Journal of the Canadian Society of Exploration Geo- Stoffa, P.L., Diebold, B.P., and Wenzel, F., 1981, Direct mapping
physicists, v. 22, p. 44-55. of seismic data to the domain of intercept time and ray param-
Kong, S.M., Phinney, R.A., and Roy-Chowdhury K., 1985, A non- eter-A plane wave decomposition: Geophysics, 46, p.
linear signal detector for enhancement of noisy seismic record 255-267.
sections: Geophysics, v. 50, p, 539-550. Taner, M.T., and Koehler, F., 1969, Velocity spectra-Digital com-
Lamer, K., Gibson, B., and Rothman, D., 1981, Trace interpola- puter derivation and applications of velocity function: Geo-
tion and the design of seismic surveys [abs.]: Geophysics, v. physics, v. 34, 859-881.
46,p.407. '
Monk, D.J., 1992, Approach to optimum slant stack and its applica- Thorson, J.R., 1984, Velocity stack and slant stack inversion meth-
tion as a seismic noise attenuator [expanded abs.]: Society of ods: Palo Alto, Calif., Stanford University, Ph.D. thesis.
Exploration Geophysicists, 62nd Annual International Meeting Thorson, J.R., and Claerbout, J.F., 1985, Velocity stack and slant
and Exposition, p. 1161-1164. stack stochastic inversion: Geophysics, v. 50. p. 2727-2741.

Published in the Central Region, Denver, Colorado

Manuscript approved for publication May 26, 1994
Edited by Richard W. Scott, Jr.
Graphics prepared by Gayle Dumonceaux
Photocompositon by Carol Quesenberry

-t..U.S. ~ PRIN'l'1N:; OFFICE: 1994-673-046/86065

 SELECTED SERIES OF U.S. GEOLOGICAL SURVEY PUBLICATIONS

Periodicals Coal Investigations Maps are geologic maps on topographic or

Earthquakes & Volcanoes (issued bimonthly). planimetric bases at various scales showing bedrock or surficial geology,
Preliminary Determination of Epicenters (issued monthly). stratigraphy, and structural relations in certain coal-resource areas.
Oil and Gas Investigations Charts show stratigraphic informa-
Technical Books and Reports tion for certain oil and gas fields and other areas having petroleum po-
tential.
Professional Papers are mainly comprehensive scientific reports
Miscellaneous Field Studies Maps are multicolor or black-and-
of wide and lasting interest and importance to professional scientists and
white maps on topographic or planimetric bases on quadrangle or irreg-
engineers. Included are reports on the results of resource studies and of
ular areas at various scales. Pre-1971 maps show bedrock geology in re-
topographic, hydrologic, and geologic investigations. They also include
lation to specific mining or mineral-deposit problems; post-1971 maps
collections of related papers addressing different aspects of a single sci-
are primarily black-and-white maps on various subjects such as environ-
entific topic.
mental studies or wilderness mineral investigations.
Bulletins contain significant data and interpretations that are of
Hydrologic Investigations Atlases are multicolored or black-and-
lasting scientific interest but are generally more limited in scope or geo-
white maps on topographic or planimetric bases presenting a wide range
graphic coverage than Professional Papers. They include the results of
of geohydrologic data of both regular and irregular areas; the principal
resource studies and of geologic and topographic investigations; as well
scale is 1:24,000, and regional studies are at 1:250,000 scale or smaller.
as collections of short papers related to a specific topic.
Water-Supply Papers are comprehensive reports that present sig- Catalogs
nificant interpretive results of hydrologic investigations of wide interest
Permanent catalogs, as well as some others, giving comprehensive
to professional geologists, hydrologists, and engineers. The series covers
investigations in all phases of hydrology, including hydrology, availabil- listings of U.S. Geological Survey publications are available under the
ity of water, quality of water, and use of water. conditions indicated below from USGS Map Distribution, Box 25286,
Building 810, Denver Federal Center, Denver, CO 80225. (See latest
Circulars present administrative information or important scien-
tific information of wide popular interest in a format designed for distri- Price and Availability List.)
bution at no cost to the public. Information is usually of short-term "Publications of the Geological Survey, 1879-1961" may be pur-
interest. chased by mail and over the counter in paperback book form and as a set
Water-Resources Investigations Reports are papers of an inter- microfiche.
pretive nature made available to the public outside the formal USGS "Publications of the Geological Survey, 1962-1970" may be pur-
publications series. Copies are reproduced on request unlike formal chased by mail and over the counter in paperback book form and as a set
USGS publications, and they are also available for public inspection at of microfiche.
depositories indicated in USGS catalogs. "Publications of the U.S. Geological Survey, 1971-1981" may be
purchased by mail and over the counter in paperback book form (two
Open-File Reports include unpublished manuscript reports, maps,
and other material that are made available for public consultation at de- volumes, publications listing and index) and as a set of microfiche.
Supplements for 1982, 1983, 1984, 1985, 1986, and for subse-
positories. They are a nonpermanent form of publication that may be cit-
quent years since the last permanent catalog may be purchased by mail
ed in other publications as sources of information.
and over the counter in paperback book form.
Maps State catalogs, "List of U.S. Geological Survey Geologic and Wa-
ter-Supply Reports and Maps For (State)," may be purchased by mail
Geologic Quadrangle Maps are multicolor geologic maps on to-
and over the counter in paperback booklet form only.
pographic bases in 7 112- or 15-minute quadrangle formats (scales main-
"Price and Availability List of U.S. Geological Survey Publica-
ly 1:24,000 or 1:62,500) showing bedrock, surficial, or engineering
tions," issued annually, is available free of charge in paperback booklet
geology. Maps generally include brief texts; some maps include structure
form only.
and columnar sections only.
Selected copies of a monthly catalog "New Publications of the
Geophysical Investigations Maps are on topographic or plani-
U.S. Geological Survey" is available free of charge by mail or may be
metric bases at various scales, they show results of surveys using geo-
obtained over the counter in paperback booklet form only. Those wish-
ph~sical techniques, such as gravity, magnetic, seismic, or radioactivity,
ing a free subscription to the monthly catalog "New Publications of the
~[Link] reflect subsurface structures that are of economic or geologic sig-
U.S. Geological Survey" should write to the U.S. Geological Survey,
mftcance. Many maps include correlations with the geology.
582 National Center, Reston, VA 22092.
Miscellaneous Investigations Series Maps are on planimetric or
topographic bases of regular and irregular areas at various scales; they Note.-Prices of Government publications listed in older catalogs,
present a wide variety of format and subject matter. The series also in- announcements, and publications may be incorrect. Therefore, the prices
cludes 7 1/2-minute quadrangle photogeologic maps on planimetric charged may differ from the prices in catalogs, announcements, and pub-
bases which show geology as interpreted from aerial photographs. The lications.
series also includes maps of Mars and the Moon.