Geophys. J. Int. (2006) 167, 1373–1384 doi: 10.1111/j.1365-246X.2006.03162.

Time domain Gauss–Newton seismic waveform inversion

in elastic media

Dong-Hoon Sheen,1,2 Kagan Tuncay,1 Chang-Eob Baag2 and Peter J. Ortoleva1

1 Laboratoryfor Computational Geodynamics, Indiana University, Bloomington, Indiana, 47405, USA
2 Schoolof Earth and Environmental Sciences, Seoul National University, San 56–1, Shillim-dong, Gwanak-gu, Seoul, 151–742 Korea.
E-mail: [email protected]

Accepted 2006 August 3. Received 2006 August 3; in original form 2005 March 30

SUMMARY
We present a seismic waveform inversion methodology based on the Gauss–Newton method
from pre-stack seismic data. The inversion employs a staggered-grid finite difference solution
of the 2-D elastic wave equation in the time domain, allowing accurate simulation of all
possible waves in elastic media. The partial derivatives for the Gauss–Newton method are
obtained from the differential equation of the wave equation in terms of model parameters.
The resulting wave equation and virtual sources from the reciprocity principle allow us to apply
the Gauss–Newton method to seismic waveform inversion. The partial derivative wavefields are
explicitly computed by convolution of forward wavefields propagated from each source with
reciprocal wavefields from each receiver. The Gauss–Newton method for seismic waveform
inversion was proposed in the 1980s but has rarely been studied. Extensive computational
and memory requirements have been principal difficulties which are addressed in this work.
We used different sizes of grids for the inversion, temporal windowing, approximation of
virtual sources, and parallelizing computations. With numerical experiments, we show that the

GJI Seismology
Gauss–Newton method has significantly higher resolving power and convergence rate over the
gradient method, and demonstrate potential applications to real seismic data.
Key words: finite difference methods, Frechét derivatives, synthetic seismograms, wave
equation, waveform inversion.

partial derivative seismic wavefields which had been used for elec-
1 I N T RO D U C T I O N
tromagnetic problem (Rodi 1976), and solved seismic waveform in-
Seismic waveform inversion can be defined as an iterative proce- version with the Newton method. The partial derivative wavefields
dure for obtaining physical properties of the Earth from pre-stack are obtained from new wave propagation simulations driven by the
seismic data. It is well known that the inversion of seismic data is virtual sources at the grid points where model parameters are deter-
a computationally demanding task. Early studies include those of mined. In other words, the number of required forward simulations
Lailly (1983) and Tarantola (1984) who presented a method of con- is equal to the number of model parameters.
structing the gradient (or steepest descent) direction for the inversion Shin et al. (2001) used an efficient method for calculating partial
of the acoustic problem without computing the partial derivatives derivative wavefields using the reciprocity relation between the vir-
explicitly. This method constructs the gradient direction by cross- tual sources and the receivers. The reciprocity theorem is proven
correlating forward propagated wavefields from a seismic source in Aki & Richards (1980) for an elastic anisotropic continuous
with backward propagated wavefields from the data residuals. As medium. This theorem allows the source–receiver locations to be
each iterative loop of the inversion requires only several forward interchanged. The recorded seismograms are identical if the sources
simulations for each seismic source, it made seismic waveform in- and receivers are located inside the domain or on its boundary
version feasible in the 1980s. Mora (1987) applied this method to (Eisner & Clayton 2001). Thus, the computation of the partial deriva-
elastic problems in the time domain, whereas Pratt (1990) and Pratt tive wavefields doesn’t depend on the number of model parameters
et al. (1998) applied it to elastic and acoustic problems, respectively, but depends on the numbers of shots and receivers.
in the frequency domain. Although developments in computer technology have been im-
Tarantola (1984) also presented the Gauss–Newton algorithm pressive, it has been impractical to make use of a Newton type
called ‘total inversion’ (Tarantola & Valette 1982), although it was method for high-resolution seismic inversions. To avoid extremely
not possible to implement when it was presented because of limited expensive computation of the Jacobian or the Hessian matrix, Hicks
computational resources. In recent years, however, it has become & Pratt (2001) proposed a two-step inversion procedure. The back-
feasible. Pratt et al. (1998) used ‘virtual source’ terms to obtain ward propagation method is used for finding reflectors consisting


C 2006 The Authors 1373
Journal compilation 
C 2006 RAS
1374 D.-H. Sheen et al.

