GeoScienceWorld
Volume

An Analysis of Least-Squares Velocity Inversion

By Fadil Santosa and William W. Symes
Edited by Raymon L. Brown

Abstract

“This book grew out of an attempt to understand the mechanisms through which band-limited reflection seismograms determine velocity distributions in elastic models of the earth's crust. The authors were especially interested in the feasibility of recovering very slowly varying (out-of-passband) velocity components from band-limited (highfrequency) reflection data. That interest was spurred by reports of successful inversions for layered media.”

  1. Page 1
    Abstract

    The output least-squares approach to inverse problems in seismology has attracted a great deal of attention in recent years. Also known as least-squares inversion, nonlinear iterative inversion, generalized linear inversion, and a variety of other names, it involves a systematic search for an earth model in some class which best fits some type of seismic data in a least-squares (root-mean-squares, rms) sense. Recent contributors include Bamberger et al. (1979 and 1982), Tarantola and Valette (1982a and 1982b), Lesselier (1982), Keys (1983), Tarantola (1984 and 1986), Lailly (1984), McAulay (1985 and 1986), Gauthier et al. (1986), Kolb et al. (1986), Canadas and Kolb (1986), Mora (1987a and 1987b), Chapman and Orcutt (1985), Shaw and Orcutt (1985), and Pan et al. (1988). For a lucid discussion and many older references consult Lines and Treitel (1984). Further discussion may be found in a new monograph by Tarantola (1987).

    There seems to exist considerable confusion regarding the sort of information about the subsurface which one might expect to extract using the least-squares approach to the inverse problem of reflection seismology, and also regarding the quality of least-squares inversion results relative to the output of conventional processing methods. The relation between least-squares inversion and migration of both stacked and unstacked data is now well understood

  2. Page 9
    Abstract

    We consider a linearly elastic isotropic medium confined to the half-space {z > 0}, subject to prescribed traction on {z = 0}. We introduce the notation

    1. x = position vector = (x1, x2, x3) or = (x, y, z) interchangeably;

    2. u = displacement vector;

    3. ρ = material density;

    4. λ, μ = Lamé parameters;

    5. τij = λ∇ ·ij + μ(∇jui + ∇iuj) = + μ(ui,j + uj,i)

    6. = (i,j)-component of stress

    7. f = source wavelet (time function).

    Then the response to a horizontally polarized shear line load with wavelet f, extended in the y− (x2−) direction, is the solution of the initial boundary value problem consisting of the equations of motion together with boundary conditions and initial conditions

    Remark. We have chosen to study the line-load case, despite certain unrealistic aspects, because:

    (i) We can, in this way, isolate the scalar SH-wave problem, which is a simpler framework in which to present the main ideas;

    (ii) It is much less expensite to simulate 2-D (i.e. line-source) seismograms and associated quantities, for the purpose of numerical illustration.

    No obstacle exists, in principle, to the extension of the results presented here to the point-load case.

    The SH line-source surface response (seismogramid, ealized common source gather) is the quantity

    To make explicit the dependence of s on the elastica nd sourcep arameters, we shall also write s(x, t) = s[p, μ , f, x0](x, t). Note that, since the coefficients of the wave equation are independent of x

  3. Page 15
    Abstract

    The most primitive version of the seismic inverse problem, specialized to the layered, μ = const. SH line-source model, is the functional equation (3.1) The right hand side of (3.1) is meant to represent the Radon transform of the measured line source SH-velocity surface trace, and we demand that the equation hold at least in the precritical part P of the (p, τ) domain.

    Generally, equation (3.1) is overdetermined, and will have no solutions c at all if G is produced by slant-stacking a noisy data set. Data G for which a solution c of equation (3.1) exists are called consistent. Generic data are inconsistent. Therefore we must modify the statement of the inverse problem.

    Before stating alternative formulations, we should make explicit the goal which any reformulation is to achieve. A modest objective might be the stability of the solution for near-consistent data. We could reasonably hope that any reasonable data set is near a consistent data set: otherwise, our model of the basic physics is erroneous. Thus suppose that

    where

    for somer easonablev elocity profile c*, and N is a “small” noiset erm (we shall be more explicit about the meaning of “small” later). Then we should consider a restatement of the inverse problem successful if it enables us to derive from G a velocity estimate c which is close to c*.

    A second useful quality of a formulation of the inverse problem would be that computation of its solution

  4. Page 17
    Abstract

    As mentioned in the introduction, a popular approach to the problem of inconsistent data is to seek a model (i.e., velocity profile c) for which the mean-square error(4.1) is as small as possible. This is the output least-squares formulation of the inverse problem, specialized to the precritical p-tau section.

    Actually an even more popular approach is to attempt to minimize the mean-square error in the full (x, t) seismogram: (4.2) The error defined in equation (4.1) is essentially identical to the error in the point source seismogram (see Appendix B; the extra “p” in the integral in (4.1) is not a misprint, but is the correct weight to bring equation (4.1) and the point source seismogram error as close as possible). The relation of the point-source seismogram error with the line source seismogram error (4.2) is more complicated. Also it is possible to introduce weights (“data covariance matrix”) to reflect the presumed structure of data errors. See Tarantola and Valette (1982a and 1982b) for details.

    As we will explain, the minimization of equation (4.1) is easier than the minimizationo f equation( 4.2), so we shallc oncentratem ostlyo n equation (4.1).

    Unfortunately, the solution of even this “precritical” output-least-squares problem is excessively sensitive to data noise, for several reasons. This is so even in the “perfect” case of impulsive data, i.e., f(t) = δ(t). In that case, a stable estimate of c

  5. Page 21
    Abstract

    In this chapter we examine the sensitivity of the least-squares inversion to data noise or perturbation. We include this well-known material for the sake of completeness; it has been treated with great care, for instance, in Lines and Treitel (1984). In the notation of the previous chapters, we want to determine the effect on a (the) minimum of caused by a perturbation δG in the data G. We shall use here the standard notation from advanced calculus “D” for derivative: as applied to the plane-wave seismogram S for a reference velocity c and a perturbation δc, Thus DS[cc (“the directional derivative of S at c in the direction δc”) is exactly the first order perturbation in S due to the perturbation δc in velocity. DS[c] itself is a linear operator, mapping velocity perturbations δc to corresponding perturbations in S.

    In Chapter 2, we definedt he formal first orders eismogramp erturbation δS as the solutiono f the perturbationabl oundaryv aluep roblem( 2.6). If the limit defining DS[c]δe xists,t hen it oughtt o be true that DS[c]δc = δS. Circumstancesu nder which this relation holdsa re discussedin Symes (1986a)a nd Sanrosaa nd Symes( 1988b).

    Now J is a real valued function of two arguments (c and G), hence has a gradient with respect to each. Let 〈 , 〉e denote a suitable inner product for velocity perturbations, possibly incorporating some weighting. A simple (but not the only) choice is

    Then the c-gradient is that

  6. Page 25
    Abstract

    The singular spectrum of the normal operator DS[c]*DS[c] is easiest to analyze in case the reference velocity c is slowly varying, i.e., smooth. In fact, this analysis underlies much of the heuristic reasoning in reflection seismology.

    We temporarily restrict our attention to a single plane-wave component, which may as well be the normal incidence component (p = 0), and write For constant reference velocity c = const., the derivative of S0 may be computed in closed form: This is the famous “Born approximation,” or convolutional model, which forms the basis of much seismic data processing. See for example Waters (1981), Cohen and Bleistein (1979), and Gray (1980). The following is also a special case of the formula derived in Appendix C: (6.1) where is the reparameterization of δc by two-way time.

    The singular spectrum of the normal operator DS[c]*DS[c] is easiest to analyze in case the reference velocity c is slowly varying, i.e., smooth. In fact, this analysis underlies much of the heuristic reasoning in reflection seismology.

    We temporarily restrict our attention to a single plane-wave component, which may as well be the normal incidence component (p = 0), and write For constant reference velocity c = const., the derivative of S0 may be computed in closed form: This is the famous “Born approximation,” or convolutional model, which forms the basis of much seismic data processing. See for example Waters (1981), Cohen

  7. Page 41
    Abstract

    In this chapter we show that the presence of interfaces, or even near-discontinuities, of significant magnitude in the reference velocity model c(z) completely changes the character of the singular spectrum of DS[c].

    Most of the spectral change is evident in the simplest example with an interface; the “layer-over-a-half-space” configuration

    In this chapter we show that the presence of interfaces, or even near-discontinuities, of significant magnitude in the reference velocity model c(z) completely changes the character of the singular spectrum of DS[c].

    Most of the spectral change is evident in the simplest example with an interface; the “layer-over-a-half-space” configuration

    In this chapter we show that the presence of interfaces, or even near-discontinuities, of significant magnitude in the reference velocity model c(z) completely changes the character of the singular spectrum of DS[c].

    Most of the spectral change is evident in the simplest example with an interface; the “layer-over-a-half-space” configuration

    In this chapter we show that the presence of interfaces, or even near-discontinuities, of significant magnitude in the reference velocity model c(z) completely changes the character of the singular spectrum of DS[c].

    Most of the spectral change is evident in the simplest example with an interface; the “layer-over-a-half-space” configuration

    In this chapter we show that the presence of interfaces, or even near-discontinuities, of significant magnitude in the reference velocity model c(z) completely changes the character of the singular spectrum of DS[c].

    Most of the spectral change is evident in the simplest example with an interface; the “layer-over-a-half-space” configuration

    In this

  8. Page 65
    Abstract

    In Chapter 5, we saw that the first goal of any formulation of the inverse problem, i.e., that the solution not be inordinately sensitive to data noise, can be met by the least-squares formulation provided that the smallest singular value of the linearized seismogram is not too small. In Chapters 6 and 7 we saw that while this smallest singular value is very small indeed for smoothly varying models, it may be quite a bit bigger for “rough” models, i.e., models with a sufficiently dense set of strong reflectors, provided that a sufficiently wide range of precritical slownesses are incorporated in the definition of the (plane wave) seismogram. Since a rough, reflector-rich, model is a more accurate approximation of (at least some) sedimentary columns than is a smoothly varying velocity profile, one might be mildly optimistic regarding the stability property.

    In this chapter we explain the implications of our analysis for the prospects of achieving the second goal mentioned in Chapter 3, namely efficiency of computation.

    In Chapter 5, we saw that the first goal of any formulation of the inverse problem, i.e., that the solution not be inordinately sensitive to data noise, can be met by the least-squares formulation provided that the smallest singular value of the linearized seismogram is not too small. In Chapters 6 and 7 we saw that while this smallest singular value is very small indeed for smoothly varying models, it

  9. Page 83
    Abstract

    The general features of nonlinear least-squares solvers suitable for use in the seismic inverse problem have been well explained in the numerous articles advocating the output-least-squares approach: see for instance Bamberger et al. (1979), Tarantola (1984), Lailly (1984), and Lines and Treitel (1984).

    Accordingly, we shall describe the implementation of our least-squares algorithm only in general terms. Such algorithms require (approximate) computation of the reference, perturbational, and adjoint fields. For the model explained in Chapters 2-4, are determined by the initial/boundary value problems (9.1)(9.2)

    The general features of nonlinear least-squares solvers suitable for use in the seismic inverse problem have been well explained in the numerous articles advocating the output-least-squares approach: see for instance Bamberger et al. (1979), Tarantola (1984), Lailly (1984), and Lines and Treitel (1984).

    Accordingly, we shall describe the implementation of our least-squares algorithm only in general terms. Such algorithms require (approximate) computation of the reference, perturbational, and adjoint fields. For the model explained in Chapters 2-4, are determined by the initial/boundary value problems

    The general features of nonlinear least-squares solvers suitable for use in the seismic inverse problem have been well explained in the numerous articles advocating the output-least-squares approach: see for instance Bamberger et al. (1979), Tarantola (1984), Lailly (1984), and Lines and Treitel (1984).

    Accordingly, we shall describe the implementation of our least-squares algorithm only in

  10. Page 91
    Abstract

    We conducted a number of numerical experiments with the code described in Chapter 9. We will detail here the results of several such tests.

    In all cases, we produced the data (p-tau section) by solving the system (9.1) numerically. As noted in Chapter 2, the assumption that plane-wave seismogram components are modeled by the plane-wave equations is only approximately true. Still, the approximation seems to be quite good, and the distinction has been ignored with apparent impunity in work on modeling (by the reflectivity method), migration and velocity analysis, and inversion. This gives us some hope that the results presented here will accord well with future tests involving precriticai projection of the full 2-D seismogram.

    We should also note that the data were produced by applying the finite-difference algorithm to models defined on somewhat finer grids than the nodal grids for the splines used to represent the reflectivity estimates in the quasi-Newton code. Thus there is no possibility that the iterative process can achieve zero residual seismogram, although the spline grid is chosen sufficiently fine (relative to the spatial wavelengths present in the fields) to allow a reasonably good fit.

    We conducted a number of numerical experiments with the code described in Chapter 9. We will detail here the results of several such tests.

    In all cases, we produced the data (p-tau section) by solving the system (9.1) numerically. As noted in Chapter 2, the assumption that plane-wave seismogram components are modeled

  11. Page 105
    Abstract

    In the foregoing chapters, we have explained both the success of least-squares inversion in extracting both passband and subpassband information about the earth model from band-limited data, and the computational difficulty inherent in the approach, in the context of the layered velocity model.

    The main point in our study has been to determine whether information about the trend of a velocity profile is contained in band-limited seismic data. We emphasize here the importance of the smoothly varying (low frequency) component of the velocity in seismic reconstruction. As mentioned in the introduction, not only are reflector locations misplaced when the trend is incorrect, but also the strength of the reflectors are miscalculated. This complication is even more severe in the case of multiparameter (soundspeed and impedance) inversion.

    We have shown that when the model attains a certain degree of “roughness”, i.e., has a sufficiently dense set of reflectors, smooth perturbations are associated to large changes in the band-limited seismogram, because they induce traveltime changes. Consequently the singular spectrum of the perturbational seismogram about a rough velocity profile has a completely different character than that about a smooth profile: in the latter case, smooth perturbations are associated with small singular values. Additionally, we have arguedt hat whent he slownesrsa nge[ i.e., the rangeo f c(z)p] is sufficiently large at each depth of interest, no very small singular values occur at all, and the linearized least-squares problem is safely convex. Thus velocity perturbations are well-determined down to a length scale determined by the

  12. Page 109

Purchase Chapters

Recommended Reading