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 released 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 

Books 

Professional Papers, Bulletins, Water-Supply Papers, Techniques of Water-Resources Investigations, Circulars, publications 

of general interest (such as leaflets, pamphlets, booklets), single 

copies of Earthquakes & Volcanoes, Preliminary Determination of 

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 

U.S. Geological Survey, Map Distribution 

Box 25286, MS 306, Federal Center 

Denver, CO 80225 

Subscriptions to periodicals (Earthquakes & Volcanoes and 

Preliminary Determination of Epicenters) can be obtained ONLY 

from the 

Superintendent of Documents 

Government Printing Office 

Washington, DC 20402 

(Check or money order must be payable to Superintendent of 

Documents.) 

Maps 

For maps, address mail orders to 

U.S. Geological Survey, Map Distribution 

Box 25286, Bldg. 810, Federal Center 

Denver, CO 80225 

Residents of Alaska may order maps from 

U.S. Geological Survey, Earth Science Information Center 

101 Twelfth Ave., Box 12 

Fairbanks, AK 99701 

OVER THE COUNTER 

Books and Maps 

Books and maps of the U.S. Geological Survey are available 

over the counter at the following U.S. Geological Survey offices, all 

of which are authorized agents of the Superintendent of Documents. 

• ANCHORAGE, Alaska-Rm. 101, 4230 University Dr. 

• LAKEWOOD, Colorado-Federal Center, Bldg. 810 

• MENLO PARK, California-Bldg. 3, Rm. 3128, 

345 Middlefield Rd. 

• RESTON, Virginia-USGS National Center, Rm. 1C402, 

12201 Sunrise Valley Dr. 

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

125 South State St. 

• SPOKANE, Washington-U.S. Post Office Bldg., Rm. 135, 

West 904 Riverside Ave. 

• WASIDNGTON, D.C.-Main Interior Bldg., Rm. 2650, 

18th and C Sts., NW. 

Maps Only 

Maps may be purchased over the counter at the following 

U.S. Geological Survey offices: 

• FAIRBANKS, Alaska-New Federal Bldg, 101 Twelfth Ave. 

• ROLLA, Missouri-1400 Independence Rd. 

• 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 velocity stack is originally defined as the summation 

of a seismic event along a hyperbolic moveout trajectory, but 

this definition is extended to include both linear and parabolic moveouts in this report. Thus, our velocity stack can be 

used for modeling events with a variety of moveouts. The 

ability to model via a velocity stack provides us a number of 

applications in the area of seismic data processing, such as 

multiple suppression, coherency filtering, trace interpolation, and velocity analysis. To take advantage of the computational efficiency in the frequency domain, only linear and 

parabolic moveouts are considered here. 

Velocity stacking can be accomplished through the use 

of either a transpose operator, such as in conventional velocity analysis, or by an inverse operator. These operators 

depend only on the source-receiver offsets of seismic traces 

and the range of velocities. The choice between the transpose 

and the inverse operator for computational purposes is immaterial because either operator need only be computed once. 

Some advantages in using the inverse operator over the transpose operator include preserving the input waveform (the 

transpose operator tends to bias the frequency content toward 

lower frequencies) and minimizing any mixing effects. These 

favorable properties of the inverse operator also allow for 

better coherency filtering and trace interpolation. 

Because velocity stacking works well at all sourcereceiver offsets, its performance in suppressing multiples is 

better than two-dimensional filtering in the frequency-wavenumber domain. When working with good-quality seismic 

data, using an inverse operator for velocity analysis has significant advantages over a transpose operator. However, as 

the quality of the data decreases with the addition of more 

noise, the use of the inverse operator is only marginally better. 

INTRODUCTION 

Normal moveout (NMO) correction and stacking are 

standard, industry, seismic-data-processing sequences. 

Combining these two processes is identical to stacking along 

a family of hyperbolic curves and is described as a velocity 

stack using a transpose operator (Thorson, 1984; Thorson 

and Claerbout, 1985). Conventional constant velocity analysis is one form of velocity stack using a transpose operator. 

It is a linear transformation and maps the data into velocity 

space. One problem in using a transpose operator is the 

smearing effect it has on the data. To better discriminate signal and noise in velocity space, Thorson ( 1984) proposed an 

inverse op~rator, which essentially deconvolves the spatial 

mixing effect of the transpose operator. 

Thorson ( 1984) described two approaches in estimating 

an inverse operator: the generalized inverse and stochastic 

inverse. The original formulation of Thorson ( 1984) requires 

inversion of a large matrix. Hampson ( 1986) showed that, 

for a typical 96-trace seismic shot of 3-s duration, the matrix 

order would be 150,000 by 150,000. Thorson and Claerbout 

(1985) suggested an iterative solution using a conjugate-gradient method, which is still time consuming but more efficient. In order to reduce computational time in estimating an 

inverse operator further, Hampson (1986) proposed a parabolic trajectory rather than a hyperbolic trajectory. The parabolic trajectory is not a physically observable traveltime 

curve in the offset domain, but it can be approximated by 

using NMO corrections with appropriate velocities. This 

approximation enables the use of the Fourier transform, and 

a computationally efficient method can be derived. 

Computing an inverse operator is more time consuming 

than computing a transpose operator. However, because both 

computations depend only on offset and number of velocities 

(or slownesses ), these operators need only be computed 

once, and, hence, the choice between the two methods based 

on computing time alone is immaterial. 

To a first-order approximation, seismic data in the offset domain can be considered as a linear superposition of 

many events whose traveltimes describe certain trajectories. The collection of events are defined as a model in the 

velocity domain. The inverse operator allows us to accurately estimate such a model from seismic data contaminated with random noise, and it provides us with a number 

of processing techniques, such as multiple removal, coherency filtering, trace interpolation, and velocity analysis. 
2 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

The performances of both the inverse and transpose operators for each application are compared, and the advantages 

of the inverse operator are discussed. 

Thorson (1984) described m!lltiple suppression by 

stacking along the primary velocity in velocity space. He utilized the property of better velocity discrimination of an 

inverse operator in his multiple suppression. Another 

approach to multiple suppression is to model the multiple 

energy by an inverse operator and subtract this multiple 

energy from the original data. This approach of multiple suppression by velocity stack is described well by Hampson 

(1986), and the actual performance for a field data set is documented by Geist and others (1990). This paper employs 

Hampson's approach and uses models to quantify the performance of multiple suppression-we emphasize the use of 

static corrections for better multiple elimination. 

The random noise in seismic data cannot be modeled as 

coherent events using a velocity stack. By transforming the 

data into model space and back to the offset domain, random 

noise will be reduced. This property can be used as a coherency-filtering or a signal-enhancement technique. Kong and 