of a large number of parameters, and then the Newton method is how the Gauss–Newton method outperforms the gradient method.
applied to background velocities with a much smaller number of We conclude by showing several numerical examples illustrating
parameters. Shin et al. (2001) took advantage of the diagonally the application of seismic waveform inversion.
dominant nature of the ‘approximate’ Hessian matrix (Pratt et al.
1998). Diagonal elements of the approximate Hessian were used as
a pre-conditioner for an iterative inversion. In all of these studies, 2 I N V E R S E P RO B L E M
the finite element method was used to solve the forward acoustic Seismic waveform inversion is the problem of finding properties
problems. of the Earth from seismic data. In order to infer a set of model
Pratt et al. (1998) and Shin et al. (2001) used a frequency domain parameters which represent the Earth, the inverse problem seeks
method to carry out the inversion with discrete frequencies. By to minimize the residuals between the model response obtained by
selecting a few frequencies, the method allows a significant re- forward simulation and the observed seismic data.
duction in computational burden. However, success of the method In general, the seismic responses d of the Earth represented by
requires a careful selection of frequencies, which highlights the ro- model parameters m would be recorded at receivers. This relation-
bustness of the time domain approach over the frequency domain ship can be expressed with the non-linear functional F:
approach (Freudenreich & Shipp 2000). Temporal windowing in the
time domain approach, like discrete frequencies in the frequency d = F(m). (1)
domain approach, is proven to have an important role in reducing The residual error is defined as the difference between the model
computation (Shipp & Singh 2002). responses and the observed data:
In this work, we apply the Gauss–Newton method to elastic wave- 4
form inversion in the time domain. The Gauss–Newton method is d = F(m) − dobs , (2)
known to show local quadratic convergence. Local minimization obs
where d is the observed data. We now introduce the least-squares
algorithm, such as a gradient or Newton method, requires us to set a
problem:
initial model close to a true model. Such inversions, therefore, may
converge to a local minimum (Gauthier et al. 1986; Mora 1987). 1 1
Sd (m) = d2 = dt d, (3)
Several approaches have been proposed to overcome this difficulty. 2 2
Traveltime analysis precedes and gives a initial background veloc- where S d means the data misfit function, the factor 1/2 allows sub-
ity model to the inversion (Shipp & Singh 2002). Decoupling the sequent simplifications, and the superscript t represents the matrix
high and low wavenumber of velocity variation has also been used transpose. The inverse problem becomes the minimization of S d .
to make the solutions to converge to the correct model (Snieder Thus, our purpose is to find a model m∗ which minimize S d (m∗ ).
et al. 1989; Symes & Carazzone 1991; Clément et al. 2001; Hicks
& Pratt 2001). They used different types of the misfit function to
avoid falling into local minima: the simplex method, the differen- 2.1 Gradient method
tial semblance optimization, the migration based traveltime, and The gradient method minimizes S d by updating model parameters
the adaptive depth stretching. Following Shipp & Singh (2002) and in the opposite direction of the gradient of S d (m) iteratively:
Sirgue & Pratt (2004), we assume that the initial model is close
to the global minimum, which can be obtained from the traveltime mn+1 = mn − α∇ Sd , (4)
analysis. where the superscripts represent iteration numbers and α is a step
Because seismic waveform inversion is the systematic fitting of length. The gradient direction can be obtained by taking partial
synthetic to observed seismograms, it is important to generate ac- derivatives of eq. (3) with respect to model parameters m:
curate seismogram which can account for subsurface elastic fea-  
tures. Thus, we use the velocity-stress staggered-grid finite differ- ∂Fi (m) t
∇ Sd = d = Jt d,
ence method (Levander 1988) to solve the elastic wave equation in ∂m j
the time domain and implement the perfectly matched layer (PML) (i = 1, . . . , NS × NR; j = 1, . . . , M), (5)
method (Berenger 1994) as an absorbing boundary condition. The
solution of the elastic wave equation allows us to use multicompo- where Jt is a transpose of the Jacobian matrix, the subscripts i and j
nent data. The major obstacle to seismic waveform inversion is the indicate indices for seismograms and model parameter, respectively,
explicit calculation of the Jacobian and the approximate Hessian ma- and NS, NR and M are the numbers of shots, receivers and model
trices. To overcome this, we utilize the reciprocity principle and the parameters, respectively.
convolution theorem. A Newton type seismic waveform inversion It is not necessary to compute the Jacobian explicitly so as to find
still requires huge amounts of memory and computation. For elastic the gradient direction. The backward propagation method requires
problems, much more resources are required than those for acous- only several forward computations to construct the gradient direc-
tics. We surmount the limitation by using different sizes of grids tion. Formal derivations for the backward propagation method are
in the spatial and time domains for the inversion stage, temporal given by Tarantola (1984) and Mora (1987).
windowing, approximation of virtual sources, and parallelizing the
method for massively parallel computers using the message pass-
2.2 Newton method & Gauss–Newton method
ing interface (MPI) approach. Integration of all these features make
it possible to implement the classical seismic waveform inversion It is well recognized that the Newton and Gauss–Newton meth-
scheme presented in Tarantola (1984). ods are effective and robust techniques for numerical optimization
This paper is organized as follows: we recap the inverse problem of non-linear problems and guarantee faster convergence rates than
and show how to use the reciprocity principle and the convolution the gradient method when the initial model is close enough to a local
theorem for calculating partial derivatives explicitly. We then present or the global minimum. In geophysical inverse problems, however,
the more detailed inversion scheme used in this work and illustrate the Newton method has not been often used because it is difficult to


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
 Gauss–Newton seismic waveform inversion 1375

compute the Hessian matrix and the cost usually surpasses the gain. A choice of λ value between 0 and infinity produces a compromise
The Hessian matrix is the analogue of the second-order derivative. result. After several computations, appropriate values of λ are cho-
Pratt et al. (1998) calculated the second-order derivative using the sen as λ 1 = 0.05 and λ 2 = 0.0005. This means that 5 and 0.05 per cent
second order virtual source which requires as many additional for- of the maximum value of diagonal elements of Jt J are added (Shin
ward simulations as the number of model parameters. In this work, et al. 2001), respectively.
we use the Gauss–Newton method which neglects the second-order The simplest way of choosing the step length is to take it as a
derivative. constant through all the iterations, which can be obtained by trial and
The Newton method (Tarantola 1987) is given by error. In this work, an optimal value of the step length is determined
by a linearized approach. The optimal value for α n is given by
mn+1 = mn − H−1 Jt d, (6)  
where H is the Hessian matrix. Each element of the Hessian matrix n ∼
[Jt gn ]t F(mn ) − dobs
α = , (14)
can be expressed in differential forms: [Jt gn ]t [Jt gn ]
∂  t  ∂Jt where
H = ∇ 2 Sd (m) = J d = Jt J + d. (7)  −1  
∂m p ∂m p gn = Jt J + λLt L Jt dn − dobs . (15)
Because the second term of eq. (7) is usually small and negligible
For more details, see Gauthier et al. (1986) and Pica et al. (1990).
(Tarantola 1987), we obtain the Gauss–Newton formula
−1
mn+1 = mn − [Jt J] Jt d = mn − [Ha ]−1 Jt d, (8)
3 PA RT I A L D E R I VAT I V E S
where H a is the approximate Hessian matrix. The elements of the ap-
proximate Hessian are obtained from the zero-lag crosscorrelation Partial derivatives of seismic waveform inversion represented by the
of the partial derivative wavefields and, specifically, the diagonal el- Jacobian matrix J can be directly obtained from a residual wave-
ements of the matrix are obtained from the zero-lag auto-correlation field from two forward simulations with and without perturbations
of the derivatives. The partial derivative wavefields are mostly un- of model parameters. For magnetotellurics data, Rodi (1976) showed
correlated with each other in the high-frequency limit. However, the that partial derivatives can be obtained from solving forward prob-
wavefields from adjacent nodes are more or less correlated because lems with virtual sources at the perturbation locations. Pratt et al.
the frequency content is finite. Thus, in general, the approximate (1998) applied this to acoustic waveform inversion. This approach
Hessian matrix will be diagonally dominant due to auto-correlation, is appropriate for the case in which there are many stations and a
and banded due to finite frequencies. It is also known that the inverse few model parameters. Rodi also proposed a more efficient method
approximate Hessian matrix operates like a sharpening or focusing when the number of model parameters exceeds the number of the
filter (Pratt et al. 1998). data using the reciprocity principle. Shin et al. (2001) adopted this
In the application of the Gauss–Newton method to geophysical approach for acoustic waveform inversion. In all previous studies,
inversion, regularization is particularly useful for stabilizing the sys- the forward problems were solved in the frequency domain. In this
tem and incorporating a priori information to the problem (Tarantola paper, we formulate the wave propagation simulation in the time
1987). The regularized misfit function S can be defined as, domain which is based on the finite difference method (Levander
1988), and use the PML boundary condition (Berenger 1994) as an
S(m) = Sd (m) + λSm (m), (9) absorbing boundary condition. The convolution theorem is used to
where S m is the model objective function that contains a priori utilize the virtual source and the reciprocity principle.
information of the model and λ is a scalar value that globally controls
the relative importance of the model objective function S m . The
model objective function can be written as combinations of discrete 3.1 Partial derivative wavefields from virtual sources
linear operator L: In 2-D Cartesian coordinates, for an isotropic, linearly elastic
1 medium, the equations of motion can be written as a set of first-
Sm (m) = Lm2 . (10) order hyperbolic equations:
2
Then the regularized Gauss–Newton formula can be written as, ρ v̇i = τi j, j + Fi ,
 −1
