Linearized Inversion of 

Reflection Traveltimes 

U.S. GEOLOGICAL SURVEY BULLETIN 1965 
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 are listed in various U.S. Geological Survey catalogs (see back inside cover) but not listed in the most recent annual "Price and 

Availability List" are no longer available. 

Prices of reports released to the open files are given in the listing "U.S. Geological Survey Open-File Reports," updated monthly, which is for sale in microfiche from the U.S. Geological Survey, Books and Open-File Reports Section, Federal Center, Box 

25425, Denver, CO 80225. 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 given 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 the foregoing series that have gone 

out of print at the Superintendent of Documents, are obtainable by mail 

from 

U.S. Geological Survey, Books and Open-File Reports 

Federal Center, Box 25425 

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, D.C. 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 

Federal Center, Box 25286 

Denver, CO 80225 

Residents of Alaska may order maps from 

Alaska Distribution Section, U.S. Geological Survey, 

New Federal Building - Box 12 

101 Twelfth Ave., Fairbanks, AK 99701 

OVER THE COUNTER 

Books 

Books of the U.S. Geological Survey are available over the 

counter at the following Geological Survey Public Inquiries Offices, all 

of which are authorized agents of the Superintendent of Documents: 

• WASHINGTON, D.C.--Main Interior Bldg., 2600 corridor, 

18th and C Sts., NW. 

• DENVER, Colorado--Federal Bldg., Rm. 169, 1961 Stout St. 

• LOS ANGELES, California--Federal Bldg., Rm. 7638, 300 N. 

Los Angeles St. 

• MENLO PARK, California--Bldg. 3 (Stop 533), Rm. 3128, 

345 Middlefield Rd. 

• RESTON, Virglnia--503 National Center, Rm. 1C402, 12201 

Sunrise Valley Dr. 

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

South State St. 

• SAN FRANCISCO, California--Customhouse, Rm. 504, 555 

Battery St. 

• SPOKANE, Washington--U.S. Courthouse, Rm. 678, West 

920 Riverside Ave .. 

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

• ANCHORAGE, Alaska--Federal Bldg, Rm. E-146, 701 C St. 

Maps 

Maps may be purchased over the counter at the U.S. Geological Survey offices where books are sold (all addresses in above list) and 

at the following Geological Survey offices: 

• ROLLA, Mlssouri--1400 Independence Rd. 

• DENVER, Colorado--Map Distribution, Bldg. 810, Federal 

Center 

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

Reflection Traveltimes 

By MYUNG W. LEE and 

WARREN F. AGENA 

U.S. GEOLOGICAL SURVEY BULLETIN 1965 
U.S. DEPARTMENT OF THE INTERIOR 

MANUEL LUJAN, JR., Secretary 

U.S. GEOLOGICAL SURVEY 

Dallas L. Peck, Director 

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. 

UNITED STATES GOVERNMENT PRINTING OFFICE: 1991 

For sale by the 

Books and Open-File Reports Section 

U.S. Geological Survey 

Federal Center 

Box 25425 

Denver, CO 80225 

Library of Congress Cataloging-in-Publication Data 

lee, Myung W. 

Linearized inversion of reflection traveltimes I by Myung W. Lee and Warren 

F. Agena. 

p. em. - (U.S. Geological Survey bulletin; 1965) 

Includes bibliographical references. 

1. Seismi.c traveltime inversion. 2. Seismic reflection 

method. 3. Seismology-Mathematical models. I. Agena, Warren 

F. II. Title. Ill. Series. 

QE75.B9 no. 1965 

[QE539 .2.S43] 

557.3 s--dc20 91-7790 