others (1985) and Monk (1992) described nonlinear, coherency-filtering techniques using the invertible Radon transform or 't-p transform. Seismic events can also be predicted 

along various moveouts in model space, and this property 

can be used in trace interpolation. An inverse operator used 

for coherency filtering and trace interpolation produces data 

with reduced trace mixing (smearing) and preserves the 

high-frequency components. 

The compact response, or reduced side lobes, of an 

inverse operator in the slowness domain is ideally suited for 

high-resolution velocity analysis (Thorson, 1984 ). However, 

noise introduced into real data sets degrades the performance 

of an inverse operator. This limitation in the presence of 

noisy data will be discussed later. 

ACKNOWLEDGMENTS 

We thank John J. Miller and William J. Stephenson 

for constructive reviews. All data were processed using 

CogniSeis DISCO software with additional programs 

developed by U.S. Geological Survey personnel. 

THEORY 

Detailed theoretical development of inverse velocity 

stacking can be found in Thorson (1984) and Hampson 

(1986). To a first-order approximation, considering only 

arrival times, seismic data can be modeled as a superposition 

of laterally coherent events. For example, a plane-wave 

response of dipping layers consists of laterally coherent 

events whose stepout times are linear with distance. Also, a 

seismic event recorded on a horizontal-receiver spread from 

a shot can be modeled as an event whose traveltime is represented by a hyperbolic curve. The only difference between 

the above two examples is the functional relationship of 

arrival times with respect to the appropriate coordinate system. In the data space or offset space, we use offset, h, and 

two-way time, t, as variables, and in the model space or 

velocity space we use slowness, p, and zero-offset time, 't. 

Therefore, the seismic data d(h,t) can be represented by: 

d (h, t) = J J u (p, 1:) () (f(t, h, 1:, p)) dpdt + n (h, t) 

(1a) 

where 

d(h,t) is the measured seismogram at offset, h, and twoway time, t, 