mn+1 = mn − α n Jt J + λLt L Jt d, (11) τ̇i j = μ (vi, j + v j,i ) + λ δi j vk,k + G i j , (16)
where α is a step length. We use the conjugate gradient method to where the symbols are given in Table 1.
invert the regularized approximate Hessian to calculate the model In this formulation, the model parameters will be density, ρ, or
updates mn+1 . If L = I (the identity matrix), eq. (11) yields the one of Lamé’s moduli, λ and μ. In the case of density as the model
damped least-squares method (Levenberg 1944; Marquardt 1963). If
L is a discrete spatial differential operator, the model objective func- Table 1. Definition of symbols.
tion controls the roughness of spatial variations among the model Symbol Definition
parameters. Sasaki (1989) used discrete 2-D Laplacian operator:
v Velocity
Pi m = (m i ) E + (m i )W + (m i ) N + (m i ) S − 4(m i ), τ Stress tensor
(12) ρ Density
λ, μ Lamé’s moduli
where the superscripts E, W , N, and S refer to the four neighbours δi j Kronecker delta,
of the ith model parameter and P i is the ith row of the Laplacian δ i j = 0 for i = j and δ i j = 1 for i = j
operator matrix whose elements are either 1, −4, or 0. In this work, Fi, Gij Body force and traction source
these model objective functions are used simultaneously: ,k Spatial derivative, ∂/∂ xk
 −1 ˙ Temporal derivative, ∂/∂t
mn+1 = mn − α n Jt J + λ1 Pt P + λ2 It I Jt d. (13)


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
1376 D.-H. Sheen et al.

This approach, however, is based on the time domain finite differ-

p q p q ence method and can only use limited frequency band of a source
function. In a sense, thus, it is a bandpass filtered partial derivative
wavefield that we can obtain from this approach. Because the seis-
mic data for digital processing and numerical analysis are discrete
F D F R and they can’t carry information above a threshold frequency, for
example, Nyquist frequency, our limited frequency band of partial
m (a) m (b) derivative wavefields still guarantee the resolution of the problem.
We can decompose the Jacobian matrix J into virtual sources and
reciprocal wavefields:
Figure 1. Schematic diagram for general reciprocal relation for partial
derivative wavefields. Partial derivative wavefields from a virtual source m
J p,q = Fpm ∗ Dqm = Fpm ∗ Rqm , (18)
(a) and a reciprocal source (b). Symbols p, q and m denote the indices for
a source, receiver and model parameter, respectively. F, D and R represent where ∗ means convolution.
the forward wavefield, differential wavefield from the virtual source and It should be noted that the reciprocal wavefields are identical only
reciprocal wavefield from the receiver, respectively. if the orientations of the source and the receiver match (Arntsen &
Carcione 2000; Sheen et al. 2004). Fig. 2 shows how to generate the
parameter, eq. (16) is differentiated with respect to a ρp : partial derivative wavefields with the reciprocity principle. Config-
uration (i) can be decomposed into (ii) and (iii). By the reciprocity
∂ v̇i ∂τi j, j ∂ρ principle, those configurations can be replicated by (iv) and (v). The
ρ = − v̇i ,
∂ρ p ∂ρ p ∂ρ p relations between configurations (iii) and (v) of Fig. 2(a) and be-
 
∂ τ̇i j ∂vi, j ∂v j,i ∂vk,k tween (ii) and (iv) of Fig. 2(b) show that the locations of source and
=μ + + λ δi j , (17)
∂ρ p ∂ρ p ∂ρ p ∂ρ p receiver should be reciprocal as well as the orientations. Thus, for
example, the reciprocal simulation of the horizontal motion from
where −(∂ρ/∂ρ p )v̇i is a virtual source for the partial derivative wave- the source can be generated by summing horizontal and vertical re-
fields ∂v i /∂ρp propagated from the location of a particular model sponses from a horizontal source (see Fig. 2a). This suggests that
parameter, ρp . The virtual source in eq. (17) is obtained from solu- eq. (18) becomes
tions of eq. (16). Partial derivative wavefields with respect to Lamé’s  m

moduli, are given in Appendix. m

J p,q = J p,q , Jm , (19)
x p,q z
Simple Earth model can be represented by three parameters, such
as density and Lamé’s moduli or density, and P and S-wave veloc- where
m

ities. Tarantola (1986) showed that the choice of model parameters J p,q x
= Fpm x ∗ Rqm x x + Fpm z ∗ Rqm x z ,
for elastic inversion of seismic data affects the performance of the
inversion. From numerical experiments, he chose density and P- and and
S-wave impedances as adequate parameters for seismic waveform m

J p,q z
= Fpm x ∗ Rqm z x + Fpm z ∗ Rqm z z ,
inversion. In this paper, we choose density and P and S-wave veloc-
ities as the parameters. Expressions for partial derivatives in terms where | x and | z represent the orientation of recorded wavefields and
of density and velocities, and density and impedances are also given x and z mean the orientation of the source component.
in Appendix.

3.3 Jacobian & approximate Hessian matrix

3.2 Seismic reciprocity principle & convolution
The gradient direction in eq. (5) can be obtained by
Fig. 1 shows how the partial derivative wavefields are generated. The
virtual source m samples the wavefield driven by the source p. Then, Jt d = J p,q
m
⊗ d p,q , (20)
the partial derivative wavefield propagates from the virtual source to
where ⊗ means zero-lag crosscorrelation and dp,q represents the
the receiver q. Using the virtual source, thus, explicit computation
residual at the receiver q of the shot gather with the source p. Before
of the partial derivative wavefields requires as many as the number
crosscorrelating the Jacobian and the residual, the residual should
of model parameters of the forward simulation.
be convolved with the source function used for the reciprocal wave-
Seismic waveform inversion is usually an underdetermined prob-
fields R. The procedure to compute the Jacobian matrix, the gradient
lem which means the number of parameters exceeds the data. A
direction and the approximate Hessian matrices for multi shots is
more efficient way of evaluating partial derivatives is to use the reci-
summarized below:
procity principle. The reciprocity principle (Aki & Richards 1980)
allows the source–receiver locations to be interchanged. Thus, we (i) Simulate all forward wave propagations from the sources and
can get the same partial derivative wavefield at the virtual source m sample forward wavefields at the model parameter m. Compute the
which is propagated from the receiver q with the forward wavefield residuals;
F recorded at the virtual source m. This can be calculated by convo- (ii) Simulate reciprocal wave propagations from the horizontal
lution of two wavefields, F and R, at the virtual source m. Then, the and vertical component of the receiver q, respectively. Sample re-
number of forward simulations for the derivatives doesn’t depend ciprocal wavefields at the model parameter m;
on the number of parameters but depends on the number of shots (iii) Convolve the source function of step (ii) and the residuals
and receivers. from step (i);
In the frequency domain, the partial derivative wavefield can be (iv) Construct the Jacobian matrix consisting of partial derivative
obtained by convolution of the forward wavefield F with the recipro- wavefields obtained from convolving wavefields from steps (i) and
cal wavefield R which is the impulse response of the delta function. (ii);


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
 Gauss–Newton seismic waveform inversion 1377

