GEOPHYSICS, VOL. 53, NO. 1 (JANUARY 1988); P. 2131, 11 FIGS.
Fullwaveform inversion of planewave seismograms in stratified acoustic media: Theory and feasibility
G. S. Pan*, R. A. Phinney*, and R. I. odomr
prescribed set of mathematical transformations to the observed data, based on an exact or approximate theory derived from conventional forward scattering theory, to model the data. For a layered medium, this requires solving some form of the GelfandLevitanMarchenko integral equation (Ware and Aki, 1969; Nilsen and Gjevik, 1978; Berryman and Greene, 1978; Berryman, 1979). Although some quite interesting and powerful algorithms of this sort have been described, they do not readily lend themselves to considering the effects of noise on the inversion. Inverse scattering work in two or three dimensions generally leads to some sort of generalization of the familiar process of wavefield migration, such as the BornWKBJ method of Clayton and Stolt (1981) or Tarantola's (1984a) formulation using linearized least squares. In inversion studies, no attempt is made to find a mathematical operator which converts the data into a model; rather, a known algorithm for modeling the data from an earth model is employed in an iterative procedure. Emphasis is placed on developing more information in the process than merely a single model which fits the data. This process includes an assessment of the uniqueness problem (i.e., what is the class of models which fit the data in some sense?) and the limits on information imposed by the presence of noise. In seismology, the most familiar of such studies have taken the leastsquares approach; the inverse theory developed for the normal modes of the earth is now a standard tool in developing models for the crust, mantle, and core, and for assessing the variance and resolving power of the results (Backus and Gilbert, 1967, 1968, 1970; Gilbert, 1971; Franklin, 1970; Jackson, 1972, 1973, 1979; Jordan and Franklin, 1971; Wiggins, 1972; Wiggins et al., 1973). In a recent series of papers, Tarantola and Valette (1982a, 1982b) and Tarantola (1984b) proposed a nonlinear leastsquares inversion method which lends itself to application in reflection seismology. In this paper we adopt the Tarantola formalism and specialize it to the problem of a medium whose density and elastic properties depend strictly on depth alone, where data collected at several sourcereceiver offsets contain the infor
ABSTRACT
We have implemented an inversion procedure for obtaining velocity and density profiles from multioffset data in a layered acoustic medium. Using an iterative modeling technique and the pr; representation, the procedure is derived from the nonlinear leastsquares formalism of Tarantola and Valette. The feasibility of this method depends upon obtaining Frechet derivatives during the modeling process and on vectorizing the ThompsonHaskell reflectivity algorithm. Test data sets for this study are gathers of two to four planewave synthetic seismograms which may contain both precritical and postcritical arrivals generated by a seismic wavelet. These examples demonstrate convergence to the correct velocity and density profiles with reasonable accuracy.
INTRODUCTION
A common objective in reflection seismology is to determine the velocity and density profiles of a horizontally layered subsurface from seismic reflection data. There are three approaches to this problem: conventional data processing, an inverse scattering method, or an inversionbased approach. Conventional processing subjects the observed data to a sequence of transformations (frequently linear) which are designed to convert observed traces to some kind of image of the section, out which do not map wavefield dynamics correctly into elastic moduli and density. Such processing has proven to be a methodology of great practical power and has broad application to twodimensional (2D) and threedimensional (3D) problems; however, major difficulties hamper the process of extracting velocity and density from the betterdetermined impedance function. An inverse scattering approach involves applying a suitably
Manuscript received by the Editor August 25, 1986; revised manuscript received June 11, 1987. *Department of Geological and Geophysical Sciences, Guyot Hall, Princeton University, Princeton, NJ 08540. tFormerly Department of Geological and Geophysical Sciences, Princeton University; presently Honeywell Inc, Marine Systems Division, 6500 Harbour Heights Parkway, Everett, W A 98204. ? 1988 Society of Exploration Geophysicists. Allrightsreserved.
21
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
22
Pan et al.
mation needed to extract these unknown depth functions. We have simplified the problem somewhat in the interest of having a fairly streamlined procedure during these early stages of testing and evaluation; the model is considered to be purely acoustic, without attenuation, and the data are considered to consist of a set of transient plane waves rather than a set of different offset geophone recordings. The description of the data in terms of plane waves is particularly suitable for the purely layered media under consideration. The procedure is formulated in terms of horizontal slowness p and intercept time 't in which each trace is the seismogram for a plane wave with a particular angle of incidence at the surface. Previously, Clayton and McMechan (1981), Phinney et al, (1981), Stoffa et al. (1981), and Treitel et al. (1982) showed that seismic reflection arrivals of the distancetime (rt) section for a stratified medium are mapped into a series of quarterellipse trajectories in the p? section, and that each trace in the p't section corresponds to the transient response of the system to an incident plane wave specified by the particular p value. By setting up the inverse problem in terms of planewave data, we seek to infer the velocity and density functions from two or more planewave reflection seismograms with different angles of incidence. We are thereby restricted to data which are adequately sampled with offset to permit the transformation of the data (a source gather with a large number of channels) into planewave seismograms. Carrion et al. (1984) performed an experiment in the p: domain using reflection coefficients in a single scattering approximation with synthetic data. They were able to recover the onedimensional (ID) density and velocity profiles correctly from only two planewave traces generated by an infinite broadband source; they were less successful with a bandlimited source. In this paper, we rederive Tarantola's theory for the specific case of horizontally stratified media and show the form which the Frechet derivatives take for plane waves. Specific examples are chosen to demonstrate the feasibility of the methodology, using lownoise synthetic data with a minimumphase bandlimited (854 Hz) source wavelet. The examples consist of two to four trace p't gathers where the noise is restricted to the numerical error of the modeling process. For the demonstration models reported here, and for more realistic models described in a companion paper (Pan and Phinney, 1987), the inversion procedure proves to be capable of extracting the properties of the model with reasonable accuracy, and it can extract them in the presence of noise.
(I)
where
Vi={~[:r(r :r)J+~a~~}
is the horizontal Laplacian operator. The Radon transform (Deans, 1983), which maps conventional offsettime (rt) solutions of the wave equation to the (p't) domain, must have a form which takes the cylindrical coordinate system into account. It can be obtained as a generalized form of the slant stack (Chapman, 1981) or by a detour into the frequency domain and an application of the FourierHankel transform (Treitel et ai., 1982). The latter approach enables one to obtain expressions for the Radon transform of the wave equation and its solutions. Applying the FourierHankel transform of order 0 to equa
·· .. · ·· .. ·· .. . ·· .. ·· .. : . : : : : ·8: : : : : G: : : : : : : : : : : . . ::::: :.:::: :.:.:.:.:.:.:.:.:.:.:.:.:.: · .. ·· .. · ·· .. ·· .. · .
? ???????? '.0 ?????????????
· · . . . .. · · . . . .. · . . . .. · · . . . . . .. .. · · . . . .. ·
. . . .. .
.. .. .. .. ..
THE ACOUSTIC WAVE EQUATION
Consider a layered acoustic medium embedded between two halfspaces (Figure 1). The material properties, bulk modulus K(Z) [or velocity c(z)], and density p(z) are taken to depend upon depth z only. The layering is regarded as a device for approximating arbitrary behavior of realistic models. Given cylindrical coordinates (r, <1>, z) with a point source at (r, z, t) = (0, 0, 0), the acoustic pressure field P and the source function s are related by
:::::::::::~::::::::::::::::::::::::::::::::::: ·· .. ·· .. ·· .. ·· .. ·· ..
FIG. 1. Geometry of the model. The model is a layered acoustic medium embedded between two halfspaces. The material properties, velocity C j and density Pi' are constant within each microlayer. The source and receivers are at the same level.
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
Multioftset Fullwaveform Inversion
23
tion (1),we get
(  ioof (kr )2 0 {   +p(Z)    K(Z) oz p(Z) OZ
= 8(z)s(00),
[I
o]}_ P(k
r"
Z
(0)
(2)
It is well known that a normalincidence seismogram can be mapped (under noisefree Circumstances) into a plot of changes in acoustic impedance (pc) with depth, and vice versa. In exactly the same fashion, a planewave seismogram at some nonzero value of p can be mapped into the depth variations of the modified impedance function
where s(oo) is the source spectrum, and k, is the horizontal wavenumber along the raxis. Note that we use a tilde and a circumflex to denote quantities in the (p, (0) domain and the (p, 1:) domain, respectively. The use of order 0 implies the wave field is produced by a verticalforce vibrator or by an explosion. We now evaluate equation (2) along constant p = kr/oo to get
( ioo)2 { K(Z)
=
p/q = p ( C Z  pZ
I
)I/Z
.
(6)
+ (oop)2 _ ~
p(Z)
oz
[_1 ~]} oz
p(Z)
P(oop
' ,
Z
(0)
8(z)s(00).
(3)
After some manipulation, equation (3) becomes a Helmholtz equation with source p(z)8(z)s(00):
 {P(Z)
~ [_1_ ~J + oo2 q2}p az p(Z) oz
p(ooP,
z, (0)
This dependence on incidence angle suggests that two different planewave seismograms can be inverted to two such modified impedance functions, from which both velocity and density can be calculated. In a similar way, in an elastic medium, three planewave seismograms should suffice to yield density, compressional velocity, and shear velocity. In practice, of course, the presence of noise and the illconditioning of the inversion for some depths indicate that more than this minimum number of data traces would be needed (Phinney et aI., 1981). A planewave section of n traces consists of n distinct experiments which sample the same structure. By the same token, n seismic traces from individual receivers in an array provide redundant experiments on a single structure.
THE INVERSE PROBLEM
= p(z)8(z)s(00).
(4)
The subscript p indicates that the horizontal slowness p is a constant parameter. We have obtained this result using the relations
In this section we review Tarantola's formulation of the inversion process to establish the notation and basic characteristics of what follows. Consider a class of models m and the class of data d, written as vectors
00 k=, c q
1
=
(
C
Z

)1/2 p2 ,
and
(7)
and
d=
where k, k z , a, and q are the wavenumber, the vertical wavenumber, the slowness, and the vertical slowness, respectively. After an inverse temporal Fourier transform, we obtain the ID acoustic wave equation in (z, 1:) for a fixed p; namely,
0 0  { p(z) ~ az p(z) az = p(z)8(z)S(1:).
[
Q>1(t1)' Q>1(tZ)' ... , Q>2 (t l), Q>z (t2)' ... , Q>.(t l ), Q>.(tz), ...
T'
(8)
[1
J [2(z)  p zJ 0 z a
2
where Q>. is the nth seismic trace and T denotes the transpose. Then
d =fCm)
}
01:
P (z, 1:)
A
p
(5)
Signals propagate in the z direction at the speed (a Z _ p2)I/Z. The p value, or ray parameter, uniquely specifies a particular raypath in a stratified medium. Since the input to our inversion procedure is a P1: section, each individual trace of the section can be described by the ID equation (5), for a particular fixed value of p. From this equation, we obtain the Green's function and Frechet derivatives needed to invert planewave data (next section). We also need the transient planewave response of a stack of layers, and base this on the elementary form of the left side of equation (5), using any of a number of familiar methods.
is the statement of the forward problem, where f is the wave equation operator. Given an observed data set do and a starting model mo chosen on the basis of ancillary information, in general do # f(m o)' We seek a model m minimizing the quadratic form
(d  do)T{:d;, I(d  do) + (m  rnof{:';;ol(rn  mo), (9)
where {:d is the covariance matrix for the data and {:m is the covariance matrix for the model. In function analytic terms, we want to find a point (d, m) that defines the shortest distance to the given point (do, mo) in the I3 sense in the combined (data, model) space, using the weighted norm equation (9). If one views the inversion prob
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
24
Pan et al.
abilistically, the goal of the inversion is to maximize the product of probability density functions X(d) and x(m) of the multivariate Gaussian distribution with covariance matrices (do and (mo' where X(d) and x(m) are X(d) = ;~i/2 x(m) =
IC'I'/2 exp [1 (d  dO)T(d;, '(d  do)] , "2
1('1'/2
i12
2:
exp
[1 "2
(m  mof(';;o'(m  mo) ,
]
show that we obtain not only the synthetic seismograms, but also the Green's function and the Frechet derivatives as byproducts of the computation. The analytical expressions of Frechet derivatives, which are essential to the computational feasibility of the method, were derived in a general form by Tarantola (1984b) and are rederived below in the form required for plane waves in a layered medium. The synthetic seismograms The well known ThompsonHaskell matrix method may be used to set up a solution for the synthetic planewave seismograms in the poo domain; the result is then converted to the pi domain by an inverse Fourier transform. The model is considered to consist of a large number of homogeneous microlayers. Let K i and Pi be the bulk modulus and density, respectively, of the ith layer (Figure 1). In the ith layer, the equation of propagation of plane waves [equation (5)] reduces to (14) In the frequency domain, equation (14) becomes
and i and j are the dimensions of the data and model spaces, respectively. Hence we minimize the exponential terms, that is, equation (9), to obtain the largest probability. This is the wellknown maximum a posteriori estimation. The inverse problem described by equations (7), (8), and (9) has been shown by Tarantola and Valette (1982a, b) to be solvable by the following iterative scheme:
mk + = mk + W {(mF[Ci' [do  !(m k ) ] ,
W=(l+(mF[(i'fk)',

(m,  m o)} , (11)
(12)
where f = a<\>/am denotes the Frechet derivative and the subscript k denotes the kth iteration. The use of the a priori model mo is to constrain the inverted model to be close to the starting model [in the sense of the norm equation (9)], thus reducing the chance of instability. The data covariance matrix (d describes the secondorder statistics of the data. Similarly, the model covariance matrix (m describes the secondorder statistics of the model; the diagonal values give a variance estimate, and the neardiagonal values characterize the resolution of the model if the model parameters are naturally ordered. In order to reduce the computational load, several simplified forms of the W operator were also suggested (Tarantola, 1984b),for example,
(::2 + ? )
002 q
P(oop, z, (0)
=
Pis(oo)8(z).
(15)
We can set up a suitable wave solution in the model by patching together homogeneous solutions to equation (15) in each layer so that the vertical displacement and acoustic pressure are continuous at each interface. The infinitesimal displacement u and the displacement potential 'I' are connected to P by
P= KY'U;
u=Y'I';
W = aI,
where a is a constant, or W = diag (l + (mf~(i'Fo)'; or
(13)
in the pt? domain, they are
P = KV 2 1jJ = poo 21jJ.
To generate the synthetic seismograms in the pt? domain, a series of Haskell matrices is set up to propagate the pressure field P and the vertical displacement ii z through a stack of layers in terms of upgoing and downgoing displacement potentials ('1'., 'I'd) by prescribing the displacement potentials to have only a downgoing part in the lower halfspace, i.e., ['1'., 'I'd]~ = [0, 1]~ (Figure 1).The equation is
alternatively, one can simply replace Fo by Fk to get two additional forms. From equation (11), it is clear that the synthetic seismograms and the Frechet derivatives must be computed at each iteration. This stresses the importance of a rapid forward algorithm.
THE FORWARD PROBLEM
[~dJ
where
~o
= EoTo'I, ...
~n2 I;\ In,En,I;'tTnl~dJ ' ~n
(16)
(17)
Inversion by iterative modeling is feasible only if the forward modeling step is fast and accurate. In the iterative scheme, using a gradientrelated method, the synthetic data and the gradient (Frechet derivative) of the current model are needed at each iteration in order to update the model for the next iteration. In this section we review the solution of the forward problem in terms of ThompsonHaskell matrices, and
and (18) In't In' the interface matrix, brings the displacement potentials across the nth interface by first converting the potentials
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
Multloffset Fullwaveform Inversion
25
into pressure and vertical displacement to match the boundary conditions at the interface and then converting them back into potentials again. En' the phase matrix, propagates the displacement potentials from the bottom to the top of the nth layer. The reflectivity function is obtained by taking the ratio of upgoing waves to downgoing waves at the receiver level. The synthetic seismograms are obtained by multiplication with the source spectrum, which is a simulated bandlimited (854 Hz) minimumphase source spectrum, followed by an inverse Fourier transform: (19)
strong pulse at the twoway delay corresponding to depth z. (An example of this phenomenon is illustrated in an example in the next section of this paper.) The actual expressions (20) and (21) are, of course, exact and include the effects of all internal multiples. Considering that the source(s) and receivers are at the same level in our model and that, in the pt? domain, convolution are replaced by multiplication and ( im), equations and (20) and (21) can be expressed in simpler forms in terms of the Green's function:
a/at
U(Zg, ro; Zs I z) =
l
  gp(z, m; zs)
( im) K(Z)
J2 s(m). Zs)J2
(22)
The Frechet derivatives In this study the source function is assumed to be known; consequently, the Frechet kernels are composed of only two parts, i.e., of/am = (OPp/OK, oPp/op). For a ID model [K = K(Z) and p = p(z)], the Frechet derivatives in the pot domain can be expressed by the following two equations using a standard perturbation analysis following Tarantola (1984b):
OPiZg' t; ZS' 0) OK(Z)
=
V(Zg, m; Zs I z) =
i ([!!...
p (z)
OZ
gp(z, ro;
(23)
The Green's function
The Green's function is defined as the response of the medium at some level Z to a unit impulse pointsource excitation at another level zo' In our application, one or the other of these responses is in the upper halfspace; owing to the reciprocity of Green's functions, we can always take Zo at the top of the layered stack. In the ThompsonHaskell reflectivity procedure, we build tables of the complex values of the upgoing and downgoing potentials within each layer in the course of calculating the synthetic seismograms. To calculate the Green's function, we simply renormalize each table. The resultant table then lists the values of the Green's function for each microlayer. The Green's function for our system is the solution of (24)
U(Zg, r ; Zs I z)
=
1
2K (z)
gp(Z, t ; Zg' 0)
(20)
oPp(Zg' t; Zs' 0) op(z)
=
V(Zg' r ; Zs I z)
= p:(Z) {[:Z gp(Z, t ; Zg' 0)
* :Z gp(Z, r ; Zs' 0) * S(t)J
Consider a small region near the source and transform equation (24) into the pt? domain:
 ( dz2
 [p2 gp(Z, r ; Zg' 0) * gp(Z, r ; Zs' 0) * ::2 S(t)J},
(21) where gp (z, r ; Zg' 0) denotes the pressure response at (z, t] due to a delta function source at (Zg, 0), sand g denote the source and the geophone, respectively, and the asterisk denotes time convolution. The basic physics of these expressions may be understood with reference to equation (20). By reciprocity, gp (z, t ; Zg' 0) = gp(Zg' 0; Z, t); therefore, the Frechet kernel of the recorded data at the geophone with respect to an infinitesimal change in K at depth Z is obtained from the convolution of the impulse response at Z to a source at the shotpoint with the impulse response at the geophone to a source at z. Since both impulse responses (that is, Green's functions) are dominated by a simple pulse with a delay, the Frechet derivative shows a
d2
+ m22) g(mp, z, ro) = qo
Poo(z).
(25)
The solution of equation (25) which satisfies the radiation condition at ± 00 is ipo . g(mp, z, ro) =   e,mqo I%I? 2mqo (26)
Throughout the stack of layers from the source level downward, the Green's function may be represented by a homogeneous solution of equation (24) such that the downward traveling wave just below the source has the amplitude ipo/2mqo' Consequently, we obtain the Green's function by rescaling the homogeneous solution produced in the course of running the ThompsonHaskell matrix product equation (16) through the layers while generating the synthetic seismogram. An inverse Fourier transform then yields the Green's function in the time (t) domain. Shown in Figure 3a is the Green's function as seen
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
26
Pan et al.
through a finite band filter, for a 36microlayer representation of the model in Figure 2, for two values of p. Microlayer thickness is chosen to be 1/4 of the shortest wavelength. Since internal multiples are weak, the Green's function is dominated by the direct signal between the source and depth, and may be regarded as a banddiagonal matrix, with suitable labeling of the coordinates. By the same token, the Frechet derivatives (Figures 3b) are banddiagonal operators, a property which can be exploited to speed convergence of the inversion. It should be clear that the practical success of the inversion will depend upon a choice of microlayer thickness which is neither too coarse to model the data well nor so fine that computational problems become important.
interfaces. Since the reflectivity function is a ratio, any rescaling of the wave functions does not affect the result. In equation (16), we revise In to be
In=(iO)[
qn
n Pn] P qn
and make the interface matrices frequencyindependent. The
p=o
20 36
p
=.636
20
36
MICIfOLAYER
NUIf8ER
Numerical considerations
To implement equation (16), it is advisable to separate the 0) terms out of the In matrix in order to use anyone of the familiar methods, such as the spectral method, the reflectivity method (Fuchs and Muller, 1971), or the slowness method (Chapman, 1978a, b). To separate the 0) terms out, instead of P and iiz, the quantities PliO) and iiz are propagated through the stack of layers. Since in the frequency domain each pass of the calculation is done at a fixed frequency, the denominator iO) does not enter into the boundary conditions at the internal
(a)
p=o
20
p
=.636
20
36
MICAOLAYER HUMBER
36
0
55
I
I
I I
~
E
I I
~110
Ii: i!l
165
220 0.75
DENSITY [gIcm'}
~
II
I
I
1.25
(al
(b)
1.00
1.50
1.75 xl0'
VELOCITY [m/s]
SLOWNESS P
.000 .125 250
375
500
[slkm]
.~~
z ....
."
.~.
ill m
.... ....
~
~
m
m
~
.4
(bl
FIG. 2. Synthetic planewave seismograms for a tenlayer model. The model has two steps in the velocity profile and one step in the density profile. This is referred to as Set 2 in the text.
FIG. 3. (a) The Green's function: The Green's function generated from a 36microlayer model which is a refined model of Set 1 at two different p values (0 and 0.636) modulated by the finite frequency band. By "refined" model, we mean that the thicknesses of the microlayers in Set 1 are redefined to have equal traveltime (4 ms) thicknesses of  6 m. The first panel is the Green's function of the normalincidence plane wave. In the second panel, the large p value is chosen to investigate the phenomena of nearly critical and critical reflection waves, refraction waves, and evanescent waves. From layer 1 to layer 5 where c = 1500 mis, the wave is traveling at an angle very close to the critical angle and the downgoing and upreflected waves are clearly visible. From layer 6 to layer 20, where c = 1550 mls, the wave passes from a critical reflection wave, to a refraction wave, to an evanescent wave. Below layer 20 where c = 1600 m/s, the decaying evanescent wave with depth is notable. The amplitudes are true magnitude in both of the panels. (b) Frechet derivatives with respect to bulk modulus at the same two p values from analytical equation (22) for the same 36microlayer model as in (a). The two panels are subjected to different scaling factors due to the fact that the amplitude of the right panel, where p = 0.636, is several orders greater than that of the left panel, where p = o.
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
MultioUset Fullwaveform Inversion
27
phase matrices En may be generated recursively through the sequence of frequencies (0, 8ro, 28ro, ...). Finally, the computation of synthetic seismograms by means of equation (16) is done by fixing p and running a sequence of frequencies as the innermost loop. This process permits vectorization and implementation on an array processor or supercomputer, where the frequency series is the basis for vectorization (Phinney et a1., 1985). Given modern multimegabyte computers, the requirement for storage of large complex matrices (dimensioned as the number of frequencies times the number of microlayers) is not a constraint.
eling program once for each microlayer and would be an intolerable computational burden in the absence of an analytic theory.
NUMERICAL EXAMPLES
Validation of the calculations In this application it is critical that the Green's functions not have any hidden errors of factors of 2 or n. The algebra and the code were checked against the case of a single layer embedded in a whole space, for which the analytical solution is easily written. In addition, the Frechet kernels, equation (20),were numerically tested by evaluating the finitedifference approximation
F[P(t),
K ]
j
~ ~P(t) ~Zj'
~Kj
(27)
which represents the column of F which displays the effect of a perturbation at depth Zj' Figure 3b shows a sample panel of Frechet kernels calculated for a 36microlayer model from the analytic expressions and verified by finitedifference evaluation. The latter method requires running the forward mod
An integrated code has been developed which combines the elements described in the previous sections. The observed conventional xt (or rt) data set has to be transformed into the p? domain. This collection of planewave traces, taken as a function defined over the domain of (p, r), can be shown (Phinney et a1., 1981) to be related to the collection of conventional waveforms, taken as a function defined by range and time (r, t), by the Radon transform and its inverse. Under certain conditions, the Radon transform can be evaluated by a simple "slant stack" in which a data point in (p, t) is generated by summing the (r, t) domain over a slant line with slope p and intercept time t. The form appropriate for cylindrical coordinates requires additional operations on the traces during transformation, and is considerably more of a computational burden (Chapman, 1981). In this study the synthetic test data were generated directly in the p: domain by assuming that the proper Radon transform had been applied to an adequate conventional rt data set. Hence the inversion can be restated as the following problem: Given a few traces of a px seismic section, invert for the density and velocity (bulk modulus) profiles. The procedure is summarized in Figure 4. It makes use of waveform matching: the "observed" planewave traces are compared with traces calculated from a model using a synthetic seismogram program. The functions of Z (velocity and density) which describe the model are adjusted in a series
PLANE WAVE REPRESENTATION
'1'(Pj,'T)
CONVENTIONAL REPRESENTATION ep(f;,t) OBSERVED DATA
1= 1 to 480
ACOUSTIC:
K=2
ELASTIC:
TRANSFORMED DATA
J = 1 to 5
RADON TRANSFORM
K=3
UPDATE MODEL ITERATIVE LOOP
COMPARE and ITERATE
COMPLETE tl~
COMPUTE COVARIANCE MATRIX
MEDIUM ELASTIC MODULI and DENSITY
COMPUTE SYNTHETIC DATA
SYNTHETIC TRANSFORMED DATA
FIG. 4. Flow chart of iterative nonlinear leastsquares inversion in the planewave domain. By using the planewave version of the data, the computationally intensive iterations avoid the Radon transform step. Also, the effort required to calculate synthetic data is minimized, since each p needs a separate pass. This inversion will run on a minicomputer. It requires input data which are densely sampled in the range r.
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
28
Pan et 81.
of iterations until a satisfactory agreement between the observed and synthetic data is achieved. Two sets of synthetic seismograms were generated in the pot domain to be used as observed data. One, which is referred to as Set 1 (Figure 5b), is from a tenmicrolayer model that has a constant density profile and two steps in the velocity profile (Figure 5a). The other, which is referred to as Set 2 (Figure 2b), is from another tenlayer model that has one step in the density profile and the same velocity profile as the previous one (Figure 2a). The models are specified by tabulating the compressional velocity and density functions at 20 m intervals, which define the microlayers of the model. The source time function (which could be treated as an unknown for the inversion) is a specified bandlimited minimumphase transient wavelet. Inversions were performed using traces selected from Set 1 and Set 2. In the first inversion, the first trace (p = 0.00) and the eleventh trace (p = 0.25) from Set 1 (Figure 5; also a magnified version appears in Figure 6a, panel B) are used as input" observed" data. The inversion parameters are chosen as follows: (1) The model covariance matrix Cmo is a tridiagonal matrix: 1 x 1018 along the diagonal and 1 x 101 7 along offdiagonals for the bulk modulus (velocity) model; 1 x 103 along the diagonal and 1 x 102 along offdiagonals for the density model. [Note that the units for the modulus covariance matrix are (kg/ms2)2 and the units for the density covariance matrix are (kg/m 3)2.] The model covariance matrix is scaled so that its diago
nal elements are unity in the hybrid units {[(g/cm3)(km2/s2)] 2}. Thus, (a) the maximum variation between the a priori bulk modulus (velocity) model and the inverted bulk modulus model is 1 if the velocity and density are in (km/s) and (g/cm3 ) ; (b) the variance of the a priori density model covariance is very small (i.e., we assume we have very good a priori information about the density function), and (c)the properties of each layer are weakly correlated with the properties of the layers immediately above and below it. (2) The data covariance matrix is a Cdo diagonal matrix with 0.025 along the diagonal. This implies that each data point is allowed to have 5 percent noise and the noise of each data point is uncorrelated with that of any other point. (3) The \y operator we implemented is [diag (! + CmE[CdEk)r 1 only, to prevent the need to invert the full \Y matrix at each iteration. (Note the Frechet derivative Ek is updated at each iteration.) (4) Each iteration, the model m, is updated with the result from the previous iteration.
55
~
L
0.75 1.00
DENSITY jglcm3J
(a)
I
I
.1
1.25
1.50
1.75
.10'
VELOCITY (mi.]
SLOWNESS P
.000
.0
,125
2SO
.375
.500
?> .;.
\slkmJ
.4
I
(b)
FIG. 5. Synthetic seismograms in the pot domain for a tenlayer model embedded between two halfspaces. The model has a constant density profile and two steps in the velocity profile. This is referred to as Set 1 in the text.
In Figure 6a we show the two planewave traces: (A) calculated from the starting (first guess) model, (B) for the observed data, (C) calculated from the final inverted model, and (D) for the residual (B) minus (C). The model velocity and density profiles of the starting model, the "true" model, and the final inverted model are shown as dashed lines, dotted lines, and solid lines, respectively, in Figure 7a. The relative error, defined as the ratio of the power of residual signals to the power of observed signals, is shown in Figure 8a versus the number of iterations . The velocity and density profiles are recovered with fairly good accuracy in this case. Ninetyfour iterations were needed in this study to reach the final inverted model. From equation (11), it is clear that a larger Cmo results in a faster search through the admissible model space and also implies that our information about the starting model is poorer. In real life, Cmo is unknown. The optimal choice of Cm is a point of o further study. The convergence rate is sped up by about a factor of four when we specify a model covariance four to five times as large as the one used. However, such a step degrades the quality of the inversion. To determine the optimal Cm for o inversion, more experiments with different parameters are needed, in order to obtain the empirical tradeoff curve between the quality of the inversion (the size of the relative error) and the convergence rate (the magnitudes of the model covariances) of the inversion. The magnitude of Cmo which is related to the step size, or in a sense corresponds to the length of the partial (inexact) line search (Fletcher, 1980), determines the search distance of each iteration in the model space. In the second inversion, all the parameters are the same as in the first one except that (a) the model covariance matrix has values 5 x 101 7 and 2 x 105 along the diagonal of the bulk modulus model and the density model, respectively; (b) the offdiagonals are assigned to be onetenth of the diagonal
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
MUltioffset Fullwaveform Inversion
29
(p
values; and (c) the two input traces are taken from the first = 0.00) and the eleventh (p = 0.25) traces of Set 2. Figure 6b shows, from (A) to (0), the starting data, the observed data, the inverted model data, and the error signal, as in Figure 6a. The model profiles of this experiment are shown at the bottom of Figure 7 with the same notations used at the top. The relative error curve is shown in the dotted line in Figure 8b. The error signals in Figure 6b, panel 0, appear to be substantially larger than those of the first experiment (Figure 6a, panel 0). We attribute the error mainly to the inaccuracy of the velocity inversion (Figure 7b), i.e., due to the inexactness of phase information of the seismogram instead of a discrepancy in the amplitude of the general waveforms. In comparison to
A.
. 0,,
B.
c.
D.
P
rs/krm
the first experiment, the accuracy of this inversion in terms of minimum relative error is degraded (Figure 8a and the dotted line, Figure 8b). In the previous inversion we used two input traces to invert for one parameter only (the velocity profile); in this inversion we used two input traces to invert for two parameters simultaneously (the velocity and density profiles). Shown in Figure 9, the acoustic impedance profile, the product of the velocity and density, is better constrained in this inversion. Whenever the inverted velocity is lower than the "true" velocity, the inverted density turns out to be higher than the "true" density, and vice versa. This relationship tends to produce the proper impedance change over the interfaces. In an attempt to improve the quality of the second experiment, an inversion was performed using four traces, the first, sixth, eleventh, and sixteenth traces (p = 0.000, 0.125, 0.250, and 0.375) from Set 2 and the same inversion parameters. The inversion converged faster and is more accurate, as expected, because more input traces provided greater redundancy for the inversion procedure. In the same notation as before, the planewave traces and the model profiles for the last experiment are shown in Figures 10 and 11, respectively. The relative error curve is shown in the vertical bars in Figure 8b. The second inversion is illconditioned with respect to the lower part of the model; by increasing the number of data traces, a
,4+++++++'(a)
0,,,.
55
110
165
220u......~~.......J._~~
........ 1.25
(a)
_'_'_'_
_'"
0.75
A.
ll .
C.
D.
P
rs/krm
1.00
DENSITY [g/cm 3 ]
.0
1T
:
0
"1
55
I
co
!
t
110
165
A+++1++1+220
I
II:
~.'i
, t!
I
I
:
:
(b)
0.75
1.00
DENSITY [g/cm3J
1.25
(b)
1.50
1.75 x 10
3
VELOCITY (m/s]
FIG. 6. (a) Inversion of two traces from Set 1: synthetic seismograms of (A) the starting model; (B) the "true" model (trace 1 and trace 11); (C) the final inverted model; and (0) the residual traces (B)(C). (b) Inversion of two traces from Set 2: the two planewave traces (A) for the starting model; (B) from the "observed data;" (C) for the best deduced model; and (D) for the error signals, i.e., traces (B)~(C).
FIG. 7. (a) Inversion of two traces from Set 1. The final inverted model is shown in the solid lines, the "true" model and starting model are shown by dotted lines and dashed lines, respectively. (b) Inversion of two traces from Set 2: the starting model, the true model, and the best deduced model are in dashed lines, dotted lines, and solid lines, respectively.
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
30
Pan et at
substantial improvement has been obtained. We remark that the illconditioning which is more or less inherent in any inversion procedure becomes noticeable deeper in the model. Consequently; the synthetic cases we have considered cannot be considered noisefree, although no simulated seismic or instrument noise has been added. Numerical errors are the limiting noise source in this case. In the companion paper (Pan and Phinney, 1987), where noise is explicitly added, the device
0
55
110
165
220 1.00 1.25 1.50
r
1.75 2.00
x
106
impedance (0)
100:1:
FIG. 9. Inversion of two traces from Set 2. The impedance profile is reconstructed with fairly good accuracy (d. Figure 7b).
75
D.
P (slkm)
I .
It
50
.ill
;::
".
25
,
.....
...... ...........................
50
75 100
)
oo......~~~~''~~~~~~~~~''~~~~....u
o
25
NUMBER OF ITERATIONS
(a)
FIG. 10. Inversion of four traces from Set 2: synthetic seismograms of (A) the starting model; (B) the "true" model (trace 1, trace 6, trace 11, and trace 16); (C) the final inverted model; and (D) the residual traces (B)(C).
100:1:
O.,~.,,
L
55
75
iIt
.. ..
It
..
50
'
:[
Ii: w
J:
110
Cl
;::
It
ill
25
165
220c.......~~~ .....................~~....L~~~....l.''
.......... ~........u
0.75
oo......~~~~''~~~~~~~~~''~~~"'"'....u
1.00
DENSITY Ig/cm'l
1.25
1.50
VELOCITY (m/s)
1.75 x 10
3
o
25
75
100
NUMBER OF ITERATIONS
FIG. 11. Inversion of four traces from Set 2. The starting model, the true model, and the final inverted model are shown in dashed lines, dotted lines, and solid lines, respectively.
FIG. 8. The relative error curve versus the number of iterations. (a) Inversion of two traces from Set 1. (b) Inversion of two traces (dotted line) and four traces (vertical bars) from Set 2.
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/
Multioffset Fullwaveform Inversion
31
of using additional data traces also successfully stabilizes the inversion. CONCLUSIONS We have cast the general inversion framework of Tarantola and Valette (1982a) into a form applicable to the strictly layered problem and implemented it as a computer program. The initial demonstrations illustrate both good conditioning and ill conditioning and suggest that the flexibility of the technique should make it quite useful in assessing the information content in a data set. In a companion paper, we have reported on the inversion of more realistic microlayered media having the appearance of a typical sonic log and explore the variance and resolving power of these examples, with the inclusion of seismic noise. Further extensions of this work would include (1) implementation for layered elastic media; (2) implementation for conventional (rt) source gathers; (3) exploration of different parameterizations of the model, including adaptive reparameterization; (4) a singular value decomposition of data and model spaces, to look deeper into conditioning; (5) experiments with improving robustness, as in casting out data points which lie on the flat tail of the error distribution; and (6) study of sequences of inversions in which data spaces and model spaces are gradually expanded. For practical application, the method looks quite promising for marine data, where shortwavelength statics are generally absent and where the commonmidpoint gather is immune to the effects of gradual lateral variation in structure. With land data, residual statics remain a limiting factor and are a part of the inverse problem. With sufficient computing power, it is conceivable that residual statics and layering might be jointly inferred from adequately redundant field data. Note that amplitude statics as well as time statics are a part of the problem.
ACKNOWLEDGMENTS This research was partially supported by the Office of Naval Research under contract NOOOI484COI51, and by research funds of the Department of Geological and Geophysical Sciences, Princeton University. Pan would like to thank Dr. K. RoyChowdhury for his constant support, encouragement, and tremendous amount of time in helping to debug the computer programs. REFERENCES Backus, G., and Gilbert, F., 1967, Numerical applications of a formalism for geophysical inverse problems: Geophys. J. Roy. Astr. Soc., 13,247 276.    1968, The resolving power of gross earth data: Geophys. 1. Roy. Astr. Soc. 16, 169205.    1970, Uniqueness in the inversion of inaccurate gross earth data: Phil. Trans. Roy. Soc. London A, 266,123192.
Berryman, 1. G., 1979, Inverse methods for elastic waves in stratified media: J. Appl. Phys., SO, 67426744. Berryman, J. G., and Greene, R. R., 1978, Discrete inverse scattering theory and the continuum limit: Phys. Lett., 65A, 1315. Carrion, Ph.M., Kuo, J. T., and Stoffa, P. L., 1984, Inversion method in the slant stack domain using amplitudes of reflection arrivals: Geophys. Prosp., 32, 375391. Chapman, C. H., 1978a, A new method for computing synthetic seismograms: Geophys. J. Roy. Astr. Soc., 54, 481518.    1978b, Body waves in seismology, in Ackerback, J. D., Ed., Modern problems in elastic wave propagation: John Wiley and Sons. 1981, Generalized Radon transforms and slant stacks: GeophysJ. Roy. Astr. Soc., 66, 445454. Clayton, R. W., and McMechan, G. A., 1981, Inversion of refraction data by wavefield continuation: Geophysics, 46,869874. Clayton, R. W., and Stolt, R. H., 1981, A BornWKBJ inversion method for acoustic reflection data: Geophysics, 46,15591567. Deans, S. R., 1983, The Radon transform and some of its applications: John Wiley and Sons. Fletcher, R., 1980,Practical methods of optimization: John Wiley and Sons, Inc. Franklin, 1. N., 1970, Wellposed stochastic extensions of illposed linear problems: J. Math. Anal. Appl., 31, 682716. Fuchs, K., and Muller, G., 1971, Computation of synthetic seismograms with the reflectivitymethod and comparison of observations: Geophys. J. Roy. Astr. Soc., 23, 417433. Gilbert, F., 1971, Ranking and winnowing gross earth data for inversion and resolution: Geophys. 1. Roy. Astr. Soc., 23, 125128. Jackson, D. D., 1972, Interpretation of inaccurate, insufficient and inconsistent data: Geophys. 1. Roy. Astr. Soc., 28, 97110.  1973, Marginal solutions of quasilinear inverse problems in geophysics: the hedgehog method: Geophys. 1. Roy. Astr. Soc., 35, 121 136.    1979, The use of a priori data to resolve nonuniqueness in linear inversion: Geophys. 1. Roy. Astr. Soc., 57,137157. Jordan, T. H., and Franklin, J. N., 1971, Optimal solutions to a linear inverse problem in geophysics: Proc., Nat. Acad. Sci.,68, 291293. Nilsen, A., and Gjevik, B., 1978, Inversion of reflection data: Geophys. Prosp., 26,421432. Pan, G. S., and Phinney, R. A., 1987, Full waveform inversion of planewave seismograms in stratified acoustic media: applicability and limitations: submitted to Geophysics. Phinney, R. A., adom, R. I., and Fryer, G., 1985, Rapid generation of synthetic seismograms in layered media by vectorization of the algorithm: Internal research report, Dept. of Geol. and Geophys. Sci., Princeton Univ. Phinney, R. A., RoyChowdhury, K., and Frazer, L. N., 1981, Transformation and analysis of record sections: 1. Geophys. Res., 86, 359377. Stoffa, P. L., Buhl, P., Diebold, J. B., and Wenzel, F., 1981, Direct mapping of seismic data to the domain of intercept time and ray parameterA planewave decomposition: Geophysics, 46, 255267. Tarantola, A., 1984a, Linearized inversion of seismic reflection data: Geophys. Prosp., 32, 9981015.    1984b, Inversion of seismic reflection data in the acoustic approximation: Geophysics, 49,12591266. Tarantola, A., and Valette, B., 1982a, Generalized nonlinear inverse problems solved using the leastsquares criterion: Rev. Geophys. Space Phys., 20, 219232.    1982b, Inverse problems  quest for information: Geophys. J. Roy. Astr. Soc., SO, 159170. Treitel, S., Gutowski, P. R., and Wagner, D. E., 1982, Planewave decomposition of seismograms: Geophysics, 47, 13751401. Ware, J. A., and Aki, K., 1969, Continuous and discrete inverse scattering problems in a stratified elastic medium: J. Acoust. Soc. Am., 45,911921. Wiggins, R., 1972,The general linear inverse problem: implications of free oscillations and surface waves for earth structure: Rev. Geophys. Space Phys., 10, 251285. Wiggins, R., McMechan, G., and Toksoz, M., 1973, Range of Earthstructure nonuniqueness implied by body wave observations: Rev. Geophys. Space Phys., 11,87113.
Downloaded 25 Aug 2011 to 99.172.42.203. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/