u(p,'t) are seismic models with slowness, p, and zero-offset time, 1:, 

() is the Dirac delta function, 

flt,h,'t,p) is the traveltime relationship between data and 

model spaces, and 

n(h,t) is the measurement noise at offset, h, and two-way 

·time, t. 

Three functional relationships, flt,h,'t,p), in equation 

1a are considered by numerous authors. Stoffa and others 

(1981) usedflt,h,'t,p) = 1:-(t-ph) as their functional relationship (this is the well-known slant stack or 1:-p operation), 

whereas Thorson (1984) used a hyperbolic relationship, 

flt,h,'t,p) = 1:-J?- p

2

h

2 , for his inverse velocity stack. 

Hampson (1986) used a parabolic approximation to a 

hyperbolic traveltime curve (which isflt,h,'t,p) = 't-(t-sh2), 

with a new parameter, s) in his study of inverse velocity 

stack for multiple removal. The meaning of s will be discussed later. In this section, it is treated as a parameter 

similar to slowness, p. 

Thorson (1984) defined the integral relationship in 

equation 1 a 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: 

T u (p, 1:) = L d (h, t) = 

f f d (h, t) () (f(t, h, 't, p)) dhdt (lb) 

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 function of 1:. Equation 1 b indicates that performing normal 

moveout correction and stacking, or stacking along the 

moveout defined by flt,h,'t,p), are methods of transforming 

data into velocity space or model space. Under what conditions is this operation invertible? In other words, what operator will correctly invert the velocity space into the data 

space? As pointed out by Thorson (1984), the transpose operator defined in equation 1 b does not work well for a truncated 

data set, such as any real data set. The purpose of this section 

is to derive an efficient inverse operator that performs better 

than the transpose operator defined in equation 1 b. 
THEORY 3 

A computationally efficient inverse operator is difficult to derive if a hyperbolic moveout is assumed in equation lb. But, for a linear moveout or parabolic moveout, an 

efficient method can be derived in the frequency domain 

(Hampson, 1986). If we use a parabolic relationship for fin 

equation 1 a, then: 

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

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

Equation 2 indicates that the seismic data can be modeled 

by summing the model along a parabolic trajectory. Equation 2 also indicates that the data space can be modeled by 

a convolution of the model u(p,t) with the parabolic trajectory in the time domain. This convolutional property using 

the parabolic or linear traveltime curve in equation 1 is useful in the frequency domain because convolution in the 

time domain can be represented by multiplication in the 

frequency domain. 

D (hp ro) 

-iros 1h: -iros2h: 

e e 

D (h2, ro) -iros h~ -iros h~ e e 

D (hN, ro) -iros h~ -iros h~ e e 

·e 

·e 

Let us take a Fourier transform of equation 2, and let 

upper-case variables represent the Fourier-transformed 

variables. Then: 

I -irot d(h, t) e dt = 

D (h, ro) = I I u (s, t) e-irote -irosh

2

dsdt + N (h, ro) 

J -irosh2 

= U (s, ro) e ds + N (h, ro) (3) 

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(ht.ro),D(h2,ro), .... ,D(hN,ro)) 

with dimension N, 

U is a column vector of (U(st.ro),U(s2,ro), .... ,U(sM,ro)) 

with dimension M, 

N is a column vectorof(N(ht.ro),N(h2,ro), .... ,N(hN,ro)), and 

L is an NxM matrix with element of Lnm = e-irosmh

2

n. 

Explicitly, equation 4 can be written: 

-irosMh: U(spro) N (hp ro) 

-irosMh~ U (s2, ro) N(h2, ro) 

+ (5) 

-irosMh~ ·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, 

from the noisy discrete data, D, while minimizing the misfit 

function or noise, N. Thorson and Claerbout (1985) demonstrated that this problem can be solved as an overdetermined linear system. The standard least-squares solution to 

equation 4 is: 

(6) 

Equation 6 shows the desirable inverse operator, G = 

(LTL )-1 L~ in the velocity stack. 

The meaning of the s in the parabolic moveout is not as 

straightforward as the p in the linear or hyperbolic moveout 

because this is not a physically observable moveout. The linear relationship between the time variables ( 't and t) makes 

possible the frequency domain operation in equation 1. 

Under which circumstances can a hyperbolic moveout be 

approximated by a parabolic moveout that is linear in 't and 

t? The hyperbolic-moveout equation can be expanded using 

a Taylor series: 

J2 22 J 22 2 't + p h = 't 1 + (p h I 't ) 

If (p2h2!t2) << 1, dropping higher order terms can be justified. The only adjustable variable that will satisfy this condition is the slowness variable, p, and this adjustment can be 

accomplished by applying a normal moveout (NMO) correction with a slowness, Pc• as shown in Hampson (1986). 

Therefore the time (t') after the NMO correction with Pc can 

be written using equation 7: 

Therefore, if Pr is very small or p is close to Pc• a parabolic 

approximation can be justified. The parameter s defined in 

equation 2 is (p,2/2't), and, strictly speaking, it is not a slowness parameter. Because the coefficient h2 in the parabolic 

approximation is a function of zero-offset time, the shape of 

a parabola for a given s is a function of zero-offset time, different from the slant stack. Also, smearing can be expected 

as events deviate from the ideal parabolic form. 

INVERSE AND TRANSPOSE OPERA TORS 

Denoting the element of inverse operator G as G ij• 

equation 6 can be written as: 

U(sp ro) G11G12 · G1M D (hp ro) 

U (s2, ro) G2I G22 · G2M D (h2, ro) 

= (9) 

GNIGN2 

For a given frequency, each element of G is a complicated 

function of h and s, which is easily seen from equation 5, 

but does not involve actual data. In other words, for the 

same field geometry (same offset) and range of slownesses, 

the same inverse operator can be applied to the data set to 

compute U or D. Therefore, we can pre-compute an inverse 

operator or a transpose operator before applying it to the 

data. Because we only have to compute the operator once, 

the fact that computing an inverse operator is more expensive that computing the transpose operator is rather trivial. 

Figure 1A shows the elements of an inverse operator, 

Gu, convolved with a zero-phase band-pass wavelet for a 

parabolic moveout, and figure 1B shows the transpose operator. 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 

seen from figure 1, the elements of the transpose operator is 

essentially a band-pass wavelet located at a trajectory of the 

parabolic moveout, whereas those of the inverse operator are 

a very complicated function of offset and slowness along the 

parabolic moveout trajectory. Notice the waveform changes, 

including polarity reversal, for the inverse operator. The 

waveforms of the transpose operator are identical for all offsets 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 different from each other, so the application of an inverse operator is similar to a weighted stack along a parabolic moveout. 

If the numerical values of D in equation 9 are identical 

for all offsets, as for a horizontal event with the same wavelets, the model or velocity space can be derived by summing the rows of the inverse operator or summing along the 

offsets. Summing along the offsets is the same operation as 

summing each slowness panel in figure 1. The output of 

model space for the operators shown in figure 1 is shown in 

figure 2A for the inverse operator and is shown in figure 2C 

for the transpose operator. The output of the inverse operator is symmetrical in time, but this is not true for the transpose operator. The amplitude of the inverse operator decays 

quickly with increasing moveout time, whereas the amplitude of the transpose operator decays slowly with increasing moveout time. The increased lateral spread of the 

transpose operator in the slowness domain causes lateral 

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 ~ 

~( ( l( • ( 

u 

0 ~ 

• 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 parabolic 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 

the transpose operator can be easily examined by summing 

U over all slownesses. The result is shown in figure 2B for 

the inverse operator and in figure 2D for the transpose operator. Equation 3 indicates that integration with respect to 

slowness in the model space with h = 0 will result in the 

zero-offset data in the offset domain. Therefore, if accurate 

U is modeled from D (data in the offset domain), a summation of the model U along the slowness will reproduce the 

input data. The result of this operation for the inverse operator, shown in figure 2B, clearly shows the waveform as identical to the input, which is a zero-phase band-pass wavelet 

with 4/6-851100 Hz, but the result of the transpose operator 

shows a bias toward lower frequency. This is one reason 

why the conventional velocity-analysis method, summing 

along the assumed velocity trajectory, shows a lack of resolution, as mentioned in Thorson (1984). For horizontal 

events, this result can be represented by: 

M M N 

D (h = 0, ro) = L u (sm, ro) = L L Gnm (10) 

m=l m=ln=l 

This equation can be used to examine whether the computed 

inverse operator, G, is correct or not. 

An application of inverse and transpose operators in 

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 · 

-60 0 60 

i 

U) 

E 

0 

0 

1 

A 

I I 

MOVEOUT TIME. IN MllliSECONDS 

-60 0 60 

! 

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 

using the inverse operator; and figure 3C shows the result of 

applying the transpose operator. As indicated in figure 2, the 

result of the transpose operator fails to model the higher frequency part of the input data. Loss of the high-frequency 

component is evident at the far offsets. Results from using an 

inverse operator, however, show an almost perfect reproduction of the input data. 

APPLICATIONS 

The ability to extract model space from a noisy data set 

can be used for a variety of seismic-data-processing techniques, and a number of applications of inverse velocity 

stack, such as velocity analysis, is discussed in Thorson 

(1984). Hampson (1986) presented an application to remove 

multiples using a parabolic approximation. In this paper, 

applications of the inverse and transpose operators will be 

discussed in the area of multiple suppression, coherency filtering, velocity analysis, and spatial interpolation. 

MULTIPLE SUPPRESSION 

The methodology of multiple suppression using the 

velocity stack can be summarized in the following five steps: 

1. Apply NMO correction using stacking velocities of 

multiples. 

2. Model the multiple event using the inverse or transpose operator. 

3. Subtract the modeled multiple events in step 2 from 

the original data. 

4. Remove the NMO correction applied in step 1. 

5. Stack with the primary velocity. 

The performance of multiple suppression outlined above is 

analyzed using models and a real data set. 

MODEL DESCRIPTION 

In order to analyze the performance of the multiple suppression technique by a velocity stack, a 24-channel common midpoint (CMP) gather, with near offset of 400 m and 

far offset of 2, 700 m, was generated. Figure 4 shows the synthetic model filtered with a 4/6-42/54-Hz symmetrical 

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 perf ormance of the multiple suppression technique, and events 

below 2.5 s represent multiple energy. The average rootmean-square (RMS) amplitude of multiple energy is very 

close to the RMS amplitude of the reference trace. The RMS 

velocity of multiples is modeled by a linearly increasing 

function with depth: 1,400 m/s at 2,500 ms and 1,600 m/s at 

4,500 ms. Four primary events with identical amplitude and 

polarity are chosen and shown in figure 4B. The amplitude 

of the primaries in figure 4B is 20 percent of the reference 

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 

3,950 ms with 2,200 m/s. 

In our model study, figure 4A was mixed with figure 4B 

for a variety of scale factors. Three synthetic models were 

generated for the analysis. 

Modell: Multiple strength 14 dB higher than primaries. 

Model 2: Multiple strength almost equal to primaries. 

Model3: Multiple strength 20 dB higher than primaries. 

Figure 5A shows model 1, and figure 5B shows model 

2. Notice that, in modell, primary events cannot be observed 

by the eye. The amplitude relationship in figure 5A is similar 

to the field data example presented later. Figure 5C is model 

1 with random noise added and will be used later. 
APPLICATIONS 7 

400 1,600 2,600 400 

5.00 

A B 

OFFSET. IN METERS 

1,600 2.600 400 

c 

1,600 2,600 

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 

The results of straight stacking for noise-free models 

are shown in the "a" parts of figure 6, and results of multiple 

suppression by a velocity stack using the inverse operator 

(inverse velocity stack) are shown in the "c" parts of figure 

6. The results for the noise-free model (the "c" parts of fig. 

6) indicate that multiple suppression by the inverse velocity 

stack suppresses the multiple energy enough to recognize the 

primaries. Also, the result shown for model2, "a" part, indicates that straight stacking will be sufficient to suppress multiples when the multiple strength is almost the same as that 

of the primaries (model 2). However, it demonstrates that, if 

the multiple strength is more than 14 dB higher than the primaries, straight stacking is not enough to enhance the primaries (fig. 6, models 1 and 3, "a" parts), and some kind of prestack multiple-suppression technique, such as an inverse 

velocity stack, is required. 
8 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

(/) 

0 

:z 

0 

(..') 

. ~ 3.00 ~ 

u.i 

~ 

i= 

4.00 

OFFSET. IN METERS 

400 1,000 1,600 2,200 400 1,000 1,600 2,200 

5.00 ---t+t++t+tt+H-H-++t++t-+t+t---i-++t++++++-H++-++1-++-<i-++-1--l--

A B 

4 

Figure 4. Synthetic model filtered with a 4/6--42/54-Hz symmetrical wavelet. A, Multiple model except the event at 2 s, which 

serves as a reference amplitude; B, Primary model with four 

events. See text for additional explanation. 

In order to see the details of multiple suppression, an 

amplitude analysis of the result of velocity stack for model 2 

(fig. 6, model 2, "a" part) was performed, and the result is 

shown in figure 7 A. This figure indicates that the velocity 

stack suppressed the multiple energy on the order of 40 dB. 

Notice the five peaks near the -10-dB line in figure 7A, indicating almost perfect primary recovery. On the other hand, a 

straight stack only suppress the multiple energy on the order 

of 12 dB. Approximately a 28-dB enhancement is achieved 

by the velocity stack. 

The above analysis is for the ideal case in which we 

know exact multiple and primary velocities. For this ideal 

case, there is neither random nor other coherent noise, and 

each channel has an identical waveform. These ideal conditions are never met with real data. 

NOISE MODEL 

In order to evaluate the performance of the velocity 

stack in a "real" data situation, noise in the form of random 

static shifts, in the range of 4 ms, and random amplitude with 

Gaussian random noise were introduced into the CMP 

gather. The static shift was chosen for our noise model for 

the following reasons: 

1. The velocity stack or other multiple-suppression technique 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. 

4. Static shifts account for a variation of streamer depth 

during data acquisition . 

The performance of velocity stack with static shifts is 

also shown in figure 6. The "b" parts of figure 6 represent 

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 

stacking (compare the "a" parts with the "b" parts of fig. 6). 

However, the performance of the inverse velocity stack is 

highly dependent on the static shift. In the "d" parts of figure 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 7 B. 

The strength of primary events in figure 7A (noise free) and 

figure 7B (static noise) are similar, but the amount of multiple suppression under the static shifts is on the order of 20 

dB (it is about 40 dB in the noise-free case). In other words, 

if the multiple energy contaminated by static shifts is more 

than approximately 20 dB higher than that of primaries, it 

is difficult to enhance primaries over multiples using the 

inverse velocity stack. 

Notice that the strength of multiple energy is almost the 

same as that of the primaries in figure 6, model 3, "d" part, 

after multiple suppression. Recalling that the multiple 

strength before the multiple suppression for model 3 is 20 dB 

higher than the primary, the signal-to-noise ratio for model3 

can be considered as a threshold for recognizing the primaries. However, in order to differentiate primaries from multiples, the amplitudes of primary events after multiple 

suppression should be a couple of decibels higher than those 

of the multiples. Figure 6, modell, "d" part, whose multiple 

strength is about 14 dB higher than the primaries, clearly 

shows the primaries after multiple suppression even though 

there are 4-ms static shifts. Therefore, in processing real data 

sets, we could not enhance the primaries if the strength of 

multiples is about 17 dB greater than that of the primaries 

(see fig. 6, model1, "d" part and fig. 6, model3, "d"·part). 

Approximately a 20-dB loss of performance of the 

velocity stack for the noise model is due to small static shifts. 

In theory, static shifts can be computed and static corrections 

can be applied before the velocity stack. One of the results of 

applying static shifts computed by non-surface-consistent 
APPLICATIONS 9 

400 1.600 2.600 400 

5.00 

A B 

OFFSET. IN METERS 

1.frdo 2.600 400 

c 

1,600 2,600 

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 

figure 8 for model3. As expected, the result (fig. 8D) is comparable to the result for the noise-free model (fig. 8B). 

The effect of random-amplitude noise in multiple suppression is analyzed using model 1 by adding Gaussian random noise with maximum amplitudes to those of primaries 

(fig. 5C). The result of multiple suppression using an 

inverse operator is shown in figure 9. Figure 9A shows the 

result of straight stacking with random noise, and figure 9B 

shows the result of straight stacking with both random 

noise and static shifts in the range of 4 ms. There is not 

much difference between figures 9A and 9B. Both are too 

noisy to recognize primary events. The result of applying 

an inverse operator to the random noise section without 

static shifts is shown in figure 9C. All primary events are 

clearly observable. The effect of both random-amplitude 
10 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

1111111 IIIII 1111111 111111 1111111 IIIII 

MODE~ .DE. 

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

(( (( 

.. ;.. 

a b c d Ill 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 

of application of static correction is shown in figure 9 E. 

Comparing with figure 6 for model 1, it is obvious that 

noise in the form of 4-ms static shifts has a more pronounced effect than an equal amount of random noise. 

The result of multiple suppression using the transpose 

operator is shown in figure 10. In comparing figures 9C and 

9E with figures IOC and IOE, the advantage of the inverse 

operator over the transpose operator is evident. But when 

there is too much noise, the results shown in figure 9D and 

figure 1 OD imply that there is not much difference between 

the inverse and transpose operators. Hence, for high-quality 

input data, the inverse operator works better than the transpose operator. But, as the amount of noise (particularly in the 

form of static time shifts) increases, the advantage of the 

inverse operator over th~ transpose operator decreases. 

REAL DATA EXAMPLE 

Figure 11 shows an original stacked section generated 

using a standard processing sequence. This data set was 

acquired in Lake Baikal using a tuned air-gun array as a 

source. Unstacked data consists of 24-fold CMP gathers 

with a near offset distance of 400 m and a far offset distance of 2, 775 m. 

Our model study indicates that correcting for static 

shifts is a very important step in suppressing multiple 

energy. Thus, one additional step was added to the steps outlined at the outset of the multiple suppression section: static 

computation and application after step 1. The static computation was performed in the CMP domain. 

Figure 12 shows the result of multiple suppression. Figure 12A is the result of applying two-dimensional velocity 

filtering, or F-K filtering, a technique similar to the one 

described by Ryu (1982), and figure 12B shows the result of 

the application of the inverse velocity stack. For the inverse 

velocity stack, the multiples are modeled with seven slownesses, and the range of parabolic moveouts is -20 ms to 40 

ms at the far offset. Both techniques performed adequately 

for the multiple suppression below about 4,600 ms. Overall 

signal-to-noise ratio, however, is better in the section produced by the application of inverse velocity stack. The 

details of figure 12 near CMP 2,500 is shown in figure 13, 

and this figure clearly demonstrates that the inverse velocity 
APPLICATIONS 11 

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

-20 

-10 -

- T 

-10 f.- T 

~ L.Ucn 

C/) C/) 

-LU 

Ocr:: 

Zo... -20 f.-

0... 

C/) 

.....J 

LU 

CD 

u 

LU 

0 

~ -30 uJ 

0 

~ 

1-

:::::i 

0... 

::2 

<( 

-40 

~ z 

.~ 0 

Ci5 

C/) 

LU 

0:: 

- a: 

~ 

C/) 

LU 

C/) 

.] 

(5 

z 

-

~ ~· ~: ~ u ... · ~-. f.-

~ .. v ~· 

I J I • I ·. A 

-50 

-60 

.. ~ 

C/) 

... .. 

~·· ·.·~ -~~~ ~ 

·.~ .. ~ ·. 

~ ~ 

I I I ~ 

-30 -

-40-

-50 -

-60--------~ .. L-----L--------L-----L~ 2 3 4 5 2 3 4 5 

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. 

Notice much stronger and clearer dipping events are imaged 

in figure 13B. 

We tested many seismic profiles using the inverse 

velocity stack versus the F-K domain approach. In almost 

every case, the inverse velocity stack worked better than 

the F-K domain approach, particularly for data sets with 

high noise and low frequency. A similar conclusion was 

reached by Geist and others (1990) in their study of longperiod (> 1 s), first-order, multiple removal over a continental slope. Better performance of the inverse velocity 

stack for the lower frequency range has an important implication for deep crustal studies. As mentioned in Hampson 

(1986), the better performance of the inverse velocity stack 

is due to the fact that inverse operators work well for all 

offsets, whereas the F-K approach leaves significant multiple energy on the near-offset traces. 

COHERENCY FILTERING 

It is well known that stacking suppresses random noise 

proportional to the root of the number of traces in the stack. 

If all reflections are horizontal, this is a good method for 

enhancing signal-to-noise ratio (this is the fundamental concept applied in CMP stacking). In the stacked section, this is 

not an appropriate method because there are numerous nonhorizontal events. In order to enhance dipping events, stacking should be done along different moveouts. Stacking along 

the linear moveouts can be effectively implemented using a 

slant stack in the 't-p domain. This kind of coherency filtering is necessary for data with a low signal-to-noise ratio. 

In this section we will discuss coherency filtering in the 

context of model extraction from noisy data. As indicated in 

equation 6, the model space can be estimated by applying the 

inverse operator to the offset-domain data set. This estimated 
12 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

C/) 

Cl 

:z 

0 

(..) 

1.00 

2.00 

~ 3.00 ~ 

u.i 

~ 

i= 

4.00 

5.00 

J 

\' 

{ 

) 

t 

fr 

( ~ 

,, 

I 

.~ 

I 

!'I 

•)•) 

H 

'i 

) ) ) ) 

A B c D 

Figure 8. Multiple suppression for model 3 using an inverse operator in the velocity stack. A, Straight stack without static shift; 

B, Multiple suppression without static shift; C, Multiple suppression 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. 

en 

0 

z 

0 

u 

1.00 

2.00 

~ 3.00 ~ 

w 

~ 

~ 

4.00 

5.00 

~ 

~ 

~ ~ 

{ 

~\· ·r; 

"' 

~ 

'1 

A 

~ ~ 

"i'((lft'ft ~ 

'( it ( 

< 

( 

;. 

. . P' 

~~~ ~ 

! 

B c D E 

Figure 9. Multiple suppression for model 1 with randomamplitude 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 random noise only; B, Straight stack with random noise and static 

shift; C, Multiple suppression with random noise only; D, Multiple suppression with random noise and static shift; E, Multiple 

suppression with random noise and static shift after static correction with "trim static" (see text) in the CMP domain. 
APPLICATIONS 13 

en 

a 

z 

0 

(._) 

1.00 

2.00 

~ 3.00 ~ 

u.i 

~ 

i= 

4.00 

5.00 

~ 

~ 

f ~ 

{ · 

('\W~i· 

~ 

tUi 

A 

~ ----- ~ 

~~~ ~~» 

~ t < ~ ~ 

<'(, :t ( . 

•}' 

"' 

~~~1~} 

B c D E 

Figure 10. Multiple suppression for model 1 with randomamplitude 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 random noise only; B, Straight stack with random noise and static 

shift; C, Multiple suppression with random noise only; D, Multiple suppression with random noise and static shift; E, Multiple 

suppression with random noise and static shift after static correction with "trim static" (see text) in the CMP domain. 

en 

0 

z 

0 

(._) 

UJ 

en / 

u.i 

~ 

i= 

2.00 

4.00 

5.00 

6.00 

1,375 

COMMON MIDPOINT 

1,975 2,575 

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 acquisition parameters. 
14 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

1,375 1,975 2,575 

COMMON MIDPOINT 

1,375 1,975 2,575 

2.00 

4.00 

Cl) 

c 

:z 

0 

u 

U.J 5.00 ·.·· Cl) 

~ 

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 twodimensional 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 

domain data set using equation 4 without the noise term. The 

presence of the noise in equation 4 indicates that the transform generated by equation 6 is not a complete space. Some 

of the noise, n(h,t), that has not been included in U is lost 

when transforming back to the offset domain. This operation 

can be written as: 

(11) 

where 

fi is the coherency-filtered output. 

Figure 14A shows the result of applying an inverse 

operator as a coherency filter to the section shown figure 

12B, and figure 14B shows the result of applying a transpose 

operator. The inverse or transpose operator was designed 

using linear moveout with 11 slowness steps and five offsets. 

The range of linear moveout at the far offset is ±24 ms (the 

offset of the center trace is 0 m). The output is estimated at 

the center trace. Both figures 14A and 14B show a reduction 

of background random noise or the enhancement of local 

coherency estimated within five traces. Details of coherency 

filtering in the time domain are shown in figure 15, and figure 16 shows the frequency-wavenumber domain representation of the data shown in figure 15. The low-frequency 

enhancement and mixing effects produced by the transpose 

operator are apparent in figures 15 and 16. 

TRACE INTERPOLATION 

The concept of coherency filtering by inverse velocity 

stack can be applied to trace interpolation. The only difference is the application of the L operator in equation 11. In 

equation 11, the element of L applying to model space U to 

transform back to the offset domain is the same element used 
APPLICATIONS 15 

en 

Cl 

z 

0 

(..) 

5.00 

2,275 

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

~ 

uJ 

~ 

i= 

7.00 

A 

2.475 2,675 

COMMON MIDPOINT 

2.275 

B 

2.475 2,675 

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 

be used for the output. In other words, a different offset 

value, not included in the input data, can be used for the 

interpolation. Assume that a trace is interpolated using two 

input traces, and the interpolation is done at the middle of 

two traces using linear moveout. Then the L used for computing G can be written as: 

(12) 

The element L* for the interpolation can be written as: 

l -iws 1h1* -iws2h 1*) L* = e e (13) 

where 

h* =(hi+ h2)/2 

An example of trace interpolation using an inverse 

operator is given in figure 17. The trace interval of figure 17 

is 25 m, and five traces were used for the interpolation. An 

inverse operator was designed to pass seismic events whose 

linear moveouts are within ±4 ms per trace. Highly dipping 

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 

the interpolated section. 

VELOCITY ANALYSIS 

Moveouts in the CMP domain depend on the interval 

velocities of the layered Earth, and the model space is estimated in the slowness domain. Thus, the inverse velocity 

stack is a logical choice for velocity analysis. When using 

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 correction with velocity, Vc (or slowness, Pc-see eq. 8). 

Parameter "sh2" in parabolic moveouts is presented as 

the moveout time at the far-offset distance; thus, a relationship is required to derive an actual velocity from the velocity 

stack using parabolic moveouts. The relationship between 

apparent and true velocity can be derived from equation 8: 

Using the parameters: 

2 

v 

1 1 

-+-

2 2 2 

v vr vc 

(14) 

Because the parameter is given by the moveout time at the 

far offset, equation 14 can be written as: 
16 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

1,375 1,975 2,575 

COMMON MIDPOINT 

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. 

v = 

where 

Mt is moveout time at the far offset, and 

hmax is the far-offset distance. 

(15) 

As indicated in equation 15, the true velocity depends on the 

zero-offset time as well as the correction velocity for NMO 

and moveout time. 

A velocity analysis for the noise-free model shown in 

figure 5B is given in figure 18 with linearly increasing velocity ofvc = 1,500 m/s at 2,500 ms and Vc = 1,700 m/s at 4,500 

ms. The range of moveout times at the far offset is -240 ms 

to 120 ms, with a 12-ms interval. The output of the application of the inverse operator, shown in figure 18A, clearly 

shows linearly increasing velocity for the multiples and four 

primary events. The moveout time of P4 is -168 ms at twoway 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. However, the output from the transpose operator or conventional 

velocity analysis, shown in figure 18B, show a smearing 

effect of the adjacent moveouts. The above analysis was 

applied to the noisy model, random noise (primary signal-torandom-noise ratio of about 2), and static shifts (in the range 

of 4 ms), and the result is shown in figure 19. Although the 

result from the inverse operator is better than the result from 

the transpose operator, the difference is not as dramatic as 

that shown in figure 18. 

An application to a real data set is shown in figure 20. 

The result using the application of an inverse operator (fig. 

20A) produces better resolution but is noisier than the result 

produced by the transpose operator (fig. 20B) 
2.275 2.475 2.675 

5.00 

6.00 en 

0 

z 

0 

u 

UJ 

en 

z 

-

L..iJ 

::E 

i= 

7.00 

A 

2.275 

B 

COMMON MIDPOINT 

2.475 2,675 2.275 2.475 2.675 

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> 

::g 

t""' 

....... 

n 

;J> 

>-3 

....... 

0 

z 

en 

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

20.00 

15.00 

q 

~· N 

li: 

~~I U.J 

:I: 

z 

;:. 10.00 

I ~~ . ~ ~ (_) 

z I U.J 

~~ :::::> ~ 

~· 0 

U.J 

a: 

u.. 

5.00 

0.00 

-0.5 0.0 0.5 -0.5 

A B 

I 

,, II II, 

Itt, 

0.0 

WAVENUMBER 

0.5 -0.5 0.0 0.5 

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 

Many conventional seismic processing methods are 

based on the transpose operator or conjugate operator 

(Claerbout, 1992). Although the transpose operator applied 

to perfect data generally gives an imperfect result, sometimes the transpose operator performs better than the 

inverse operator because the conjugate operator tolerates 

imperfections of input data owing to the limited aperture, 

truncation, nonstationary noise, missing data, etc. On the 

other hand, an inverse operator gives a perfect result for a 

perfect data set but has the possibility of noise amplification when using actual data and requires more information 

to do a better job. As can be seen from equation 6, the 

transpose operator is the first step of inversion, and, if 

(LTL)-1 is close to an identity matrix, the inverse operator 

can be approximated by a transpose operator. 

To reduce noise amplification of the inverse operator 

applied to an imperfect data set or to stabilize the inverse 

operator, a random noise factor, A, such as the one used for 

spiking deconvolution, can be applied. When implementing 

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

the computation of the inverse operator. All examples in this 

paper are generated by adding 1 percent of each diagonal element to the diagonal terms. The amount of noise can be 

adjusted to tune the performance of the inverse operator. 

Thorson and Claerbout (1985) defined this form of inverse 

as "stochastic inverse" (Aki and Richards, 1980) and suggest 

using the ratio of noise variance to variance of the points in 

the model space for A. The method proposed by Thorson 

and Claerbout (1985) may provide an optimum inverse operator for each data set, but the operator cannot be pre-computed because A cannot be determined without data. 

Two basic assumptions in inverse velocity stacking are 

that, along various moveouts, events are aligned and have 

uniform amplitudes. Multiple suppression using real data 

sets shows improved results after static time corrections are 

applied to improve alignment. However, the improvements 

were marginal in most cases. One example is shown in figure 

21. This suggests that simple static corrections alone will not 

suffice in certain cases. This is because our waveform 

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, preprocessing, such as amplitude balancing and waveform 

deconvolution, will be beneficial to the performance of the 

inverse operator. 

When inverse velocity stacking is applied for irregular 

offset gathers similar to land data or in a marine data set, 

whose CMP fold is less than the number of channels, it is 

impractical to compute different inverse operators for each 

offset gather. For example, 24-fold marine data acquired 

with 96 channels have 4 sets of offsets in the CMP domain, 

and 4 sets of inverse operators are required. To handle irregular geometries, one effective method is to apply differential 

NMO correction to a reference offset gather. 

The estimation of model space from the noisy input 

data by means of the velocity stack is defined as a process of 

coherency filtering. Coherency filtering improves the signalto-noise ratio, as shown in figure 15, and helps subsequent 

processing steps, such as migration. Most signal-enhancement techniques use weighted velocity stacks. Kong and 

others (1985) proposed a nonlinear signal estimator by multiplying semblance (Taner and Koehler, 1969) to the slantstacked data and inverse-transforming the weighted slantstack data. Monk (1992) presented a diversity slant stack, 

which is a weighted slant stack based on the power in the 

slant-stack domain. Both Kong and others (1985) and Monk 

(1992) applied a transpose operator in equation 11. In order 

to preserve the original frequency content, which is altered 

by applying a transpose operator, Monk (1992) and Kong 

and others (1985) applied an appropriate Rho filter (Claerbout, 1984) to the output. The main purpose of weighting, 

proposed by Monk (1992), is to reduce the smearing effect 

of the transpose operator and to modify the transpose operator, more or less, to look similar to an inverse operator. A 

weighted slant stack using an inverse operator provides an 

alternative to the conventional approach to signal enhancement. The price we have to pay for a weighted slant stack for 

coherency filtering in this formulation is that U should be 

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'ir rn•-'E-f'r'!Jlnr~r'rlr ,~ 

~cP> )-, ~f 