(a) Horizontal receiver (b) Vertical receiver

(i) (i)

(ii) (iii) (ii) (iii)

(iv) (v) (iv) (v)

Figure 2. Reciprocal reconstructions of horizontal response (a) and vertical response (b). Star and triangle denote seismic source and receiver, respectively.
Black arrows represent source orientations and white arrows receiver orientations, respectively.

(v) Calculate the gradient direction by zero-lag crosscorrelation

of the wavefields from step (iv) and (iii). 3.4 Computational aspects
(vi) Calculate the approximate Hessian matrix by zero-lag cross- Although eq. (11) has been proposed several decades ago, it is still
correlation of the Jacobian matrix from step (iv). a challenging work to solve this problem. The algorithm presented
here is implemented with the MPI on an IBM cluster which contains
The explicit calculation of the Jacobian requires NS + 2NR simula- two PowerPC970 2.2 GHz processors and 2 GB of memory in a node
tions. and a total of 968 processors. Shipp & Singh (2002) parallelized the
The drawback of the Gauss–Newton method, exclusive of the backward propagation method so that each simulation runs on differ-
computational demands, is a huge memory requirement to store ent processors simultaneously. In this work, domain decomposition
the Jacobian and the approximate Hessian matrices. Furthermore, is used for parallelizing the forward and reciprocal simulation. Each
the procedure for constructing the Jacobian already costs huge processor solves the problem within its small subdomain, commu-
amounts of the memory [see step (i), above]. However, it turns out nicates with neighbouring processors to update wavefields at each
that the Jacobian matrix does not have to be computed and stored time step and samples forward or reciprocal wavefields within target
as a whole. Because we can divide the Jacobian matrix into sub- area in which model parameters are updated.
matrices, it is possible to compute the approximate Hessian matrix as The resolution required by an accurate forward simulation and
follows: the resolution that can be achieved by an inversion are quite dif-
⎡ ⎤ ferent. Numerical stability criteria for wave propagation simula-
J1 tions require smaller grid sizes in the spatial and time domains than
 ⎢ ⎥  t
NR
Ha = Jt J = Jt1 · · · JtNR ⎣ ... ⎦ = JpJp. (21) those for the inversion. Thus, we use different sizes of grids in
p the spatial and time domains: a fine grid for forward and recipro-
J NR
cal simulations, and a coarse one for seismic inversion. After both
simulations, a virtual source for the inversion is generated by sum-
Therefore, the full Jacobian matrix does not have to exist at one ming wavefields recorded at grid points inside the block for a model
time. The approximate Hessian matrix can be computed by sum- parameter.
ming a series of elements. Each element in H a is a zero lag cross- Because of block parametrization, partial derivatives of each vir-
correlation of two partial derivatives at the receiver locations. At tual source consist of a large number of partial derivatives inside a
the high-frequency limit, these wavefields would be perfectly un- block. As the number of grid points describing a virtual source in-
correlated with each other. However, the partial derivatives from creases, the accuracy of the partial derivative seismogram increases
adjacent nodes are correlated to some extent at finite frequencies. as well. However, there are no significant discrepancies, although
Thus H a is diagonally dominant, due to the auto-correlation on derivatives at only 4 grid points are used. Therefore, memory and
the diagonal and banded due to the finite frequency effects (Pratt computation for the Jacobian matrix can be reduced by more than
et al. 1998). In addition to this, the approximate Hessian ma- 90 per cent.
trix is symmetric so that it’s not necessary for all elements to be In the spatial domain, a model parameter is taken as a block and
computed. in the time domain, longer sampling interval and time window of
Because the full Jacobian matrix is not used, Jt gn in eq. (14) can wavefields are used. Required time window for wave propagation
be obtained by a finite difference approximation: simulations is determined to receive all possible reflected waves
from the subsurface. The simulation over the entire time window is
1 not necessary for the derivatives but for the residuals because the
Jt gn ∼
= [F(mn + gn ) − F(mn )] , (22) derivatives are obtained from convolution of waves from a source

with those from a receiver. Thus we determine the time window
length (NT p ) for the partial derivative by
where is a sufficiently small value. Then each iteration step of the
inversion requires 2NS + 2NR simulations. N T p = N T /2 + (DT + 2/C F)/dt + TP,


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
1378 D.-H. Sheen et al.

0.0 0.0 0.0

(a) (b) (c)

0.4 0.4 0.4
Depth (km)

0.8 0.8 0.8

1.2 1.2 1.2

Gradient method Gauss-Newton method
True model True model True model
Initial model Iteration 1 Iteration 1
1.6 1.6 1.6
Iteration 100 Iteration 10

3.2 3.6 4.0 4.4 4.8 5.2 3.2 3.6 4.0 4.4 4.8 5.2 3.2 3.6 4.0 4.4 4.8 5.2
P velocity (km s−1) P velocity (km s−1) P velocity (km s−1)
Figure 3. Inversion results from 1-D model: (a) P-wave velocity structures of the true model and the initial model and (b and c) updated velocity structures
from the gradient method (b) and the Gauss–Newton method (c).