[551.2'2] CIP 
CONTENTS 

Abstract 1 

Introduction 1 

Acknowledgments 2 

Theory 2 

Computational aspect 4 

Comparison of methods 5 

Accuracy of inversion 6 

Field configuration 9 

CMP domain inversion 11 

Mode conversion 13 

Real data examples 15 

Discussion 19 

Conclusions 19 

Appendix A 20 

Resolution matrix 20 

References cited 21 

FIGURES 

1. Schematic ray-path diagram for a reflection from the Nth interface 3 

2. Ray tracings for a split-spread shot gather for the shot domain inversion 

test 7 

3. Graph showing the root-mean-square (RMS) error of inversion 9 

4. Graph showing smallest singular values for a one-layer model 9 

5-7. Ray tracings for: 

5. Shot gathers with geophone group interval of 0.2 km with 20 

channels 11 

6. Common midpoint (CMP) gathers for CMP domain inversion 13 

7. Shot domain inversion including mode conversion at the reflecting 

boundary 15 

8. Graphs showing measured and computed arrival times for marine multichannel 

data 17 

9. Migrated stacked section of sign-bit data 18 

10. Graph of measured and computed arrival times for shot domain inversion 20 

TABLES 

1. Layer parameters for model 1 and initial estimate for traveltime inversion 6 

2. Inversion results using split-spread shooting geometry 8 

3. Comparison of inversion results using different random noise content 10 

4. Comparison of inversion results between split-spread and off-end shooting 

geometries 12 

5. Inversion results and Dix result for common midpoint shooting geometry and 

variable amounts of noise contamination 14 

6. Inversion results including mode conversions for split-spread shooting 

geometry 16 

7. Comparison of inversion results using real marine seismic data 17 

Contents Ill 

Linearized Inversion of Reflection Traveltimes 

By Myung W. Lee and Warren F. Agena 

Abstract 

A computationally accurate and efficient traveltime 

inversion method for a plane-layered earth model can be 

formulated using a generalized linear-inversion technique. 

Sui~able parameterization enables us to estimate layer 

vanables such as velocity, dip, and depth simultaneously for 

all layers by iteratively solving the reduced-size matrix. 

Although computationally similar to the layer-stripping 

method, our method offers advantages in terms of convergence 

and consistency of the traveltime inversion. 

This method was applied both in the shot and common 

midpoint (CMP) domains, and results indicate that the shot 

domain inversion is more sensitive to measurement noise 

compared to the CMP domain inversion. We used a 

regularization method to overcome the adverse effects of 

measurement noise, but the method often resulted in 

inconsistent parameter estimation. This suggests that care must 

be taken when implementing regularization. 

The formulas presented are also applicable to the traveltime inversion of mode-converted waves when longitudinalor shear-wave velocity is known. 

INTRODUCTION 

An accurate estimation of interval velocities and 

reflector geometry enables us to perform detailed geological 

analysis of seismic data. Interval velocities are extremely 

important in constraining lithologic interpretations of deep 

crustal seismic data. Conventional velocity estimation in 

reflection seismology is performed in the common midpoint 

(CMP) domain under the assumption of hyperbolic moveout 

and negligible dip angle. This method works well when 

offset distance is small compared to target depth and when 

reflections have moderate dips. In order to overcome these 

limitations, a linearized inversion method in both the CMP 

domain and shot domain is presented in this paper. 

Traveltime inversion is one of many non-linear geophysical inversion problems and has been applied for a long 

time to observed earthquake data in order to investigate the 

Manuscript approved for publication November 30, 1990. 

internal structure of the earth (Aki and Richards, 1980). 

Generalized linear-inversion theory has also been used to 

solve the non-linear traveltime inversion problem (Backus 

and Gilbert, 1970; Aki and Lee, 1976). Aki and Richards 

(1980), Tarantola (1987), and Lines (1988) have written 

papers that contain excellent discussions of generalized 

inversion theory. Lines (1988) also compiled significant 

recent papers related to geophysical inverse problems. 

Traveltime inversion can be formulated either in a 

distance-time (x-t) domain or in an intercept time-horizontal 

slowness ('t-p) domain. Diebold and Stoffa (1981) 

developed a 't-p traveltime inversion method valid for both 

reflections and refractions and Diebold and others (1988) 

applied this method to expanding-spread profile (ESP) data 

with offsets of 6.5 to 14 km. Diebold and others (1988) 

analyzed wide angle reflection and refraction data and 

computed interval velocities for intervals as deep as 30 km 

in the area of the Baltimore Canyon trough on the east coast 

of the United States. 

Liu and others (1986) applied the x-t domain 

inversion method to ESP data acquired in the Sevier desert 

of west-central Utah and estimated interval velocities in 

intervals as deep as 30 km. 

Inversion of reflection times in three dimensions 

using surface seismic data was discussed by Gjfi'Systdal and 

Ursin (1981). They estimated interval velocities and 

reflection interfaces using a layer-stripping method. Bishop 

and others ( 1985) applied a tomographic inversion 

technique in order to estimate interval velocity and reflector 

depth in a laterally varying medium using surface seismic 

profiles. A tomographic inversion technique to obtain twoand three-dimensional velocity structures of the subsurface 

using well logs, vertical seismic profiles (VSP's), and 

surface measurements was investigated by Chiu and Stewart 

( 1987). They also extended the inversion method to 

transversely isotropic media. 

In this paper, a traveltime inversion method in the x-t 

domain is formulated by suitable parameterization and by 

explicit expression of the Jacobian matrix. The formula was 

developed assuming iso-velocity layers separated by plane 

boundaries within the earth's crust. This paper posits that 

the general inverse problem of simultaneously estimating 

Introduction 
layer parameters is equivalent to iteratively computing the 

reduced-size matrix. This method is computationally similar 

to the layer-stripping method. In the case of a model with N 

layers containing M observations for each layer, 

conventional simultaneous parameter estimation involves 

the inversion of a 3N x NM matrix, whereas our method 

involves the inversion of a 3 x M matrix N times. By 

introducing explicit derivatives and suitable 

parameterization, the inversion method becomes computationally efficient and accurate. 

The method is suitable to the traveltime inversion of 

surface shot or CMP gathers and can easily accommodate 

the inversion of mode-converted waves. We tested the 

method on real marine and land multichannel seismic data 

to determine its feasibility and limitations. Traveltime 

inversion in terms of the domain of application, field 

configuration, resolution, and regularization are also 

discussed. 

ACKNOWLEDGMENTS 

We thank John J. Miller and William A. Schneider, 

Jr. for their constructive reviews. 

THEORY 

Many researchers have extensively investigated the 

application of linear inverse theory to traveltime inversion 

and other geophysical problems. Theoretical development 

for geophysical inversion has been intensive (Aki and 

Richards, 1980; Tarantola, 1987; Lines and Treitel, 1984). 

The details of traveltime inversion using the generalized 

inverse are given by Lee (1990). Thus, only brief discussion 

of theory pertinent to the understanding of this paper is 

given. 

The solution of nonlinear inverse problems generally 

begins with a Taylor series expansion of the computed 

model response for a given parameter distribution. 

Parameters are updated iteratively assuming convergence to 

a global minimum. Thus, traveltime can be approximated by 

the following series expansion that retains only the first 

order terms: 

where 

tN is the traveltime from the Nth reflector with 

dimension M, 

2 Linearized Inversion of Reflection Traveltimes 

P n+l is the (n+ 1)th estimation of parameters with 

dimension L, and 

5P1 is the parameter update given by: 

The least squares solution of 'OP1 in equation 1 for a 

given observation oj (i=1, 2, .... ' M) is attained by 

minimizing the mean-squared error (E): 

M 2 

E = L [ oi-t j (P n + 1

) ] 

i=1 

The formal solution of equation 2, as given by the least 

squares method, can be written as: 

where 

(3) 

o-l is the generalized inverse of theM XL Jacobian 

matrix, 

0 is the M-dimensional column vector of the 

observations, 

tn is the M-dimensional column vector of the 

modeled traveltime at the nth iteration, and 

P is the £-dimensional column vector of layer 

parameters. 

The conventional matrix inverse ( o-1

) exists only if G is 

square (M = L) and nonsingular. However, M >> L for most 

inverse problems, so the matrix inverse presented in this 

paper is identified with the general inverse case (see, e.g., 

Aki and Richards, 1980, chapter 12). As shown in equations 

2 and 3, the essence of traveltime inversion is the evaluation 

of the Jacobian matrix and its generalized inverse. Using the 

singular value decomposition method for the generalized 

inverse, the solution of equation 3 can be written as: 

(4) 
where 

V is an L x L matrix of orthonormal eigenvectors, 

U is an M x M matrix of orthonormal eigenvectors, and 

A is a diagonal matrix of L singular values. 

the accuracy and speed of the computation for the layeredearth model, an explicit expression of the Jacobian matrix is 

derived here. 

It is apparent from equation 3 that accurate estimation 

of the Jacobian matrix is important in order to solve 

nonlinear inversion problems. Generally, a finite difference 

approximation can be used for the evaluation of the Jacobian matrix (Lines and others, 1984). In order to increase 

Figure 1 shows a simple ray-path diagram for the 

derivation of the Jacobian matrix including mode 

conversion. We define the coordinate of the intersection of 

the downgoing ray with ith interface as (Xi, Zi) and for the 

up going ray as (J(, ~). The traveltime from the source to the 

kth receiver via the reflection at the Nth interface (T k) can 

SOURCE RECEIVER 

----------~~--------------------------------------~----------~x 

\ 

\ 

\ 

\ 

\ I 

\ I 

1st layer 

\ I --~---~------------/--------- (i-1 )th layer 

(Xi, Zi) 

\ 

\ 

\ 

\ 

\ 

\ 

\ I 

\I 

z 

I 

I 

I 

I 

ith layer 

Nth layer 

Figure 1. Schematic ray-path diagram for a reflection from the.-v~ interface. (X;, Z;) represents 

the intersection of a downgoing wave with the layer interface; (X& Z;) represents the intersection 

for upgoing waves. 5; is the slowness ofthe downgoing wave and~; is the slowness of the upgoing 

wave. 

Theory 3 
be written by: 

N N 

T~= L [siQ~,;+siQ~,i] =L s;Q~.;+U (5) ~I ~1 

where 

s; is the slowness of the downgoing wave in ith 

medium, 

s; is the slowness of the upgoing wave in the ith 

medium, and 

tl)(Nk . = XNk . -XNk ·-1 ,l ,l 'l 

t;Z,Nk . = zNk . - zNk ·-I ,l ,l ,l 

The tilde indicates a quantity corresponding to the upgoing 

ray path and is denoted as "U'' from here on for notational 

simplicity. As indicated in equation 5, if there is no mode 

conversion at the reflecting boundary, s; = s;. 

As shown in equation 5, Q ~. i represents the length 

of ray segment between (i-1 )th and ith interface (or medium 

i) for the ray reflected from the Nth interface and received 

at the kth geophone. The dipping layer is defined using 

slope a; and intercept h; by: 

Thus, Z ~. i can be written as: 

zNk ·=a-XNk . +h· ,l l ,l l (6) 

Using the above notation, the required derivatives shown in 

equation 2 may be written as follows: 

(7) 

4 Linearized Inversion of Reflection Traveltimes 

Q k -1 A '7 k X k u 

Sj+ I ( N,j+ 1) l...l.L N,j+ 1 N,j+ 1 + 

In all cases in equation 7, j :::;; N. A similar derivation 

of the derivatives of traveltime with respect to layer 

parameters without mode conversion is given by Lee 

(1990). 

COMPUTATIONAL ASPECT 

If we assume N reflections measured with M geophones and no mode conversion at the reflecting boundary, 

then equation 3 can be written in the following matrix form: 

where 

(8) 

OP; is the ith layer parameter update defined as a 

column vector of (Os;. Oa;, Oh;), and 

Om; is a column vector of order M consisting of the 

difference between observed and computed arrival 

times from the ith layer. 

In equation 8, the matrix elements G; j are themselves 

matrices defined as: 

G;j = 

As indicated above, the matrix elements of G; j consist 

of the derivatives of the traveltime reflected from the jth 

interface with respect to the ith model parameters. By 

defining the parameter vector as shown in equation 8, the 

generalized inverse of equation 8 can be computed easily 

because matrix G is in a lower triangular form. The iterative 

solution of the general inverse is given as follows: 
(9) 

8PN = (GNN) -I ( 8mN- ~ GiN 8P) <=1 

Thus, instead of solving the 3N x MN matrix shown 

in equation 8, this arrangement of parameters leads to the 

inverse of a 3 x M matrix N times. This reduced matrix 

dimension reduces computation time and increases the 

numerical accuracy of the matrix inversion. We call this 

method the iterative whole-layer inversion (IWI) method. 

The layer -stripping method used by Lines and others ( 1984) 

is very similar to this approach from a computational point 

of view. The main advantages of the IWI method are: 

1. Computationally, it is equivalent to the layerstripping method. The only computational overhead 

is the additional matrix multiplication, shown in 

equation 9. 

2. Unlike the layer-stripping method, the errors of 

the upper layers do not propagate to the lower layers. 

3. In the case of high noise contamination with a 

poor initial estimate, the IWI method converges to the 

expected solution better than the layer-stripping 

method. This will be shown later in the model study. 

4. Because it involves a reduced size matrix, 

numerical accuracy and computational efficiency can 

be achieved. 

Some disadvantages of this approach compared to the 

layer-stripping technique are the slow convergence for a 

poorly defined initial estimate and more cross talk between 

parameters in an unconstrained problem. 

In practice, we do not use equation 9 directly because 

the iteration may not converge due to the non-linearity of 

the inverse problem. To circumvent this problem, we use a 

formula from Aki and Richards (1980): 

where 

0<~~1. 

8P= ~ G-1 8m (10) 

Numerical simulation indicates that ~ is a very 

important parameter when the initial estimate is poor. A 

suitable choice for ~ improves the result of inversion and 

will be discussed later. 

COMPARISON OF METHODS 

In the previous section, we briefly discussed the 

advantage of the IWI method over the layer-stripping 

method. In this section, the performance of the IWI method 

will be examined through numerical simulations (models). 

Table 1 shows the layer parameters used for model 1 

and the parameters for the initial estimate. The ray-path 

diagram for a shot gather with 0.2 km group intervals and a 

split-spread geometry is shown in figure 2A. The ray path 

for the initial estimate is shown in figure 2B. As can be seen 

from figure 2, the initial estimate is poorly defined for this 

model. The initial root-mean-square (RMS) error between 

the arrival times of the model and the initial estimate was 

1.757 s. Generalized inversion using the IWI method 

reduced the RMS error to the expected numerical accuracy 

of 3x1o-5 s in seven iterations. The layer-stripping method 

also produced the same result. 

In the case of low noise contamination, where traveltime data were contaminated with Gaussian random noise 

having a standard deviation of less than 8 ms (crt< 8 ms), 

both methods converged to similar numerical values. In the 

case of high noise contamination (0t = 16 ms), inversion 

results from both methods were quite different for certain 

layers even though both methods converged to the expected 

RMS error. Inversion results with crt = 16 ms are shown in 

table 2. Notice the differences in parameters for layer 4. The 

layer-stripping method provided a much lower interval 

velocity of 2.42 km/s compared with the true velocity of 4.8 

km/s and calculated a thickness of 30 m instead of the true 

value of 100m. The inversion result from the IWI method 

produced an interval velocity of 4.86 km/s and a layer 

thickness of 1 00 m. Except for this discrepancy, inversion 

results from both methods are quite similar. 

With slightly different initial conditions, the layerstripping method converged to different numerical values in 

many cases, but the IWI method converged to consistent 

numerical values in most of the cases we tried. For example, 

with ~ = 0.5 in equation 10 (instead of ~ = 0.42, which was 

used for the result shown in table 2), the inversion results by 

the layer-stripping method were very similar to those 

determined by the IWI method shown in table 2. 

In order to examine the effects of ~ on convergence to 

a solution and to test the consistency of the solution, the IWI 

method was applied to the model shown in figure 2A. The 

RMS error, with respect to iteration number and~. for 0t = 

16 ms is shown in figure 3. The computer program 

initializes ~ to 1.0, and, in succeeding iterations, if the RMS 

error is larger than the error of the previous iteration, ~ is 

adjusted in order to prevent divergence. Figure 3 shows an 

increase in RMS error from iteration 1 to 2. Thus, ~ was 

adjusted to the values shown in figure 3 to allow 

convergence. 

Figure 3 indicates that if ~ is too large or too small (~ 

= 0.2 or ~ = 1.0), the IWI method does not converge to the 

Comparison of Methods 5 
Table 1. Layer parameters for model 1 and initial estimate for traveltime inversion 

Velocity 

Layer (kilometers 12er second} Dip Depth 

number P-wave S-wave (degrees) (kilometers) 

Modell parameters 

1 1.5 1.0 7.0 0.7 

2 2.5 1.25 -5.0 1.3 

3 3.8 1.9 0.0 2.0 

4 4.8 2.5 0.0 2.1 

5 3.8 2.0 -10.0 3.0 

6 4.2 1.9 15.0 5.0 

7 3.5 1.6 0.0 7.0 

Initial estimate for P-P refl~tion inversion 

1 1.5 

2 2.0 

3 2.5 

4 3.0 

5 3.5 

6 4.0 

7 4.0 

solution. Also, the converged solution indicates consistent 

estimates of layer parameters. This figure also shows that 

the rate of convergence depends on the value of fl For 

example, the IWI method converged in five iterations for ~ 

= 0.6, but in seven iterations for ~ = 0.8. There seems to be 

an optimum value of~ for a given problem, but currently we 

do not know how to choose the optimum value at each 

iteration. 

The consistent numerical results achieved by the IWI 

method indicate a distinct advantage over results obtained 

by the layer-stripping method. This example demonstrates 

the flexibility of the IWI method in choosing the ~ 

parameter. 

ACCURACY OF INVERSION 

The accuracy of traveltime inversion depends 

primarily on the error in observation in relation to the 

shooting geometry. Numerical simulations using the model 

parameters shown in table 1 indicate that, for the noise-free 

case, the inversion converged to the true layer parameters 

even though the length of surface spread (L) was much 

smaller than the target depth (h). For example, using the 

initial estimates shown in table 1, the traveltime data with a 

0.05 km group interval and 20 channels converged to the 

true solution (this corresponds to L/h = 0.15 for the deepest 

target). However, with a small error in the observation (O't = 

2 ms ), the inversion result showed a large deviation from 

the exact layer parameters even though it converged to the 

6 Linearized Inversion of Reflection Traveltimes 

0.0 0.3 

0.0 2.0 

0.0 3.8 

0.0 5.0 

0.0 6.0 

0.0 6.5 

0.0 7.0 

expected accuracy (RMS error= 1.8 ms). The velocity, dip, 

and depth of layer 7 were calculated to be 3.18 km/s, 

-22.71°, and 10.22 km when crt= 2 ms. 

The quantitative aspect of inversion error due to 

errors in observation can be analyzed using equation 3. Let 

us assume that there is no measurement error and the 

inversion converged to the solution in the nth iteration. 

Then, (OP)n = (G-1

)n (Om)n = 0, since Om (the difference 

between observed and computed arrival times) is zero. But, 

if there is uncorrelated Gaussian random measurement 

error, the variance of the model parameters can be written 

as: 

where 

Var (OP) = Var (a-t om) 

Var represents the variance, and 

crt represents the standard deviation of the measurement error. 

The variance, from equation (4) and Jackson (1972), 

can be written as: 

where 

Var (OP) = L vt (~t) 2 (11) 

j 

A.j are the singular values of the Jacobian Matrix G, 

and 

vij are the orthonormal eigenvectors. 
;> 

(') 

(') 

5; 

~ '< 

a 

:; 

~ 

ii1 

cs· 

::s 

...... 

(J) 

a: 

w 

fw 

~ 

0 

OFFSET, IN KILOMETERS 

-3 -2 -1 0 2 3 -3 -2 -1 0 2 3 

0 1t i i l I· 1i_ i I I Llllliiiilk J J 1 IJ 2

1 

J -'i i -i i 1 tn l I I 0 

@ 

2 

=---- I \mmlllllmm I 

2 

~ 

g 1~ U\t' 1'tlAim .. ··xy',1 ' ( 1',fW,,!W'\,,I,\'i\'.'\jl I i r

1

'/'l r l \\\\\~\\\\\~\\\IV1tMJ#IIIIiiiiiii\\WifiiiVKIJIUII/Ii/ll//ll \---4 z 

- ~-~"""= ~z I 1 y i' I , ~x ' ' ' '}' {. \i'\li' '\1 \I i , I .. I ! I I 1 ~- ~· u/vl't. , v

1

1, ,, ·,fU, \\I\' 11''J\\I. i!: @ a... 

w 

0 

I 

.\\\ .. ~~~~rmititt'' \l',~·,.··r~·,,· .~ 1

1 I i r I \ j' iii!/ I , ,t ' 1

·

1 \! , ! i 1 \ l I .\:.r 1

111! :1 lh I II ,'::,:,k._,l I! 

§ 

I'---" '1 '1,l/~Url!llillli 1il'/ if J \~~~~~~~/ I 

6-~ \ ii\AI',!'ri!l!lli i,l::,,r, '1i,, r 11 1 6 

@ 

@ 

A B 8_J ________________________________________________ _L_ _J---------------------------------------------------L-8 

Figure 2. Ray tracings for a split-spread shot gather for the shot domain inversion test. The geophone group interval is 0.2 km and near-offset distance 

is 0.1 km. P-wave velocities for each layer (in km/s), from table 1, are circled. A, Model whose parameters are shown in table 1; 8, Initial estimate whose 

parameters are shown in table 1 . 
Table 2. Inversion results using split-spread shooting geometry 

[Shooting geometry shown in figure 2. Random noise content ( cr1

) equals 16 milliseconds; group 

interval equals 200 meters; near-offset distance equals 152 meters; 20 channels] 

Velocity 

Layer (kilometers per second) Dip Depth 

number P-wave (degrees) (kilometers) 

Layer-stripping method 

1 1.53 7.46 0.72 

2 2.52 -4.79 1.31 

3 3.96 -0.12 2.05 

4 2.42 0.58 2.08 

5 3.70 -10.35 2.97 

6 4.50 17.52 5.14 

7 6.90 -14.00 9.19 

Iterative whole-layer oonversion method (IWI) 

1 1.53 

2 2.52 

3 3.96 

4 4.86 

5 3.46 

6 4.46 

7 6.93 

Detailed derivation and discussions on measurement error 

may be found in van der Sluis and van der Vorst ( 1987). 

Equation 11 indicates that the error of the layer 

parameters depends on the singular values of the Jacobian 

matrix: as Aj becomes smaller, the effect on the accuracy 

becomes more pronounced. 

In order to get some insight into inversion accuracy in 

relation to the singular values, a simple one-layered model 

with velocity (v) = 2 km/s, dip (JZj) = 20°, and target depth 

(h) = 2 km was used. The smallest singular value (As) for 

this model with various L/h ratios is shown in figure 4. 

When L!h is equal to 1, As:::: 0.04. When L/h is less than 1, 

As decreases very rapidly with respect to decreasing L/h. 

When L/h is greater than 1, the variation of As with respect 

to the increasing L/h is flattened out. Thus, the inversion 

result of the previous example with 0.05 km group intervals 

could have a large uncertainty for the deepest layer where 

L/h:::: 0.15. . 

The smallest singular values (As) for the single-layer 

model considered previously ( v = 1 km/s) are indicated in 

figure 4 as squares at L/h= 1 and L/h = 2; As values for v = 

4 km/s are indicated as crosses. The general trend for v = 1 

km/s or v = 4 km/s is very similar to that of v = 2 km/s. The 

smallest singular value for the higher interval velocity is 

less than that for the lower interval velocity. This figure 

implies that the inversion results for smaller interval 

velocities have less variance than those for large interval 

velocities with the same measurement error. The analysis of 

8 Linearized Inversion of Reflection Traveltimes 

7.46 0.72 

-4.79 1.31 

-0.12 2.04 

1.44 2.14 

-8.71 2.95 

18.85 5.13 

-13.01 9.13 

singular values for different models leads us to conclude 

that the key factor controlling the accuracy of inversion is 

the arrival time difference between near- and far-offset 

distance. 

Table 3 shows a series of inversion results for model 

1 (table 1, fig. 1) with respect to variable amounts of noise 

contamination. When crt < 8 ms, errors of inversion for the · 

upper five layers are less than about 5 percent. The L/h ratio 

for the fifth layer is about 1.3, so L/h values for the upper 

five layers are greater than 1. Thus, the variances of layer 

parameters for the upper five layers are small and the 

inversion results are reliable. The results shown in table 3 

demonstrate the increasing error with respect to increasing 

measurement error and increasing depth. The largest error 

occurs for layer 7. In this case, even for a small error of crt 

= 4 ms, the velocity error is about 20 percent. 

In any inversion problem, very small singular values 

have a dramatic effect on the solution. In order to 

circumvent this problem, a regularization is sometimes 

employed in the inversion process. In the singular value 

decomposition method, regularization can be implemented 

by dropping the solutions corresponding to small singular 

values (van der Sluis and van der Vorst, 1987), but this 

introduces another kind of error called regularization error. 

As mentioned by van der Sluis and van der Vorst (1987), 

regularization reduces the variance of the solution but gives 
(/) 0.1 

0 

z 

0 

u 

w 

(/) 

~ 

a;' 

0 

a; 

a; 

w 

(/) 

~ 

a; 0.01 

2 4 6 8 10 12 

ITERATION NUMBER 

Figure 3. Graph showing the root-mean-square (RMS) error of 

inversion result when at = 16 ms versus iteration number and 

various damping factors, B. 

a biased solution. The regularization is closely related to the 

resolution of the solution and will be discussed later. 

Using the initial estimates from table 1 and the 

geometry shown in figure 2B, with at = 8 ms, inversion 

results for layer 7 produced parameter estimates of v = 3.49 

km/s, j2j = 1.81°, and h = 7.04 km when we dropped the 

solutions corresponding to A,s less than 0.03. Compare this 

result to the estimation without regularization shown in 

table 3. The inversion result for layer 7 with regularization 

is very close to the true solution, but the inversion result 

incorporating regularization is highly dependent on the 

initial estimate: it provides a non-unique solution. When we 

use h = 8 km as an initial estimate instead of the 7 km shown 

in table 1, the inversion results with the same regularization 

are v = 5.11 km/s, j2j = 5.68°, and h = 7.98 km. The standard 

deviation for layer 7 without regularization is the column 

vector of (31 O't, 30at, 220at), and with regularization it is the 

column vector of (0.08at, 1.11 at, 0.140'1

). Hence, there is a 

significant reduction in the variance of the solution with 

regularization. As indicated by this example, the problem of 

regularization is the inconsistency of the solution. Thus, 

unless there is a good estimate of the initial parameters, it is 

difficult to apply regularization to the inversion problem. 

All of the inversion results for the models shown in the 

tables in this paper were estimated without regularization. 

0.1 

0 

X 

:< 0.01 

0.001 

0.00010 

0 

X 

EXPlANATION 

• 2 kilometers per second 

X 4 kilometers per second 

0 1 kilometer per second 

2 3 

L/h 

Figure 4. Graph showing smallest singular values of the 

Jacobian matrix (A) in the singular value decomposition with 

respect to the ratio of spread length to target depth (ljh) for a 

one-layer model with the dip angle of 20°. Continuous line 

connecting dots is for velocity (v) = 2 km/s, squares for v = 1 

km/s, and crosses for v = 4 km/s. 

FIELD CONFIGURATION 

In marine data acquisition, an off-end shooting 

geometry is always employed, but, in land data acquisition, 

either off-end or split-spread geometry may be used. In 

order to increase the accuracy of the interval velocity 

estimation using the stacking velocity, off-end geometry is 

sometimes used in a land survey, because, with the same 

group interval and number of channels, longer offsets can be 

obtained for velocity analysis. 

In order to examine the effect of the offset distance on 

the traveltime inversion in the shot domain, the geometry 

shown in figure 5 is used. Figure 5 shows the ray paths for 

two off-end geometries and one split-spread geometry. The 

layer parameters for this model are identical to those shown 

in table 1 except for the first layer. Parameters for the first 

Field Configuration 9 
Table 3. Comparison of inversion results using different random noise content 

[Shooting geometry shown in figure 2. Random noise content, crt; group interval equals 200 meters; 

near-offset distance equals 152 meters; 20 channels] 

Velocity 

Layer (kilometers per second) Dip Depth 

number P-wave (degrees) (kilometers) 

Random noise rontent (ut) = 4 milliS«onds 

1 1.51 7.11 0.71 

2 2.51 -4.95 1.30 

3 3.84 -0.03 2.01 

4 4.87 0.36 2.11 

5 3.71 -9.69 2.99 

6 4.30 16.16 5.05 

7 4.23 -2.21 7.46 

Random noise tontent (ut) = 8 milliS«onds 

1 1.52 7.23 0.71 

2 2.51 -4.90 1.30 

3 3.88 -0.06 2.02 

4 4.91 -0.73 2.12 

5 3.64 -9.36 2.98 

6 4.34 16.97 5.07 

7 5.04 -5.31 7.94 

Random noise rontent (ut) = 12 millisetonds 

1 1.52 7.35 0.72 

2 2.52 -4.85 1.31 

3 3.91 -0.09 2.03 

4 4.91 1.10 2.13 

5 3.55 -9.03 2.96 

6 4.40 17.89 5.10 

7 5.92 -8.76 8.49 

Random noise tontent (ut) = 16 millisetonds 

1 1.53 7.46 0.72 

2 2.52 -4.79 1.31 

3 3.96 -0.12 2.04 

4 4.86 1.44 2.14 

5 3.46 -8.71 2.95 

6 4.46 18.85 5.13 

7 6.93 -13.01 9.13 

Random noise rontent (ut) = 20 milliS«<nds 

1 1.54 7.58 0.73 

2 2.53 -4.73 1.31 

3 4.00 -0.15 2.05 

4 4.80 1.74 2.15 

5 3.38 -8.42 2.93 

6 4.51 19.73 5.16 

7 8.12 -18.70 9.97 

1 0 Linearized Inversion of Reflection Traveltimes 
OFFSET, IN KILOMETERS 

-3 -1 0 2 3 -3 -2 -1 0 2 -2 -1 0 2 3 

(f) 

a: 

UJ 

f--

UJ 

::2: 

0 

...J 

2 

~ 4 

~ 

Ifa.. 

UJ 

0 

6 

A 8~------------------------_L 

0 

2 

4 

6 

B c __j__ ________________________ __j_B 

Figure 5. Ray tracings for shot gathers with geophone group interval of 0.2 km with 20 channels. Model parameters are identical 

to those shown in table 1 except for the first layer; parameters for the first layer are: velocity= 2 km/s, dip= 2°, and depth = 0.4 

km. A, Off-end geometry with a source at -2 km; 8, Split-spread geometry with a source at 0 km; C, Off-end geometry with a source 

at 2 km. 

layer are: v = 2.0 km/s, S(j = 2°, and h = 0.4 km. As can be 

seen in figure 5, identical group interval spacings and 

number of channels are employed. The far-offset distance 

for the off-end geometry is twice that of the split-spread 

configuration. The inversion results using the initial 

estimate shown in table 1 with crt = 8 ms are shown in table 

4. The results shown in table 4 demonstrate that there is no 

apparent advantage for employing an off-end geometry in 

order to increase the offset distance. These results also 

indicate that split-spread geometry is better than off-end 

geometry for estimating dip angles. 

In the case of horizontal layers, inversion results 

indicate that longer offset data are better for estimating layer 

parameters. Based on this model study, we conclude that the 

best way to increase the accuracy of inversion is to use the 

longest possible group intervals with a split-spread 

geometry for a given number of channels. 

CMP DOMAIN INVERSION 

The most common velocity analysis method in 

reflection seismology is performed in the CMP domain. 

Interval velocity can be computed by Dix's formula (Dix, 

1955) if it is assumed that the stacking velocity is 

approximately equal to the RMS velocity. When dip angles 

are low and the offset distance is small compared to target 

depth, this method provides a reliable interval velocity 

estimate. The interval velocity may also be estimated in the 

CMP domain by the linear-inversion method. 

In order to compare the shot domain inversion with 

the CMP domain inversion, the ray tracing for the CMP 

geometry based on model 1 was computed and the result 

shown in figure 6A. Figure 6B shows the CMP ray tracing 

result for model 1 with a low-velocity thin layer (v = 2.8 

km/s) instead of a high-velocity layer (v = 4.8 km/s) for 

layer 4. In principle, we could include dip estimation in the 

CMP domain inversion process. However, our model study 

indicates that the inclusion of dip estimation hampers the 

convergence of the solution and produces a large variance 

for the estimated dip angle. The subsurface ray coverage for 

CMP geometry is very small (fig. 6) compared to the 

subsurface ray coverage for shot domain geometry (fig. 5). 

Thus, the estimation of dip angle using CMP geometry is 

not practical. Also, in the CMP domain, numerical 

simulation indicates that the smallest singular values for the 

inverse of the Jacobian matrix for model 1 is very small (on 

the order of 1 0-6). Therefore, the best practical method for 

inverting CMP domain data requires estimating only 

interval velocities and depths. We can incorporate dip 

information from other sources, if available, into the 

inversion process. 

The inversion results for the CMP domain data with 

the same geometry as shown in figure 6 are shown in table 

5. When crt= 0, layer parameters were estimated within 4 

percent. This error is introduced because of the horizontal 

layer assumption of the CMP domain inversion. In this 

example, where the maximum dip is on the order of 15°, the 

CMP domain inversion provided reliable layer parameters. 

CMP Domain Inversion 11 
Table 4. Comparison of inversion results between split-spread and off-end shooting 

geometries 

[Random noise content (0'1) equals 8 milliseconds; group interval equals 200 meters; 20 channels] 

Velocity 

Layer (kilometers per second) Dip Depth 

number P-wave (degrees) (kilometers) 

Oft'-end shooting source = -2.0 kilometers 

1 2.02 5.24 0.50 

2 2.57 -4.76 1.33 

3 3.97 0.41 2.07 

4 4.89 0.65 2.17 

5 3.47 -10.40 2.94 

6 4.34 16.08 5.01 

7 5.39 -3.8 8.19 

Split-spread shooting source = 0 kilometers 

1 2.02 2.40 0.41 

2 2.51 -4.91 1.30 

3 3.89 -0.04 2.02 

4 4.91 -0.72 2.12 

5 3.60 -9.28 2.97 

6 4.34 17.08 5.06 

7 5.09 -5.42 7.96 

Oft'-end shooting source= 2.0 kilometers 

1 2.02 

2 2.58 

3 3.99 

4 4.88 

5 3.74 

6 4.25 

7 4.99 

Table 5 also contains the inversion result when 

measurement errors are introduced. Comparing these results 

with the shot domain inversion, it is obvious that the CMP 

domain inversion is not quite as sensitive to measurement 

noise. Particularly remarkable in this example is that the 

depth error of layer 7, even when crt = 32 ms, is less than a 

few percent. A comparison is also made in table 5 for the 

interval velocity estimates by a least-squares fitting of a 

hyperbola to the traveltime data and the Dix formula. The 

result for crt = 16 ms is shown in table 5 and shows that the 

result by CMP domain inversion is better than Dix' s result 

for shallow reflections. This is possibly due to the 

hyperbolic assumption of traveltime. The interval velocity 

estimation by the Dix formula is also very insensitive to 

measurement noise. 

The CMP domain inversion result for figure 6B 

shows that the layer parameters ·of the. low-velocity thin 

layer have more error than those of the high-velocity thin 

layer. For example, when O't = 16 ms, the interval velocity is 

12 Linearized Inversion of Reflection Traveltimes 

-0.72 0.48 

-5.98 1.35 

-0.46 2.08 

-0.61 2.18 

-8.56 3.01 

16.22 5.07 

-7.83 8.00 

estimated to be 1.88 km/s, which is about a 33 percent error; 

thickness is estimated to be 60 m, which is 40 percent less 

than the true thickness. On the other hand, for the highvelocity layer, the velocity error is less than 4 percent and 

the thickness error is about 20 percent. When using Dix 's 

formula for the low-velocity layer, the velocity is estimated 

to be 1.03 km/s, an error of more than 100 percent, and 

thickness is estimated to be less than 60 percent of the true 

thickness. This result is inferior to the CMP domain 

inversion result. However, in the shot domain, inversion 

errors for the thin layer were similar to each other regardless 

of velocity. The reason for the greater error in low-velocity, 

thin-layer inversion can be understood by examining the 

ray-path diagram shown in figure 6. The rays are focused 

around the CMP point for the low-velocity layer (fig. 6B) 

while the rays spread out farther laterally for the highvelocity layer (fig. 6A). Therefore, the arrival time 

difference between near-offset and far-offset for the lowvelocity thin bed is much smaller than that for the high-
OFFSET, IN KILOMETERS 

-3 -2 -1 0 2 3 -3 -2 -1 0 2 3 

0--+----~ I 

2 2 

(j) 

0: 

UJ 

f--

UJ 

~ 

0 

...J 

~ 4 4 

~ 

I. 

f--

a.. 

UJ 

0 

6 6 

A B 

8~--------------------------------------~ -L--------------------------------------~8 

Figure 6. Ray tracings for common midpoint (CMP) gathers for CMP domain inversion. The CMP interval is 0.1 km, and 

model parameters are shown in table 1. A, High-velocity thin layer (velocity (v) = 4.8 km/s); 8, Low-velocity thin layer (v = 2.8 

km/s). 

velocity thin bed. In fact, the arrival time difference is 178 

ms for the low-velocity layer and 532 ms for the highvelocity layer. Thus, it may be more difficult to resolve the 

parameters of the low-velocity thin layer using CMP 

domain geometry. 

MODE CONVERSION 

Traveltime inversion including mode conversion at 

the reflecting boundary can be easily implemented using 

different velocities for the upgoing and downgoing ray 

paths. However, we cannot estimate all four layerparameters (P-wave velocity, S-wave velocity, dip, and 

depth) simultaneously: we must know either the P-wave or 

S-wave velocity. This can be accomplished by first using 

P-P inversion (or S-S inversion), and then using the 

estimated P-wave (or S-wave) velocity to estimate other 

layer parameters using arrival times of the converted wave. 

The ray-path diagram for the converted waves is 

shown in figure 7 using model 1. Figure 7 A shows ray paths 

for the P-S converted waves and figure 7 B shows ray paths 

for the S-P converted waves. The ray paths shown in figure 

7 indicate that the rays of P-S converted waves spread out 

farther laterally than the S-P converted waves. 

Inversion results for converted waves with crt = 8 ms 

are shown in table 6. Both estimated P-velocity (or S-velocity) using P-P inversion (or S-S inversion) and true 

velocity were used for comparison. Overall, there is no 

significant difference in the inversion result whether we use 

estimated or true velocity. Also, there is no significant 

difference in the inversion result between P-S and S-P 

converted waves. For layer 7, however, table 6 indicates that 

S-P data is slightly better than P-S data at estimating target 

depth. Because of the broad lateral coverage of rays from 

the P-S data shown in figure 7, intuition leads one to 

suspect that P-S data performs better than S-P data when 

estimating layer parameters. However, traveltime data show 

that the arrival time difference between near- and far-offsets 

is greater in S-P data than in P-S data. For this reason, the 

S-P data provides better results when estimating parameters 

of the deeper reflectors in this example. 

These results demonstrate that the estimation of layer 

parameters using converted waves is feasible in principle. 

However, a multicomponent seismic data set using threecomponent geophones and different kinds of source 

polarization (Squire and others, 1988) is needed to utilize 

converted waves for traveltime inversion. 

Mode Conversion 13 
Table 5. Inversion results and Dix result for common midpoint shooting geometry 

and variable amounts of noise contamination 

[Shooting geometry shown in figure 6. Common midpoint (CMP) interval equals 100 meters; 20 

channels] 

Velocity 

Layer (kilometers per second) Dip Depth 

number P-wave (degrees) (kilometers) 

Inversion result; random noise rontent (ua) = 0 millisuonds 

1 1.51 0.0 0.70 

2 2.55 0.0 1.32 

3 3.84 0.0 2.03 

4 4.76 0.0 2.14 

5 3.98 0.0 3.03 

6 4.39 0.0 5.02 

7 3.37 0.0 7.05 

Inversion result; random noise rontent (ua) = 4 millisuonds 

1 1.51 0.0 0.70 

2 2.54 0.0 1.32 

3 3.84 0.0 2.03 

4 4.76 0.0 2.12 

5 4.00 0.0 3.03 

6 4.34 0.0 5.00 

7 3.37 0.0 7.03 

Inversion result; random noise rontent (ua) = 16 milliseronds 

1 1.51 0.0 0.69 

2 2.53 0.0 1.31 

3 3.84 0.0 2.02 

4 4.65 0.0 2.10 

5 4.02 0.0 3.03 

6 4.19 0.0 4.92 

7 3.38 0.0 7.01 

Inversion result; random noise tontent (u1

) = 32 millisuonds 

1 1.50 0.0 0.68 

2 2.52 0.0 1.30 

3 3.83 0.0 2.01 

4 3.83 0.0 2.07 

5 4.08 0.0 3.02 

6 4.00 0.0 4.82 

7 3.39 0.0 6.86 

Dix result; random noise tontent (u1

) = 16 milliseronds 

1 1.51 0.0 0.70 

2 2.77 0.0 1.39 

3 4.10 0.0 2.13 

4 4.88 0.0 2.22 

5 3.90 0.0 3.10 

6 4.15 0.0 4.97 

7 3.37 0.0 7.00 

14 Linearized Inversion of Reflection Traveltimes 
OFFSET, IN KILOMETERS 

-3 -2 -1 0 2 3 

0-r------~~~~~~ 

2 2 

(f) 

a: 

w 

!--

w 

~ 

0 

....J 

:;;2 4 4 ~ 

:i 

!--

a.. 

w 

Cl 

6 6 

A 

8~--------------------------------------J_ 

I B ~---------------------------------------~8 

Figure 7. Ray tracings for shot domain inversion including mode conversion at the reflecting boundary. The geophone group 

interval is 0.2 km, and model parameters are shown in table 1. A, Longitudinal-shear (P-S) conversion; B, Shear-longitudinal (S-P) 

conversion. 

REAL DATA EXAMPLES 

In an attempt to apply the IWI inversion method to a 

real data set, a shot and CMP gather using real marine 

seismic data were used~ 

Figure 8A shows the observed arrival times (squares) 

and the computed arrival times (crosses) based on the shot 

domain inversion results shown in table 7. The RMS error 

between the observed and computed arrival times is about 

1.6 ms. The corresponding ray-path diagram (not shown) 

for the inversion results indicates that there are crossing 

layers within the ray paths (layers 2 and 3 cross at a depth 

of about 1.5 km). Currently, the inversion program cannot 

handle crossing layers within a ray path. Therefore, offsets 

greater than 1.85 km were omitted for the deeper layer(s). 

The results produced by CMP domain inversion are 

shown in figure 8B and table 7. The RMS error between 

observed and computed arrival times is about 2.2 ms. Also 

included are the results by fitting the arrival times by a 

hyperbola using the least squares method (table 7, fig. 8C); 

the RMS error between observed and computed arrival 

times is about 3.7 ms. This large RMS error was introduced 

by the mismatch between computed and measured arrival 

times at farther offsets. Notice that the computed arrival 

times are systematically earlier for the far offsets of layer 5 

in figure 8C. 

The inversion result in the CMP domain is very 

similar to the conventional method based on hyperbolic 

moveout even though a slight advantage of the inversion 

method can be justified based on figures 8B and 8C. 

Shot domain and CMP domain results varied greatly. 

During the inversion of the CMP gather, the process 

converged to consistent estimates over wide ranges of initial 

conditions, regularization, and ~·s (see equation 10), while 

shot domain inversion provided many inconsistent results. 

For example, in the range of RMS error = 1. 7 ms, the 

velocity of the second layer was estimated by shot domain 

inversion at 0.74 km/s with different regularization. 

The problem in the shot domain inversion was in the 

second layer. The resolution matrix (see appendix A) of the 

second layer for the result shown for the shot domain 

inversion in table 7 is: 

[o.m 0.08 0.11] 

= 0.08 0.99 0.0 

0.11 0.0 0.99 

Real Data Examples 15 
Table 6. Inversion results including mode conversions for split-spread shooting geometry 

[Shooting geometry shown in figure 7. Group interval equals 200 meters; near-offset distance equals 152 meters; 20 channels] 

Velocity 

Layer {kilometers ~er second} Dip Depth 

number P-wave 8-wave (degrees) (kilometers) 

P-S tonversion used estimated P-wave velocity 

1 1.52 1.00 7.10 0.71 

2 2.51 1.25 -4.89 1.30 

3 3.88 1.92 0.03 2.01 

4 4.91 2.57 0.41 2.12 

5 3.63 1.94 -9.50 2.98 

6 4.34 I 

1.96 16.49 5.06 

7 5.04 2.00 -5.31 7.67 

P-S tonversion used true P-wave velocity 

1 1.50 1.02 7.32 0.71 

2 2.50 1.26 -4.87 1.31 

3 3.80 1.97 0.22 2.03 

4 4.80 2.65 0.61 2.13 

5 3.80 1.87 -9.19 2.98 

6 4.20 1.93 16.43 5.03 

7 3.50 2.40 -4.38 7.62 

S-P ronversion used estimated S-wave velocity 

1 1.52 1.01 7.25 0.71 

2 2.51 1.25 -4.85 1.31 

3 3.85 1.92 -0.16 2.01 

4 4.98 2.54 -0.70 2.12 

5 3.64 1.95 -9.55 2.98 

6 4.31 1.94 16.59 5.05 

7 4.90 1.94 -3.53 7.57 

S-P tonversion used true S-wave velocity 

1 1.52 

2 2.51 

3 3.88 

4 5.03 

5 3.58 

6 4.29 

7 5.52 

As mentioned by Tarantola (1987) and Jackson 

~1972), the resolution matrix of well-resolved parameters is 

close to the identity matrix. The resolution matrix of layer 2 

has a very low number for the slowness. This means that we 

may not resolve the velocity of the second layer with the 

current data set and the assumptions of our current model 

(i.e., assuming that the earth's crust consists of iso-velocity 

layers separated by plane interfaces). 

Although dip directions estimated by shot domain 

inversions are consistent with the seismic data, this example 

shows that shot domain inversion has many difficulties. 

16 Linearized Inversion of Reflection Traveltimes 

1.00 7.35 0.71 

1.25 -4.94 1.30 

1.90 0.13 2.01 

2.50 0.73 2.11 

2.00 -9.59 2.99 

1.90 16.28 5.02 

1.60 -2.77 7.28 

As mentioned earlier, the main advantage of shot 

domain inversion over the CMP domain inversion is its 

ability to estimate dips. Both shot domain and CMP domain 

inversion were applied to a 1 024-channel sign-bit data set. 

The migrated portion of the line is shown in figure 9, and 

the inversion analysis was done in the vicinity of CMP 

1650. The measured arrival times are shown as squares and 

the computed arrival times are shown as crosses in figure 

10. The average RMS error between observed and 

computed arrival times was 14 ms; most of the errors 

occurred at layer 6. 
(/) 

0 

z 

0 

u 

L.U 

(/) 

~ 

w 

~ 

f= 

Table 7. Comparison of inversion results using real marine seismic data 

[RMS, root-mean-square; CMP, common midpoint] 

Layer 

number 

Velocity 

(kilometers per second) 

P-wave 

Dip 

(degrees) 

Shot domain inversion; RMS error = 1.6 milliS«onds 

1 1.47 -3.34 

2 1.12 -3.00 

3 1.80 -5.96 

4 1.24 -6.44 

5 2.40 0.86 

CMP domain inversion; RMS error = 2.2 milliS«onds 

1 1.52 0.0 

2 1.21 0.0 

3 2.06 0.0 

4 1.66 0.0 

5 2.43 0.0 

Depth 

(kilometers) 

0.73 

0.81 

1.01 

155 

0.62 

0.74 

0.80 

1.08 

1.65 

CMP domain hyperbolk fitting; RMS error = 3. 7 milliS«onds 

1 1.51 0.0 0.61 

2 1.22 0.0 0.74 

3 2.19 0.0 0.81 

4 1.64 0.0 1.08 

5 2.50 0.0 1.68 

OFFSET, IN KILOMETERS 

0.500,_-------,.-------_,2 ________ ~3 o.--------,.-------_,2 __________ 3 or--------,.-------_,2 ________ ~3 

1.00 

1.50 

2.00 

A 0 0 

B 

Figure 8. Graphs showing measured and computed arrival times for marine multichannel data. Squares indicate measured arrival 

times and crosses indicate estimated arrival times. A, Shot domain inversion results with root-mean-square (RMS) error of 1.6 ms; 

8, Common midpoint (CMP) domain inversion results with RMS error of 2.2 ms; C, Interval velocity estimated assuming hyperbolic 

moveout times and Dix's formula with RMS error of 3.7 ms. 

Real Data Examples 17 
-l 

= 

r-

:r 

~ 

~ 

N' 

~ 

Q. 

S" 

< 

~ 

... 

(I) 

6' 

= 0 

-;11::1 

~ 

:=! 

~ 

0. 

6' 

= -1 

... 

Ill 

< 

~ 

3' 

~ (I) 

(/) 

0 

z 

0 

u 

L.U 

(/) 

z 

w 

~ 

f= 

1900 

0 

2 

3 

CMP 

1800 1700 1600 1500 1400 

Figure 9. Migrated stacked section of sign-bit data. Estimated interval velocities and dips by a shot domain inversion are shown near common midpoint (CMP) 1650. 
The overall dips of the inversion result of a single 

shot agree well with the migrated section except for layer 2. 

The dip of layer 2 in the migrated section is similar to the 

dip of layer 3, but the inversion result indicates negligible 

dip for layer 2. This is possibly due to the error in picking 

arrival times for the inversion. The interval velocity based 

on the stacking velocity is about 40 to 60 percent higher 

than the inversion result shown in figure 9. However, the 

migration velocity analysis agrees well with the inversion 

result. 

DISCUSSION 

We have described a linear traveltime inversion 

method which is computationally efficient and accurate. 

This computational method was derived by suitable 

estimation of model parameters in such a way that the 

Jacobian matrix is in a lower triangular form. However, the 

application of this method to real data suggests some 

problems. 

The chief problem may be caused by the iso-velocity 

and plane-layer assumptions in our crustal model. 

Particularly in the shot domain, estimation of dip angle is 

closely related to other layer parameters such as velocity 

and depth. The elements of the correlation coefficient 

matrix (Tarantola, 1987) are close to 1.0: this indicates that 

the error in one parameter estimate will always affect the 

other layer parameters. Numerous tests on models and real 

data indicate that the dip angle is most difficult to handle in 

the shot domain inversion. With a given field geometry, 

slight changes in dip angles significantly change the other 

layer parameters. 

In the shot domain, dip estimation may be affected 

significantly due to the lateral heterogeneity of the media. 

Thus, with real data, the result of the shot domain inversion 

frequently shows inconsistent estimates. In order to deal 

with this difficulty, we need a more complex model than the 

iso-velocity and plane-layer model considered in this paper. 

For the case of a more complex crustal model, the Jacobian 

matrix cannot be computed as shown in equation 7 but can 

be computed by a finite-difference method. The variablemode velocity assumption by Zhu and Brown (1987) is one 

way of introducing lateral heterogeneity. Even in this case, 

the iterative whole-layer inversion method described here 

can be implemented by the parameterization shown in this 

paper. 

In the CMP domain, because we ignored dip angle 

during inversion, the estimated layer parameters were more 

consistent than the layer parameters estimated in the shot 

domain. As mentioned by Diebold and Stoffa (1981 ), the 

CMP domain inversion averages out lateral heterogeneity, 

and, as a consequence, the iso-velocity assumption seems to 

work better in this domain. If we want to investigate a 

detailed velocity structure of the lower crust where the 

two-way traveltime is in the range of 5 to 10 s, an 

expanding-spread geometry could have an advantage over a 

long-offset shot gather. 

Problems may arise in picking arrival times in the 

shot domain. The error in dip estimation for layer 2, using 

the sign-bit data shown in figure 9, was caused by picking 

a crossing event. The original data showed many crossing 

events in the upper part of the section making it difficult to 

identify true reflections. In cases where there is much 

coherent noise, extensive shot domain pre-processing is 

necessary. Therefore, in such cases, implementing 

conventional velocity analysis based on hyperbolic moveout 

is much easier. 

The RMS error between observed and computed 

arrival times serves as an indicator of convergence. 

However, because of the non-uniqueness of the inversion, 

sound geological judgment should be employed when 

interpreting inversion results. As indicated in the model 

study, convergence depends on the initial estimate as well as 

on the inversion method used. For example, the RMS error 

of one of the inversion results for the sign-bit data shown in 

figure 10 was 7.3 ms, which is about one-half of the RMS 

error of the result shown in figure 9. The inversion result 

with an RMS error of 7.3 ms was identical to the result 

shown in figure 9 except for layer 6. The velocity of layer 

6 was estimated as 9.8 kft/s (kilofeet per second) with an 8° 

dip. Both inversion results reveal a low-velocity zone for 

layer 6. Based on a priori geologic knowledge of the area, 

we chose the 13.5 kft/s velocity layer over the 9.8 kft/s layer 

even though it had a higher RMS error value. In this 

example, because we do not have other independent sources 

of velocity information such as well logs, we are not sure 

which value is closer to the true interval velocity. 

CONCLUSIONS 

Based on model analysis and real data examples, the 

following conclusions can be drawn: 

1. The iterative whole-layer inversion method is 

superior to the layer-stripping method for estimating 

layer parameters such as velocity, dip, and depth in 

the shot domain. However, in the CMP domain, the 

two methods yield similar results. 

2. A split-spread geometry gives better results than 

off-end geometry for the shot domain inversion, 

particularly for dip estimation. 

3. CMP domain inversion is less sensitive to 

measurement noise compared to shot domain 

inversion. 

4. In order to calculate a reliable layer parameter, the 

surface-spread length should approximately equal the 

target depth of interest. 

Conclusions 19 
OFFSET, IN THOUSANDS OF FEET 

-15 -10 -5 0 5 10 15 

0.0 

I I I I I 

ifl Ill~ l'l®fillil 

~, 

0.50 t-- l'l ~ ® @ fil 

l'l l'l 

1il l!l® 

1.00 1-- 19 j;l til lil til !! IS l!l l'l I'! lill!l 

(iJ II lil til 19 tillol lil lil fil (/) til 0 15 z 

0 

u 

~ 1.50 t-- ~ ~ R j;j !l l'l ~ 1!1 ~ ~ ~ !! til u..i ~ ~ :2: 

f= 

2.00 t--

2.50-

3.00~------------------------------------------------------------------------~ 

Figure 10. Graph of measured arrival times (squares) and computed arrival times (crosses) resulting from a shot domain 

inversion near common midpoint (CMP) 1650 in figure 9. The root-mean-square (RMS) error of inversion is about 14 ms. 

5. Regularization is a good method of reducing the 

variance of layer parameters, but, owing to the 

inconsistency of the solution, care should be taken 

when implementing regularization. 

6. CMP domain inversion has an advantage over 

conventional velocity analysis based on hyperbolic 

moveout for longer offset data. The conventional 

method has an advantage in implementation because 

it does not require the digitization of arrival times. 

7. The inversion of converted waves may be 

feasible using a real data set. 

20 Linearized Inversion of Reflection Traveltimes 

APPENDIX A 

Resolution Matrix 

The resolution matrix R can be defined in a singular 

value decomposition of the Jacobian matrix (equation 4) by 

the following formula (Jackson, 1972): 
Because V is an orthonormal eigenvector, R is the 

identity matrix (1). In the case when solutions corresponding 

to smaller singular values are eliminated in order to reduce 

the variance of the solution, the resolution matrix can be 

written as: 

q 

R ij = L vik vjk 

k=I 

(Al) 

where 

q represents the number of singular values kept for 

the solution. 

The resolution matrix given in the main text was 

evaluated using equation A 1 with q = 2. 

The physical meaning of the resolution matrix can be 

illustrated by a formula given by Tarantola (1987). He 

presented the following formula for the resolution matrix: 

where 

(A2) 

<M> is the inverse solution of ~odel parameters, 

Mt is the true model parameter, and 

MP is the prior model parameter. 

Equation A2 indicates that we are able to see the truth only 

through a filter R. In a favorable case, R~ I. 

Thus, 

The more R differs from the identity matrix, the more 

the estimation differs from the true model parameters. 

REFERENCES CITED 

Aid, K., and Lee, W.H.K., 1976, Determination of threedimensional velocity anomalies under a seismic array using 

first P arrival times from local earthquakes; 1, A 

homogeneous initial model: Journal of Geophysical 

Research, v. 81, no. 23, p. 4381-4399. 

Aki, K., and Richards, P.G., 1980, Quantitative seismologytheory and methods, 2: San Francisco, W.H. Freeman, 932 

p. 

Backus, G., and Gilbert, F., 1970, Uniqueness in the inversion of 

inaccurate gross Earth data: Philosophical Transactions of 

the Royal Society of London, v. 266, no. 1173, p. 123-192. 

*U.S. GOVERNMENT PRINTING OFFICE: 1991-573-049/46004 

Bishop, T.N., Bube, K.P., Cutler, R.T., Lanegan, R.T., Love, P.L., 

Resnick, J.R., Shuey, R.T., Spindler, D.A., and Wyld, H.W., 

1985, Tomographic determination of velocity and depth in 

laterally varying media: Geophysics, v. 50, no. 6, 

p. 903-923. 

Chiu, S.K.L., and Stewart, D.R., 1987, Tomographic 

determination of three-dimensional seismic velocity 

structure using well logs, vertical seismic profiles, and 

surface seismic data: Geophysics, v. 52, no. 8, 

p. 1085-1098. 

Diebold, J.B., and Stoffa, P.L., 1981, The traveltime equation, 

tau-p mapping, and inversion of common midpoint data: 

Geophysics, v. 46, no. 3, p. 238-254. 

Diebold, J.B., Stoffa, P.L., and the LASE study group, 1988, A 

large aperture seismic experiment in the Baltimore Canyon 

trough in Sheridan, R.E., and Grow, J.A., eds., The Atlantic 

continental margin; U.S., The geology of North America, v. 

1-2: The Geological Society of America, p. 387-398. 

Dix, C.H., 1955, Seismic velocity from surface measurements: 

Geophysics, v. 20, no. 1, p. 68-86. 

Gj!Z)ystdal, Havar, and Ursin, Bj!Z)rn, 1981, Inversion of reflection 

times in three dimensions: Geophysics, v. 46, no. 7, p. 

972-983. 

Jackson, D.D., 1972, Interpretation of inaccurate, insufficient, and 

inconsistent data: Geophysical Journal of the Royal 

Astronomical Society, v. 28, p. 97-109. 

Lee, M.W., 1990, Traveltime inversion using transmitted waves of 

offset VSP data: Geophysics, v. 55, no. 8, p. 1089-1097. 

Lines, L.R., ed., 1988, Inversion of geophysical data: Geophysics 

reprint series 9, Society of Exploration Geophysicists, 

543 p. 

Lines, L.R., Bourgeois, A., and Covey, J.D., 1984, Traveltime 

inversion of offset vertical seismic profiles-a feasibility 

study: Geophysics, v. 49, no. 3, p. 250-264. 

Lines, L.R., and Treitel, S., 1984, Tutorial: A review of leastsquares inversion and its application to geophysical 

problems: Geophysical Prospecting, v. 32, no. 2, 

p. 159-186. 

Liu, Char-Shine; Zhu, Tianfei; Farmer, Harlow; and Brown, Larry, 

1986, An expanding-spread experiment during COCORP's 

field operation in Utah, in Barazangi, M., and Brown, L., 

eds., Reflection seismology: A global perspective: 

American Geophysical Union, Geodynamics Series, v. 13, 

p. 237-246. 

Sluis, A. van der, and Vorst, H.A. van der, 1987, Numerical 

solution of large, sparse linear algebraic systems arising 

from tomographic problems in Nolet, G., ed., Seismic 

tomography with applications in global seismology and 

exploration geophysics: Dordrecht, D. Reidel Publishing 

Company, p. 49-83. 

Squire, S.G., Kim, C.D.Y., and Kim, D.Y., 1988, Interpretation of 

total wave field data over Last Hills Fields, Kern County, 

California: 58th Annual International Meeting of Society of 

Exploration Geophysicists, Expanded Abstracts, p. 

1137-1141. 

Tarantola, A.T., 1987, Inverse problem theory: Amsterdam, 

Elsevier, 613 p. 

Zhu, T., and Brown, L.D., 1987, Two-dimensional velocity 

inversion and synthetic seismogram computation: Geophysics, v. 52, no. 1, p. 37-50. 

References Cited 21 



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 hydrogeology, 

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 tn- 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 geophysical 

techniques, such as gravity, magnetic, seismic, or radioactivity, which 

reflect subsurface structures that are of economic or geologic significance. 

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 tn.-minute quadrangle photo geologic maps on planimetric bases 

which show geology as interpreted from aerial photographs. 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-andwhitemaps on topographic or planimetric bases presenting a wide range 

of geohydrologic data of both regular and irregular areas; 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 the U.S. Geological Survey, Books 

and Open-File Reports Section, Federal Center, Box 25425, 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 of 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-SupplyReportsandMapsFor(State),"maybepurchasedbymail 

and over the counter in paperback booklet form only 

"Price and A vallabillty 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" available free of charge by mail or may be obtained 

over the counter in paperback booklet form only. lnose 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.