en 

Cl 

z 

0 

u 

UJ 

en 

~ 

ui 

~ 

i= 

A 

MULTIPLE 

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 

the velocity stack implicitly uses a local coherency and dip 

estimation (or slowness in this paper) and is subject to the 

aliasing problem when interpolation distance is too large. 

The formulation of interpolation is cast as a prediction in 

the offset domain in this paper, but the essential operation 

is similar to interpolating the amplitude along the direction 

of lateral coherence of each event (Barden, 1987; Lamer 

and others, 1981). This approach differs from Spitz (1991), 

who derived a method based on the fact that linear events 

may be interpolated exactly, regardless of spatial interval, 

without any information regarding their true dip. He interpolated traces using forward-backward, one-step prediction 

filters in the frequency domain. However, the proposed 

method has the ability to keep high-frequency components 

(as opposed to other interpolation methods based on the 

transpose operator). 

One problem for the velocity analysis of real data using 

an inverse operator is low signal-to-noise ratio compared to 

the output by a transpose operator. One way to overcome this 

difficulty is to adjust the random noise component when 

deriving the inverse operator. More random noise in the 

design of the inverse operator increases the signal-to-noise 

ratio but also increases lateral smearing. An optimum tradeoff between resolution and smearing can be found by trial 

and error. Figure 22 shows a velocity analysis with an 