where NT and dt are the number of time step, the time increment resolution is decreasing with depth. As a rule, deeper part of the
for wave propagation simulations, respectively, DT and CF repre- model can be resolved after shallower part has been resolved. The
sent a delay time and central frequency of a source time function, low-frequency content of the model is already given by the initial
respectively, and TP is a cosine tapering length. Consequently the re- model. In both inversions, the high-frequency content represent-
ciprocal simulations require about 70 per cent of computation of the ing the discontinuity converges more rapidly than that of the inter-
forward simulations, which reduces the computational and memory mediate frequency representing the layer. As expected, the Gauss–
burden significantly. Newton method shows a good convergence rate and the gradient
method need much more iterations to resolve the deepest part of the
model.
4 N U M E R I C A L R E S U LT S
Pre-conditioning of the gradient can improve the convergence
In this section we present some inversion results of synthetic seismic rate of the gradient method. However, the choice of pre-conditioner
data to examine the resolving power of the Gauss–Newton seismic significantly affects its efficiency. Here, a simple pre-conditioner is
waveform inversion. In all examples, the source is assumed to be considered in which we multiply the gradient by an arbitrary power
known and the synthetic data have been generated from the same of depth. From several numerical experiments, the pre-conditioner
kernel as that used for the inversion. improves the convergence rate as much as about 20 per cent of the
gradient method. However, the Gauss–Newton method still shows
much faster convergence rate than the gradient method.
4.1 Example 1
Fig. 4 shows the seismograms and the residuals. The scale of the
First, we use a simple layered model in which P-wave velocity varies final residuals (Figs 4d and e) is twice as the others. The seismo-
only with depth whereas density and S-wave velocity are assumed to gram from the true model (Fig. 4a) contains significant reflections
be constant and known. P-wave velocity of a initial model linearly and multiples, all of which contribute to the inversion, whereas the
increases with depth (see Fig. 3). Only the vertical components of seismogram from the initial model (Fig. 4b) doesn’t show signifi-
seismogram are used as the observed data. Both source and receiver cant events. Thus, the initial residual (Fig. 4c) shows most of events
are located at the free surface. The synthetic data has 5 shot gathers from the true model. In the final residuals from the gradient method
which are generated by a point body force of the first derivative of (Fig. 4d), there are still intermediate frequency discrepancies gen-
Gaussian function with a dominant frequency of 30 Hz. Each shot erated from the intermediate spatial frequency content at shallow
gather has 20 receivers with a spacing of 100 m. The minimum and and deep region of the model. The excellent recovery by the Gauss–
maximum offsets are 100 and 1800 m, respectively. The grid spacing Newton method is shown in the final residuals (Fig. 4e). All events
and the time step for wave propagation simulation are 4 m and 0.4 are well recovered and, consequently, the final model is very close
ms, respectively, whereas, for the inversion, vertical and horizontal to the true model (see Fig. 3c).
grid spacings are different which are 20 and 40 m, respectively, and The inverse of the approximate Hessian matrix is known as a fo-
the time step is 4 ms. cusing filter or a normalizing pre-conditioner. The artefacts in the
Fig. 3 shows the inversion results from the gradient and the Gauss– model update from the gradient method can be effectively removed
Newton methods. The result from the gradient method is obtained by the Gauss–Newton method (Pratt et al. 1998). Shin et al. (2001)
after 100 iterations and the residual error decreases to 1.0 per cent of showed that a diagonal approximation of the approximate Hessian
its initial value whereas the residual error from the Gauss–Newton normalizes the gradient direction so that it improves a poorly scaled
method, even only after 10 iterations, achieves 0.4 per cent reduc- image. Fig. 5 shows the model updates after the first iteration of
tion of the initial. As shown in Kolb et al. (1986), the model tends the gradient and the Gauss–Newton methods. The update of the
to converge to the true model from the surface downward and the gradient method discriminates the shallow region more clearly than


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
 Gauss–Newton seismic waveform inversion 1379

Distance (km) Distance (km) Distance (km) Distance (km) Distance (km)
0.0 1.8 0.0 1.8 0.0 1.8 0.0 1.8 0.0 1.8
0.0
True model Initial model Initial Gradient Gauss-Newton
seismogram seismogram residuals residuals residuals

0.2

0.4
Time (s)

0.6

0.8

(a) (b) (c) (d) (e)

Figure 4. Comparison of vertical component seismograms and residuals from a shot-gather of five shot-gathers: (a) seismogram from the true model, (b)
seismogram from an initial model, (c) initial residuals, (d) final residuals from the gradient method (100 iterations) and (e) final residuals from the Gauss–Newton
method (10 iterations). For clarity, direct arrivals of (a) and (b) are removed. The scaling factor of the final residuals, (d) and (e), is 2 relative to the others.

Figure 5. Model updates from the first iterations by the gradient (a), and the Gauss–Newton method (b and c) corresponding to the gradient direction
(Jt d) and the filtered gradients ([Jt J + λIt I]−1 Jt d and [Jt J + λ 1 Pt P + λ 2 It I]−1 Jt d), respectively. All images are scaled by its own maximum
value.

the deeper region whereas that of the Gauss–Newton method dis- rate and the resolving ability, the Gauss–Newton method is more
criminates the whole region. It is clearly shown in Fig. 3(c) that, attractive than the gradient method or the backward propagation
even after the first iteration, the Gauss–Newton method reveals the method.
velocity structure at deep region of the model. Also the former
mostly accounts for the high spatial frequency contents of the up- 4.2 Example 2
date information but the latter contains the intermediate frequency
contents where is partly because of the Laplacian regularization op- The second example has complex geological structures consisting
erator. The damping regularization operator stabilizes the inversion of a fold, discontinuity and dipping layers (Fig. 6a). The model has
but artefacts from shortage of source–receiver coverage still remain variations on density and S-wave velocity which were assumed to
in the update (Fig. 5b). This is reduced by the Laplacian operator be constant in the previous example. The relationships (Castagna
(Fig. 5c). et al. 1993) between density, S-wave velocity and P-wave velocity
At each iteration, in this example, 5 simulations for the for- are assumed to be known as
ward simulation, 21 simulations for the reciprocal simulation with β (km s−1 ) = −0.055 α 2 + 1.017 α − 1.030,
a source at the receiver location and 5 simulations for finding an op-
ρ (g cm−3 ) = 1.5 α 0.225 . (23)
timal step length were required. The backward propagation method
would only need 15 simulations at each iteration because it requires After updating P-wave velocities, the other two elastic parameters
the solution of three times as many forward simulations as the shots are updated using these relationships.
(i.e. one to generate synthetic data, another to make the gradient, and The configurations, such as source properties, grid spacing and
the third to compute the step length). Considering the convergence time step, for the simulation and the inversion are same as before,


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
1380 D.-H. Sheen et al.

Figure 6. Inversion results of example 2: (a) The true P-wave velocity model and (b) inversion result after 5 iterations. Note that the images are at the same
scale. Remnant of the initial velocity model can be seen from outsides of updated area of (b) which has vertical velocity variation monotonically increasing
with depth.

-0.4 km 0.0 km 0.4 km

0.0 0.0 0.0

(a) (b) (c)

0.4 0.4 0.4
Depth (km)

0.8 0.8 0.8

1.2 1.2 1.2

True model True model True model

Initial model Initial model Initial model
1.6 Final model 1.6 Final model 1.6 Final model

3.2 3.6 4.0 4.4 4.8 5.2 3.2 3.6 4.0 4.4 4.8 5.2 3.2 3.6 4.0 4.4 4.8 5.2
P velocity (km s−1) P velocity (km s−1) P velocity (km s−1)
Figure 7. P-wave velocity profiles of example 2 at three locations (see Fig. 6). The locations of profiles are at −0.4 km (a), 0.0 km (b) and 0.4 km (c) of the
model. The final models are obtained after 5 iterations.

including the initial model. However, the inversion requires more In this example, 4200 blocks of model parameters are consid-
information on the model than that on the previous model because ered. To solve this problem, 9 forward simulations, 98 reciprocal
the level of complexity is significantly increased. Therefore, we simulations and 9 simulations for the step length are required at
increase the number of observed data. We use 9 shot gathers and a each iteration which takes about 14 min on the IBM cluster, de-
shot gather has 46 or 47 traces with a spacing of 40 m. The minimum scribed previously, with 48 CPUs. Forward and reciprocal simula-
and maximum offsets are 80 and 1920 m, respectively. Both of the tions take 3 and 27 per cent of total computing time, respectively, and
horizontal and vertical components of the observed data are used constructing the Jacobian and approximate Hessian matrix spends
for the inversion. 66 per cent. Note that the reciprocal simulation requires about
Fig. 6(b) shows that the structures are correctly resolved within 70 per cent of computation of the forward simulation. The prin-
the region of interest. The velocity profiles of the model are plotted cipal memory requirements are about 1.0 GB for saving the for-
in Fig. 7 which compares these to the true model and the initial ward wavefields, Fpm in eq. (18), about 75 MB for the Jacobian
model, showing that the velocities also agree with the true model. matrix, J p in eq. (21), and about 67 MB for the approximate Hes-
Thus, it may be concluded that the inversion with the Gauss–Newton sian matrix, Jt J in eq. (21), respectively. By virtue of the approxi-
method can identify a 2-D velocity structure very well and also mation of virtual sources, the largest memory requirement, 1.0 GB
the properties of the model can be accurately estimated after a few for the forward wavefields, is reduced to 78 MB without any loss of
iterations. In this example, the inversion using both of the horizontal accuracy.
and vertical components does not show a notable contribution to In order to invert the noise-corrupted data, we added 5 per cent
the resolving power. Although we do not show the result from the of Gaussian noise to the synthetic data. Fig. 8 shows seismograms
inversion considering only the vertical component of the observed without and with Gaussian noises. In example 1, density and S-wave
data, the velocities are also recovered accurately. velocity are assumed to be constant so that only P-wave arrivals are


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
 Gauss–Newton seismic waveform inversion 1381

As shown by Tarantola (1986), the choice of model parameters is

Distance (km) Distance (km) not neutral in the inversion. From numerical experiments, variation
0.0 1.9 0.0 1.9 in density doesn’t affect the inversion much. In order to compare
0.0
convergent rates of P- and S-wave velocities, these velocities are
updated independently of each other, but density is fixed during the
inversion. Noise-corrupted data is used as observed data. Fig. 10
shows P and S-wave velocity profiles of the model. P-wave velocity
0.2 of the model is recovered well up to the last reflector, whereas the
inversion couldn’t recover S-wave velocity profile at depths larger
than 1 km as the data has limited information of S-wave velocity
beyond this depth.

0.4
Time (s)

4.3 Example 3
In this example, we consider a more realistic case instead of apply-
ing this method to real data. The data consists of 20 shot gathers
0.6 collected with a total of 125 receivers and has 50 per cent of Gaus-
sian noise (Fig. 11). The shot and receiver spacings are 600 and
120 m, respectively. The maximum offset of each shot gather is
6 km. The time step for wave propagation simulation is 0.001 s,
total record length is 2.25 s, and a dominant frequency of 8 Hz of
0.8 source is used. The length and depth of region of interest is 14 and
3.5 km, which consists of 9800 blocks of model parameters. The re-
lationships between density, S- and P-wave velocities are assumed
(a) Noise-free data (b) Noise-corrupted to follow eq. (23).
data The inversion results are obtained after 3 iterations and shown
in Fig. 12. The result from noise-free data (Fig. 12b) shows ex-
Figure 8. Vertical component seismograms of a shot-gather from example 2: cellent recovery of the true model and that from noise-corrupted
(a) original seismogram and (b) noise-corrupted seismogram. Direct arrivals data (Fig. 12c) also shows good recovery. Higher contamination
are removed. in deeper region can be explained by relatively low signal to
noise ratio of late arrivals. Fig. 13 shows the velocity profiles
of the model. The agreement with the true velocity profile is
observed. In this example, we consider all parameters have varia- quite satisfactory, considering the crudeness of the noise-corrupted
tions and, thus, it shows diverse P and S-wave arrivals. The inversion data.
of the noise-corrupted data was performed in the same way as be- The inversion is performed with 48 CPUs and each iteration takes
fore. The result of the inversion, Fig. 9, shows that the noise slightly 3 hr mainly because of construction of the Jacobian and approximate
contaminates the model, but still yields an accurate identification of Hessian matrices, taking 83 per cent of computing time. The for-
the model. ward and reciprocal simulations take 3 and 12 per cent, respectively.

-0.4 km 0.0 km 0.4 km

0 0 0
(a) (b) (c)
0.4 0.4 0.4
Depth (km)

0.8 0.8 0.8

1.2 1.2 1.2

True model True model True model

Noisy data Noisy data Noisy data
1.6 Clean data 1.6 Clean data 1.6 Clean data

3.2 3.6 4 4.4 4.8 5.2 3.2 3.6 4 4.4 4.8 5.2 3.2 3.6 4 4.4 4.8 5.2
P velocity (km s−1) P velocity (km s−1) P velocity (km s−1)

Figure 9. P-wave velocity profiles using noise-corrupted data of example 2 at three locations. The locations of profiles are at −0.4 km (a), 0.0 km (b) and
0.4 km (c) of the model. Density and S-wave velocity are updated by using eq. (23).


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
1382 D.-H. Sheen et al.

-0.4 km 0.0 km 0.4 km

0 0 0
(a) (b) (c)
0.4 0.4 0.4
Depth (km)

0.8 0.8 0.8

1.2 1.2 1.2

True model True model True model

Initial model Initial model Initial model
1.6 Final model 1.6 Final model 1.6 Final model

3.2 3.6 4.0 4.4 4.8 5.2 3.2 3.6 4.0 4.4 4.8 5.2 3.2 3.6 4.0 4.4 4.8 5.2
P velocity (km s−1) P velocity (km s−1) P velocity (km s−1)

0 0 0
(d) (e) (f)
0.4 0.4 0.4
Depth (km)

0.8 0.8 0.8

1.2 1.2 1.2

True model True model True model

Initial model Initial model Initial model
1.6 Final model 1.6 Final model 1.6 Final model