inverse operator using 2 percent noise for A, instead of 1 percent used for figure 20A with static corrections. The output 
(I) 

c 

z 

8 

LU 

Cf) 

~ 

u.i 

~ 

i= 

-240 

A 

DISCUSSION 

MOVEOUT TIME. IN MILLISECONDS 

120 -240 

B 

21 

120 

< i 

I ) } ~ ~ 

< < L< ~ ~ r 1 J 1 

''~( I •• rr'''''''f ((~~~~({ f 

: 

~ 

• 

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 

en 

c 

z 

0 

(_) 

2.00 

2.50 

~ 3.00 ~ 

w 

~ 

i= 

en 

c 

z 

0 

(_) 

3.50 

4.00 

2.00 

2.50 

~ 3.00 ~ 

w 

~ 

i= 

3.50 

4.00 

-100 

A 

B 

APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

-50 

MOVEOUT TIME. IN MILLISECONDS 

0 50 100 

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. 
1-

z 

0 

0... 

0 

g 

M 

~g zM 

0 

~ 

~ 

0 

(..) 

0 

0 

M 

0 

0 

N 

DISCUSSION 

SONOJ3S Nl '3V\IIl 

0 

0 

Lri 

23 
24 APPLICATIONS OF VELOCITY -STACK METHODS TO SEISMIC DATA PROCESSING 

-100 -50 

MOVEOUT TIME. IN MILLISECONDS 

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 

different from the result shown in figure 20B. However, the 

result shown in figure 22 demonstrates that the improved 

performance of the inverse operator can be obtained with 

static corrections and by tuning the inverse operator with 

varied random noise. 

CONCLUSIONS 

A number of applications of velocity stack, such as 

multiple suppression, coherency filtering, trace interpolation, and velocity analysis, are described in this paper. The 

advantage of using an inverse operator over a transpose 

operator is in its ability to accurately estimate the model 

space from noisy seismic data. A transpose operator has a 

bias toward lower frequencies and more spread in the slowness domain. These properties result in trace mixing or 

smearing in the velocity stack. However, the advantage of 

using an inverse operator over the transpose operator 

decreases as noise increases. Based on this study, the following conclusions can be derived. 

1. Because of the effectiveness of velocity stack at all 

offsets using an inverse operator, multiple suppression by velocity stack is better than the frequencywavenumber domain approach. 

2. Because velocity stack relies on the alignment of seismic events along the moveouts, static correction of 

the CMP gathers is important in multiple suppression. 

Model studies indicate that about 4-ms static shifts in 

a gather degrades multiple suppression about 20 dB. 

3. Coherency filtering or trace interpolation by an 

inverse operator can be considered as a high-resolution method and minimizes trace mixing inherent to 

output by the transpose operator. 

4. Velocity analysis by an inverse operator shows high 

discrimination ability of close velocities, but the 

amount of noise, such as static shifts, controls the performance rather dramatically. To obtain high-resolution velocity analysis for field data, pre-processing, 

such as static correction and amplitude adjustment, 

should be applied. 

REFERENCES CITED 

Aki, K., and Richards, P.G., 1980, Quantitative Seismology, Theory 

and Methods: San Francisco, W.H. Freeman Company, 911 p. 

Barden, V., 1987, Trace interpolation in seismic data processing: 

Geophysical Prospecting, v. 35, p. 343-358. 

Claerbout, J.F., 1984, Imaging the Earth's Interior: Blackwell Scientific Publications, 398 p. 
REFERENCES CITED 25 

---1992, Earth Sounding Analysis, Processing Versus Inversion: Blackwell Scientific Publications, 304 p. 

Geist, E.L., Hart, P.E., and Fisher, M.A., 1990, A comparison of 

multiple suppression methods across a continental slope: Eos, 

Transactions, American Geophysical Union, v. 71, p. 1589. 

Hampson, Dan, 1986, Inverse velocity stacking for multiple elimination: Journal of the Canadian Society of Exploration Geophysicists, v. 22, p. 44-55. 

Kong, S.M., Phinney, R.A., and Roy-Chowdhury K., 1985, A nonlinear signal detector for enhancement of noisy seismic record 

sections: Geophysics, v. 50, p, 539-550. 

Lamer, K., Gibson, B., and Rothman, D., 1981, Trace interpolation and the design of seismic surveys [abs.]: Geophysics, v. 

46,p.407. ' 

Monk, D.J., 1992, Approach to optimum slant stack and its application as a seismic noise attenuator [expanded abs.]: Society of 

Exploration Geophysicists, 62nd Annual International Meeting 

and Exposition, p. 1161-1164. 

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 

Ryu, J.V., 1982, Decomposition (DECOM) approach applied to 

wave field analysis with seismic reflection records: Geophysics, v. 47, p. 869-883. 

Spitz, S., 1991, Seismic trace interpolation in the F-X domain: Geophysics, v. 56, p. 785-794. 

Stoffa, P.L., Diebold, B.P., and Wenzel, F., 1981, Direct mapping 

of seismic data to the domain of intercept time and ray parameter-A plane wave decomposition: Geophysics, 46, p. 

255-267. 

Taner, M.T., and Koehler, F., 1969, Velocity spectra-Digital computer derivation and applications of velocity function: Geophysics, v. 34, 859-881. 

Thorson, J.R., 1984, Velocity stack and slant stack inversion methods: Palo Alto, Calif., Stanford University, Ph.D. thesis. 

Thorson, J.R., and Claerbout, J.F., 1985, Velocity stack and slant 

stack stochastic inversion: Geophysics, v. 50. p. 2727-2741. 



SELECTED SERIES OF U.S. GEOLOGICAL SURVEY PUBLICATIONS 

Periodicals 

Earthquakes & Volcanoes (issued bimonthly). 

Preliminary Determination of Epicenters (issued monthly). 

Technical Books and Reports 

Professional Papers are mainly comprehensive scientific reports 

of wide and lasting interest and importance to professional scientists and 

engineers. Included are reports on the results of resource studies and of 

topographic, hydrologic, and geologic investigations. They also include 

collections of related papers addressing different aspects of a single scientific topic. 

Bulletins contain significant data and interpretations that are of 

lasting scientific interest but are generally more limited in scope or geographic coverage than Professional Papers. They include the results of 

resource studies and of geologic and topographic investigations; as well 

as collections of short papers related to a specific topic. 

Water-Supply Papers are comprehensive reports that present significant interpretive results of hydrologic investigations of wide interest 

to professional geologists, hydrologists, and engineers. The series covers 

investigations in all phases of hydrology, including hydrology, availability of water, quality of water, and use of water. 

Circulars present administrative information or important scientific information of wide popular interest in a format designed for distribution at no cost to the public. Information is usually of short-term 

interest. 

Water-Resources Investigations Reports are papers of an interpretive nature made available to the public outside the formal USGS 

publications series. Copies are reproduced on request unlike formal 

USGS publications, and they are also available for public inspection at 

depositories indicated in USGS catalogs. 

Open-File Reports include unpublished manuscript reports, maps, 

and other material that are made available for public consultation at depositories. They are a nonpermanent form of publication that may be cited in other publications as sources of information. 

Maps 

Geologic Quadrangle Maps are multicolor geologic maps on topographic bases in 7 112- or 15-minute quadrangle formats (scales mainly 1:24,000 or 1 :62,500) showing bedrock, surficial, or engineering 

geology. Maps generally include brief texts; some maps include structure 

and columnar sections only. 

Geophysical Investigations Maps are on topographic or planimetric bases at various scales, they show results of surveys using geoph~sical techniques, such as gravity, magnetic, seismic, or radioactivity, 

~h.tch reflect subsurface structures that are of economic or geologic sigmftcance. Many maps include correlations with the geology. 

Miscellaneous Investigations Series Maps are on planimetric or 

topographic bases of regular and irregular areas at various scales; they 

present a wide variety of format and subject matter. The series also includes 7 1/2-minute quadrangle photogeologic maps on planimetric 

bases which show geology as interpreted from aerial photographs. The 

series also includes maps of Mars and the Moon. 

Coal Investigations Maps are geologic maps on topographic or 

planimetric bases at various scales showing bedrock or surficial geology, 

stratigraphy, and structural relations in certain coal-resource areas. 

Oil and Gas Investigations Charts show stratigraphic information for certain oil and gas fields and other areas having petroleum potential. 

Miscellaneous Field Studies Maps are multicolor or black-andwhite maps on topographic or planimetric bases on quadrangle or irregular areas at various scales. Pre-1971 maps show bedrock geology in relation to specific mining or mineral-deposit problems; post-1971 maps 

are primarily black-and-white maps on various subjects such as environmental studies or wilderness mineral investigations. 

Hydrologic Investigations Atlases are multicolored or black-andwhite maps on topographic or planimetric bases presenting a wide range 

of geohydrologic data of both regular and irregular areas; the principal 

scale is 1:24,000, and regional studies are at 1:250,000 scale or smaller. 

Catalogs 

Permanent catalogs, as well as some others, giving comprehensive 

listings of U.S. Geological Survey publications are available under the 

conditions indicated below from USGS Map Distribution, Box 25286, 

Building 810, Denver Federal Center, Denver, CO 80225. (See latest 

Price and Availability List.) 

"Publications of the Geological Survey, 1879-1961" may be purchased by mail and over the counter in paperback book form and as a set 

microfiche. 

"Publications of the Geological Survey, 1962-1970" may be purchased by mail and over the counter in paperback book form and as a set 

of microfiche. 

"Publications of the U.S. Geological Survey, 1971-1981" may be 

purchased by mail and over the counter in paperback book form (two 

volumes, publications listing and index) and as a set of microfiche. 

Supplements for 1982, 1983, 1984, 1985, 1986, and for subsequent years since the last permanent catalog may be purchased by mail 

and over the counter in paperback book form. 

State catalogs, "List of U.S. Geological Survey Geologic and Water-Supply Reports and Maps For (State)," may be purchased by mail 

and over the counter in paperback booklet form only. 

"Price and Availability List of U.S. Geological Survey Publications," issued annually, is available free of charge in paperback booklet 

form only. 

Selected copies of a monthly catalog "New Publications of the 

U.S. Geological Survey" is available free of charge by mail or may be 

obtained over the counter in paperback booklet form only. Those wishing a free subscription to the monthly catalog "New Publications of the 

U.S. Geological Survey" should write to the U.S. Geological Survey, 

582 National Center, Reston, VA 22092. 

Note.-Prices of Government publications listed in older catalogs, 

announcements, and publications may be incorrect. Therefore, the prices 

charged may differ from the prices in catalogs, announcements, and publications.