1.8 2.0 2.2 2.4 2.6 2.8 1.8 2.0 2.2 2.4 2.6 2.8 1.8 2.0 2.2 2.4 2.6 2.8
S velocity (km s−1) S velocity (km s−1) S velocity (km s−1)

Figure 10. P- and S-wave velocity profiles using noise-corrupted data of example 2: (a)–(c) P-wave profiles, and (d)–(f) S-wave profiles. S-wave velocity is
updated independently of P-wave velocity but density is forced to be constant. The locations of profiles are at −0.4 km (a and d), 0.0 km (b and e), and 0.4 km
(c and f) of the model.

Comparison with the previous example indicates that the percent- the approximate Hessian matrices have been carried out by the
age increase of the computation time for the matrices is due to the reciprocity principle and the convolution theorem, which signifi-
increase of number of model parameter. cantly reduces the number of forward simulations. To overcome the
computational and memory limitation, we used different sizes of
grids in the spatial and time domains for the inversion, temporal
windowing for the partial derivative and approximation of virtual
5 C O N C LU S I O N S
sources, and parallelized all these approaches on massively parallel
Huge computational and memory requirements have prevented computers.
solving the Gauss–Newton seismic waveform inversion (Tarantola In numerical experiments, the Gauss–Newton seismic waveform
1984). Unfortunately, it is still a difficult problem. In this pa- inversion has been successfully applied to resolve the high and in-
per, we have implemented a seismic waveform inversion scheme termediate spatial frequency content of the model. In particular,
based on the Gauss–Newton method. We used the velocity-stress the results show that the Gauss–Newton method has faster conver-
staggered-grid finite difference method (Levander 1988) for sim- gence rate than the gradient method for seismic waveform inver-
ulating seismic wave propagation in elastic media and the PML sion. Furthermore, not only does the method recover the structure
method (Berenger 1994) as an absorbing boundary condition. This correctly, but also estimates velocities accurately. Further work
approach allows us to accurately simulate seismic wavefields, use is required to study the effect of the choices of model param-
multicomponent data, and resolve subsurface elastic features cor- eters and, most importantly, to test the method for real seismic
rectly and efficiently. Explicit computations of the Jacobian and data.


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
 Gauss–Newton seismic waveform inversion 1383

-1 km 1 km
0
(a) (b)
1

Depth (km)
2

4 True model True model

Initial model Initial model
Final model Final model
5
3.2 3.6 4 4.4 4.8 5.2 3.2 3.6 4 4.4 4.8 5.2
P velocity (km s−1) P velocity (km s−1)

Figure 13. P-wave velocity profiles of example 3: The locations of profiles

are at −1 km (a) and 1 km (b).

REFERENCES
Aki, K. & Richards, P.G., 1980. Quantitative Seismology: Theory and Meth-
ods, Vol. 1, W.H. Freeman and Company, San Francisco, California.
Arntsen, B. & Carcione, J.M., 2000. A new insight into the reciprocity prin-
ciple, Geophysics, 65, 1604–1612.
Berenger, J.P., 1994. A perfectly matched layer for the absorption of elec-
Figure 11. Vertical component seismograms for example 3: (a) original seis- tromagnetic waves, J. Comput. Phys., 114, 185–200.
mogram and (b) noise-corrupted seismogram. Direct arrivals are removed. Castagna, J.P., Batzle, M.L. & Kan, T.K., 1993. Rock physics—the link be-
tween rock properties and avo response, in Offset-dependent reflectivity–
Theory and Practice of AVO Anomalies, pp. 135–171, eds Castagna, J.P. &
Backus, M.M. no. 8 in Investigations in Geophysics, Soc. Expl. Geophys.,
Tulsa.
Clément, F., Chavent, G. & Gómez, S., 2001. Migration-based traveltime
waveform inversion of 2-d simple structures: A synthetic example, Geo-
physics, 66, 845–860.
Eisner, L. & Clayton, R.W., 2001. A reciprocity method for multiple-source
simulations, Bull. seism. Soc. Am., 91, 553–560.
Freudenreich, Y. & Shipp, R., 2000. Full waveform inversion of seismic
data: frequency versus time domain, Lithos Science Report, 2, 25–30.
Gauthier, O., Virieux, J. & Tarantola, A., 1986. Two-dimensional nonlinear
inversion of seismic waveforms: numerical results, Geophysics, 51, 1387–
1403.
Hicks, G.J. & Pratt, R.G., 2001. Reflection waveform inversion using lo-
cal descent methods: estimating attenuation and velocity over a gas-sand
Figure 12. Inversion results of example 3: (a) The true P-wave velocity deposit, Geophysics, 66, 598–612.
model, (b) inversion result from noise-free data, and (c) from noise-corrupted Kolb, P., Collino, F. & Lailly, P., 1986. Pre-stack inversion of a 1-d medium,
data. Both inversion results are obtained after three iterations. Note that the Proceedings of the IEEE, 74, 498–508.
images are at the same scale. Lailly, P., 1983. The seismic inverse problem as a sequence of before stack
migrations, in Conference on Inverse Scattering: Theory and Application,
pp. 206–220, eds Bednar, J.B., Redner, R., Robinson, E. & Weglein, A.,
AC K N OW L E D G M E N T S Soc. Ind. Apl. Math., Philadelphia.
Levander, A.R., 1988. Fourth-order finite-difference p-sv seismograms,
This work was partially supported by a grant from the Office of Sci- Geophysics, 53, 1425–1436.
ence (# DE-F602-91ER14175), a contract from the Office of Fossil Levenberg, K., 1944. A method for the solutions of certain nonlinear prob-
Energy (# DE-RA26-99FT40160) of the United States Department lems in least squares, Quart. Appl. Math., 2, 164–168.
of Energy, the BK21 program through the School of Earth and En- Marquardt, D.W., 1963. An algorithm for least-squares estimation of non-
vironmental Sciences, Seoul National University, and a project of linear parameters, J. Soc. Indust. Appl. Math, 11, 431–441.
Mora, P., 1987. Nonlinear two-dimensional elastic inversion of multioffset
development of fundamental techniques for evaluating earthquake
seismic data, Geophysics, 52, 1211–1228.
hazard from the Korea Meteorological Administration. We would Pica, A., Diet, J.P. & Tarantola, A., 1990. Nonlinear inversion of seismic
like to acknowledge Mary Papakhian (University Information Tech- reflection data in a laterally invariant medium, Geophysics, 55, 284–292.
nology Services, Indiana University) for reserving CPUs on the IBM Pratt, R.G., 1990. Inverse theory applied to multi-source cross-hole tomog-
SP cluster. This work was also supported in part by facilities pro- raphy. part 2: elastic wave-equation method, Geophys. Prospect., 38, 311–
vided through Shared University Research grants from IBM, Inc. 329.


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS
1384 D.-H. Sheen et al.

Pratt, R.G., Shin, C. & Hicks, G.J., 1998. Gauss-Newton and full newton where eq. (A1) is equvalent to eq. (17). For model parameters of
methods in frequency-space seismic waveform inversion, Geophys. J. Int., density, ρ, and P- and S-wave velocities, α and β, the equations are
133, 341–362. given by
Rodi, W.L., 1976. A technique for improving the accuracy of finite element
solutions for magnetotelluric data, Geophys. J. R. astr. Soc., 44, 483– ∂ v̇i ∂τi j, j ∂ρ
ρ = − v̇i ,
506. ∂ρ p ∂ρ p ∂ρ p
Sasaki, Y., 1989. Two-dimensional joint inversion of magnetotelluric and  
dipole-dipole resistivity data, Geophysics, 54, 254–262. ∂ τ̇i j ∂vi, j ∂v j,i ∂vk,k 1 ∂ρ (A4)
=μ + + λ δi j + τ̇i j ,
Sheen, D.-H., Tuncay, K., Baag, C.-E. & Ortoleva, P.J., 2004. Efficient fi- ∂ρ p ∂ρ p ∂ρ p ∂ρ p ρ ∂ρ p
nite difference calculation of partial derivative seismic wavefield using
reciprocity and convolution, in 74th Annual meeting, SEG. ∂ v̇i ∂τi j, j
Shin, C., Yoon, K., Marfurt, K.J., Park, K., Yang, D., Lim, H.Y., Chung, S. & ρ = ,
Shin, S., 2001. Efficient calculation of a partial-derivative wavefield using
∂α p ∂α p
reciprocity for seismic imaging and inversion, Geophysics, 66, 1856–  
∂ τ̇i j ∂vi, j ∂v j,i ∂vk,k
1863. =μ + + λ δi j
Shipp, R.M. & Singh, S.C., 2002. Two-dimensional full wavefield inversion ∂α p ∂α p ∂α p ∂α p
of wide-aperture marine seismic streamer data, Geophys. J. Int., 151, 325–
∂α (A5)
344. +2 ρ α δi j vk,k ,
Sirgue, L. & Pratt, R.G., 2004. Efficient waveform inversion and imaging: a ∂α p
strategy for selecting temporal frequencies, Geophysics, 69, 231–248. and
Snieder, R., Xie, M.Y., Pica, A. & Tarantola, A., 1989. Retrieving both the
∂ v̇i ∂τi j, j
impedance contrast and background velocity: a global strategy for the ρ = ,
seismic reflection problem, Geophysics, 54, 991–1000. ∂β p ∂β p
Symes, W.W. & Carazzone, J.J., 1991. Velocity inversion by differential
∂ τ̇x x ∂vx,x ∂vz,z ∂β
sembrance optimization, Geophysics, 56, 654–663. = (λ + 2μ) +λ − 4ρ β vz,z ,
Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic ap- ∂β p ∂β p ∂β p ∂β p
proximation, Geophysics, 49, 1259–1266.
Tarantola, A., 1986. A strategy for nonlinear elastic inversion of seismic ∂ τ̇zz ∂vx,x ∂vz,z ∂β
=λ + (λ +2μ) − 4ρ β vx,x ,
reflection data, Geophysics, 51, 1893–1903. ∂β p ∂β p ∂β p ∂β p
Tarantola, A., 1987. Inverse Problem Theory: Methods for Data Fitting and  
Parameter Estimations, Elsevier Science Publ. Co., New York. ∂ τ̇x z ∂vx,z ∂vz,x 2 ∂β (A6)
=μ + + τ̇x z .
Tarantola, A. & Valette, B., 1982. Generalized non linear inverse prob- ∂β p ∂β p ∂β p β ∂β p
lems solved using the least-squares criterion, Rev. of Geophys. and Space
Physics, 20, 219–232. For model parameters of density, ρ, and P- and S-wave impedances,
IP and IS, the equations are given by
∂ v̇i ∂τi j, j ∂ρ
A P P E N D I X A : PA RT I A L D E R I VAT I V E S ρ = − v̇i ,
∂ρ p ∂ρ p ∂ρ p
In this Appendix we present the differentiations of the wave equa-  
∂ τ̇i j ∂vi, j ∂v j,i ∂vk,k 1 ∂ρ (A7)
tion with respect to model parameters for different choices of pa- =μ + + λ δi j − τ̇i j ,
rameters. In the case of model parameters of density, ρ, and one of ∂ρ p ∂ρ p ∂ρ p ∂ρ p ρ ∂ρ p
Lamé’s moduli, λ and μ, the differential equations are given by
∂ v̇i ∂τi j, j
∂ v̇i ∂τi j, j ∂ρ ρ = ,
ρ = − v̇i , ∂ I Pp ∂ IPp
∂ρ p ∂ρ p ∂ρ p  
  ∂ τ̇i j ∂vi, j ∂v j,i ∂vk,k
∂ τ̇i j ∂vi, j ∂v j,i ∂vk,k =μ + + λδi j
=μ + + λ δi j , (A1) ∂ IPp ∂ IPp ∂ IPp ∂ IPp
∂ρ p ∂ρ p ∂ρ p ∂ρ p
IP ∂ IP (A8)
∂ v̇i ∂τi j, j +2 δi j vk,k ,
ρ = , ρ ∂ IPp
∂λ p ∂λ p
  and
∂ τ˙i j ∂vi, j ∂v j,i ∂vk,k ∂λ
=μ + + λ δi j + δi j vk,k , (A2) ∂ v̇i ∂τi j, j
∂λ p ∂λ p ∂λ p ∂λ p ∂λ p ρ = ,
∂ IS p ∂ IS p
and
∂ τ̇x x ∂vx,x ∂vz,z IS ∂ IS
∂ v̇i ∂τi j, j = (λ +2μ) +λ −4 vz,z ,
ρ = , ∂ IS p ∂ IS p ∂ IS p ρ ∂ IS p
∂μ p ∂μ p
  ∂ τ̇zz ∂vx,x ∂vz,z IS ∂ IS
∂ τ˙i j ∂vi, j ∂v j,i ∂vk,k =λ + (λ +2μ) −4 vx,x ,
=μ + + λ δi j ∂ IS p ∂ IS p ∂ IS p ρ ∂ IS p
∂μ p ∂μ p ∂μ p ∂μ p
 
∂μ ∂ τ̇x z ∂vx,z ∂vz,x 2 ∂ IS (A9)
+ (vi, j + v j,i ). (A3) =μ + + τ̇x z .
∂μ p ∂ IS p ∂ IS p ∂ IS p IS ∂ IS p


C 2006 The Authors, GJI, 167, 1373–1384
Journal compilation 
C 2006 RAS