Multiparameter Inversion: Cramer's Rule for Pseudodifferential Operators

Abstract

Linearized multiparameter inversion is a model-driven variant of amplitude-versus-offset analysis, which seeks to separately account for the influences of several model parameters on the seismic response. Previous approaches to this class of problems have included geometric optics-based (Kirchhoff, GRT) inversion and iterative methods suitable for large linear systems. In this paper, we suggest an approach based on the mathematical nature of the normal operator of linearized inversion—it is a scaling operator in phase space—and on a very old idea from linear algebra, namely, Cramer's rule for computing the inverse of a matrix. The approximate solution of the linearized multiparameter problem so produced involves no ray theory computations. It may be sufficiently accurate for some purposes; for others, it can serve as a preconditioner to enhance the convergence of standard iterative methods.

1. Introduction

The linearized inverse problem in reflection seismology aims at recovering model perturbations from data perturbations, assuming known reference or background model. Both reference model and model perturbations consist of material parameter fields (functions of position x) that characterize wave propagation. As the background model is considered known, we will suppress it from the notation and use the word “model” and the symbol 𝑚 and similar symbols to refer to model perturbations. The linearized scattering operator 𝐹 maps a model perturbation 𝑚 to the corresponding perturbation 𝐹𝑚 of the predicted data.

With these conventions, the inverse problem may be stated as follows: given observed data 𝑑, find a model perturbation 𝑚 so that 𝐹𝑚≈𝑑. (1) Interpreting (1) in a least squares sense yields the normal equations𝐹^∗𝐹𝑚=𝐹^∗𝑑. (2)𝐹^∗ is the adjoint (transpose) of 𝐹: it is a migration operator (to be precise, the result of reverse time migration, properly construed). The migration operator maps the data space to model space. We will refer to the migration output 𝑚_mig=𝐹^∗𝑑 as the migrated image, acknowledging its typical role in seismic processing. 𝑁∶=𝐹^∗𝐹 is the normal operator, or Hessian (we will use these terms interchangeably).

An example of this setup is linear acoustics; we will treat this example explicitly in this paper. The model 𝑚 consists of the velocity and density (or impedance and density, or any other equivalent combination) perturbation fields 𝑚=(𝛿𝑣,𝛿𝜌), for instance. The forward map samples the solution of the linearized (Born approximation) acoustic wave equation for a source (right-hand side) indexed by source position, at receiver positions, possibly with additional filtering to emulate antenna patterns induced by source and receiver arrays.

In simulating seismic data generation with models such as linear acoustics, waves typically propagate over hundreds of wavelengths, and model fields must resolve features on the wavelength scale. The normal equations thus represent millions, or billions, of equations in the same number of unknowns, ruling out the possibility of solving them by means of direct matrix methods such as Gaussian elimination. This paper presumes that the computations implementing the application of the normal operator are carried out by “wave equation” methods, that is, finite difference or finite element simulation. Thus each application is an expensive, large-scale computational procedure, a fact which places practical limits on the number of steps taken in an iterative scheme to solve the normal equations (2). For example, Krylov subspace methods (such as the conjugate gradient method) require at least one application of the normal operator per iteration.

We present an efficient method to approximate the solution of the normal equations that requires a few applications of the normal operator (one for one-parameter inversion, two for two-parameter inversion, and six for three-parameter inversion). This method leverages the properties of the normal operator: under some conditions (background model parameter fields slowly varying on the wavelength scale, diving wave energy eliminated from data, data polarized by propagating phases), it is a special type of matrix-valued spatially varying filter [1–9]. This particular type of spatially varying filter, called a pseudodifferential operator in the mathematical literature [10], preserves the discontinuities (events) in the model 𝑚 to which it is applied while effecting a dip- and frequency-dependent scaling of the amplitudes. Its application also mixes events when several model parameters (density, 𝑝- and 𝑠-wave velocities, etc.) are present: thus a density event will show up in the velocity component of acoustic migration output, and viceversa, for example.

In this work, we shall show how to separate the events corresponding to various model parameters by means of several applications of the normal operator to permuted image vectors, so that the result differs from an inversion of the data by an overall spatially dependent filter, common to all components. We have previously solved the problem of estimating and correcting for such a filter [11, 12]: it is the same problem as occurs in the solution of a single-parameter linearized inverse problem. Combining the two techniques, we recover an accurate approximate inversion for all parameters, within constraints which we shall illustrate.

We review our approach to single-parameter inversion (based on constant density acoustic modeling, for instance) in Appendix B. It belongs to the genre of scaling methods [13–16], which approximate the Hessian or its inverse by estimating a filter of some sort from a single application of the normal operator to a migrated image or some other test image. In previous work on single-parameter inversion, the second author pointed out that the pseudodifferential nature of the Hessian makes possible, a very simple, and accurate approximation, provided that the reflectivity dip is unique at every imaging location [16]. Our approach to single-parameter inversion without assuming uniqueness of dip (Appendix B) carries this idea further; it is similar to that described by Guilton [15], but with different constraints on filter structure motivated by the theory of the pseudodifferential normal operator cited above. The core of our method is the Pseudodifferential Operator Algorithm [17], an efficient algorithm for application of filters of the required type, which we review in Appendix A. Several other authors have explored similar pseudodifferential scaling methods based on other fast application algorithms [18, 19].

Note that data-adaptive scaling methods such as those cited in the previous paragraph differ in an essential way from scaling by an approximate diagonal of the Hessian [20, 21]. As noted by Rickett [14], a dip-independent scaling cannot well approximate an inverse Hessian in general, as the appropriate amplitude correction (from migrated image to inversion) depends on dip. On the other hand, these methods also differ from the so-called Kirchhoff or Generalized Radon Transform or Ray-Born inversion (e.g., [22–25]): both are asymptotic and approximate (either implicitly or explicitly) the normal operator as a pseudodifferential operator, but Generalized Radon Transform inversion involves explicit geometric optics computations, while scaling methods do not. Migration deconvolution methods [26] also approximate the Hessian, or its inverse, but differ from scaling methods in a similar way: migration deconvolution filters are constructed from approximate (layered) Green's function computations, rather than by data-adaptive fitting of an a priori limited filter class. Several authors have explored construction of localized filter representation of the Hessian or inverse Hessian using a more general class of Green's functions [27, 28]; again, scaling methods do not require precomputation of a Green's function database.

It has long been recognized that the offset (or scattering angle) dependence of a reflection event encodes material parameter changes across a reflector. Approximations to this relation motivated by the analysis of reflection and transmission in layered media (the Zoeppritz equations and linearization thereof, [29, 30]) lead to amplitude-versus-offset (“AVO”) analysis techniques [31]. Linearized inversion for one or more material parameters might be viewed as a quantitative counterpart, that is, as AVO inversion, and attacked via iterative linear system solvers [32–34]. Herrmann [18] uses a single-parameter inversion approach resembling ours to precondition iterative inversion for faster convergence; we envision a similar use of our techniques. Note that migration deconvolution has also been used as an aid to quantitative AVO [35]. As already mentioned, migration deconvolution relies on layered modeling to obtain approximate Green's functions. Pseudodifferential scaling is limited in other ways but does not require modeled wave dynamics used to construct the approximate Hessian to be close to those of layered media (and does not require the construction of Green's functions).

Several studies have analyzed the conditioning, or error propagation properties, of multiparameter inversion [34, 36–38]. We will not discuss this important issue except in that characteristic conditioning behavior will be evident in our examples: for linearized acoustic inversion from surface data, velocity or impedance perturbations are considerably better resolved than are density perturbations.

We have used only synthetic data in the work reported here and rather simple synthetic data at that. In fact, the very few published examples of inversion (in the sense of noise-level data fit) of exploration-scale data [39–41], strongly suggest that the success of linearized multiparameter inversion is sensitive to other aspects of the inverse problem: source pulse and radiation pattern estimation, background field approximation, modeling the physics accurately (including elasticity, attenuation, etc.), and to very careful data preprocessing. In this work, we chose to avoid these issues and focus only on the mathematical/computational issue of multiparameter inversion, given model-consistent data.

2. Theory and Methods

Recall that the aim of this paper is to solve the normal equations 𝑁𝑚=𝑚_mig=𝐹^∗𝑑 (3) for a model vector 𝑚=(𝑚₁,…,𝑚_𝑛) consisting of 𝑛 component material parameter perturbation fields 𝑚_𝑖. The migration output 𝑚_mig is a vector of the same type as the model, that is, an 𝑛-vector of components, which we shall call images, as they typically contain visual evidence of reflectors.

Since both input and output of the normal operator consist of 𝑛-vectors of material parameter fields, it is natural to represent 𝑁 as an 𝑛×𝑛 block matrix of operators 𝑁_𝑖𝑗 mapping material parameter fields of type 𝑖 to those of type 𝑗. Since material parameter fields are positive, it is possible to use relative perturbations to parametrize the model (thus, 𝛿𝜌/𝜌 for instance), in which, case 𝑁 and its blocks 𝑁_𝑖𝑗 are nondimensional. We chose not to nondimensionalize the problem as the techniques to be described below automatically produce results with appropriate units.

The Hessian may be treated as a matrix of spatially varying filters. The filter coefficient 𝐴_𝑖𝑗 defining the (𝑖,𝑗) component filter 𝑁_𝑖𝑗 depends on both the spatial position 𝑥 and the wavevector k, 𝑁_𝑖𝑗𝑚_𝑗(x)=𝐴_𝑖𝑗(x,𝑘)𝑚_𝑗(𝑘)𝑒^𝑖x⋅k𝑑k, (4) in which 𝑚_𝑗 denotes the spatial Fourier transform of the 𝑗th model component 𝑚_𝑗. Denote by 𝑁=op(𝐴) the operator defined by the matrix 𝐴=(𝐴_𝑖𝑗) of filter coefficients via (4).

Under certain conditions, the normal operator is a matrix of spatially varying filters of a special type, known as pseudodifferential, to be described below. By Beylkin [1] and Rakesh [2] established the first results of this type, which were systematically extended by Nolan and Symes [4], Ten Kroode et al. [6], De Hoop and Bleistein [5], Stolk [42], Stolk and De Hoop [43] and others. We summarize the theory as follows: the Hessian is well approximated by a pseudodifferential operator provided (i)the material parameters in the background model vary smoothly on the scale of a wavelength (since the theory is asymptotic in frequency, the technical assumption is that they are smooth, that is, infinitely differentiable, but the practical meaning is as stated here); (ii)diving wave energy is not present in the data or has been muted or dip-filtered out; (iii)the data has been polarized into propagating phase components.

These conditions are likely to be essential: either major reflectors in the background model or nonpolarized multiphase data lead to Hessians which produce nonphysical reflector images at some distance from their sources, hence cannot be well-approximated by their near-diagonal behavior. Since pseudodifferential operators are nearly local, in a sense to be made precise below, Hessians producing major reflector shifts cannot be approximated by them. For that matter, no near-diagonal approximation to the Hessian could be accurate in such cases, so these same limitations would seem to apply to all scaling methods.

Pseudodifferential operators are distinguished from other types of spatially varying filters by strong constraints on the filter coefficients 𝐴_𝑖𝑗, which must exhibit polynomial-like growth in the wavenumber 𝑘=√k⋅k, as 𝑘→∞. more precisely, require that a real 𝜇 exist so that the mixed partial derivative of order 𝛼=𝛼₁+𝛼₂+𝛼₃ grows like 𝑘^𝜇−𝛼 as 𝑘→∞: 𝜕^𝛼𝜕𝑘^𝛼1₁𝜕𝑘^𝛼2₂𝜕𝑘^𝛼3₃𝐴_𝑖𝑗(x,k)≤const.×𝑘^𝜇−𝛼. (5) That is, 𝐴_𝑖𝑗 grows like 𝑘^𝜇 for large 𝑘, and every mixed partial derivative in k decreases the order of growth by the order of the derivative. It is also required that spatial derivatives do not increase the order of growth of the coefficient—we will not write down precisely how this additional constraint is imposed, referring the interested reader to Taylor [10] for a more complete account.

Condition (5) also describes the growth behavior of polynomials in k, which are the filter coefficients of differential operators. Other functions of k, for instance 𝑘 itself, also satisfy these conditions, so the type of filter described by (5) is more general than differential operators. These considerations led to the name “pseudodifferential operator” being applied to this type of filter. In fact, it is possible to approximate any pseudodifferential operator arbitrarily well by the composition of a differential operator and a power of the Laplacian (i.e., a filter of the form 𝑘^𝜇). In this sense, pseudodifferential operators are nearly local, that is, diagonal.

By convention, filter coefficients obeying the growth rules described here are called symbols, and we shall use this terminology as well to distinguish these special filter coefficients from the Fourier representation of general filters. The number 𝜇 figuring in the constraint (5) is the order of the operator; it determines the extent to which application of the operator changes the growth (or decay) rate for in the Fourier domain, for large wavenumber. Note that order is inclusive, rather than precise: if an operator is of order 𝜇, it is of order 𝜈 for any 𝜈≥𝜇. However, order does differentiate pseudodifferential operators by the size of their effect on oscillatory data, just as is the case for differential operators. The work cited at the beginning of this section shows how the order of the Hessian depends on the dimension of the model and the source-receiver geometry. For example, for constant density acoustics in 2D, the order is 1; for 3D and a full range (i.e., 𝜋) of azimuths, the order is 2.

Equation (5) is a very strict constraint on the structure of the pseudodifferential class of spatially varying filters, so it should not be a surprise that the behavior of such filters is also highly constrained. Very important in the sequel are the composition properties: if 𝐴 is a (𝑛×𝑛 matrix) pseudodifferential operator of order 𝜇, 𝐵 is another of order 𝜈, then (1)the product (in other words, composition) 𝐴𝐵 is again a pseudodifferential operator of order 𝜇+𝜈, and the symbol of the product differs from the product of the symbols by a symbol of order 𝜇+𝜈−1;(2)consequently, if 𝑛=1 (scalar operators), then the commutator 𝐴𝐵−𝐵𝐴 is of order 𝜇+𝜈−1.

That is, to leading order in wavenumber, the product of pseudodifferential operators is simply obtained by multiplying their symbols, and scalar operators commute. These properties single out pseudodifferential operators amongst general spatially varying filters. Proofs of these and other facts about pseudodifferential operators may be found in Taylor [10], for example.

The theory cited above also showed that the normal operator is “partly invertible”, that is, at each point in the subsurface, the Hessian scales spatial Fourier components with a certain range of dips by a positive multiple of 𝑘^𝜇. It follows from point (2) above that an operator whose symbol is the reciprocal of that of the normal operator for this range of dips, and simply set to zero with suitable tapering outside this range will function as an approximate inverse Hessian. This construction gives an approximate inversion for the single-parameter inverse Hessian (e.g., for constant-density acoustics), illustrates the utility of these very special filters, and is explained briefly in Appendix B, and in more detail in [11, 12].

Multiparameter scattering results in a matrix Hessian, mixing influences between various parameters. In some cases, the Hessian is a matrix of pseudodifferential operators. For example, the Hessian for acoustic scattering or for polarized elastic scattering (𝑃-𝑃,𝑆-𝑆) has this property. Acoustic scattering has only one (𝑝-wave) phase, so it has this property [44], as does elastic scattering with suitably polarized data [43]. The close relation between the matrix 𝐴=(𝐴_𝑖𝑗) of symbols and the normal operator 𝑁 which it defines again allows us to calculate an approximate inverse via symbol computations.

To this end, recall the definition of the adjugate, or classical adjoint, of an 𝑛×𝑛 matrix A, denoted by Adj(𝐴): it is the transpose of the matrix of cofactors of 𝐴. Explicitly, Adj(𝐴)_𝑖𝑗=(−1)^𝑖+𝑗det𝑀_𝑗𝑖 (6) in which the (𝑛−1)×(𝑛−1) minor 𝑀_𝑖𝑗 is obtained from 𝐴 by dropping the 𝑖th row and 𝑗th column. The 𝑗𝑖 cofactor is the number on the right-hand side of (6). In our case, the transpose of the cofactor matrix may be ignored, as the matrix 𝐴 and thus its adjugate are symmetric (since the Hessian 𝑁 is symmetric).

The significance of the adjugate is this: when the matrix 𝐴 is invertible, then the adjugate is proportional to the inverse, 𝐴⁻¹1∶=(𝐴)(𝐴)detAdj. (7) This relation is Cramer's rule and is the critical observation leading to our approximation method. Note that Adj(𝐴)𝐴=𝐴Adj(𝐴)=det(𝐴)𝐼, (8) where 𝐼 is the identity matrix. If we define the adjugate of 𝑁 to be Adj(𝑁)=op(Adj(𝐴)) (i.e., the matrix pseudodifferential operator with symbol matrix Adj(𝐴)), with slight abuse of notation, then the multiplicative order properties of pseudodifferential operators described above imply that Adj(𝑁)𝑁≈𝑁Adj(𝑁)≈det(𝑁)𝐼. (9) The equation above features another abuse of notation, with det(𝑁)∶=op(det(𝐴)). The power of (8) is revealed when applied to (3) Adj(𝑁)𝑚_mig=Adj(𝑁)𝑁𝑚≈det(𝑁)𝑚. (10) Equation (10) reduces the inversion of 𝑁 to two steps: (i) application of Adj(𝑁), followed by (ii) inversion of the scalar operator det(𝑁), after the application of the adjugate.

Note also that while 𝑚_mig may be dimensionally heterogeneous, Adj(𝑁)𝑚_mig has redimensionalized all components so that their units all bear the same relation to the units of the components of 𝑚. Inversion of det(𝑁) removes this common unit factor and correctly dimensionalizes the model estimate.

For inversion of det(𝑁), we resort to a method similar to the one we previously developed for 𝑛=1 (reviewed in Appendix B). First apply the normal operator again to the left-hand side of (10), to form 𝑁Adj(𝑁)𝑚_mig(𝑁)(𝑁)(𝑁)𝑚=𝑁det𝑚≈det𝑁𝑚=det_mig. (11) Here, we have used the fact that scalar pseudodifferential operators approximately commute with matrices of pseudodifferential operators to commute 𝑁 and det(𝑁).

Equation (11) shows how to find a pseudodifferential approximate inverse 𝑐 to det𝑁, by minimizing an objective function, all parts of which are readily computable 𝑚_mig−𝑐𝑁Adj(𝑁)𝑚_mig², (12) over a suitable class of pseudodifferential operators 𝑐. From (9), a solution 𝑐 will approximately invert det(𝑁) and is scalar, since det(𝑁) is. Since 𝑐 removes from the left-hand side of (10) a scalar filter, it plays a similar role to the scaling factors encountered in the various scaling methods mentioned in the introduction, so we will also refer to 𝑐 as a scaling factor, although it is a pseudodifferential operator.

The reader will recall that 𝑁 is only partially invertible for a range of model wavevectors k, the symbol determinant det(𝐴) is asymptotic to a positive multiple of 𝑘^𝑛𝜇 and to a lower power outside of this “illuminated” cone of model wavenumbers. The method explained in Appendix B constructs 𝑐 as an approximation to the reciprocal of the symbol of det(𝑁) in the illuminated wavenumber cone, tapered to zero outside of this cone.

The algorithm explained in Appendix B is iterative and requires at least one application of the trial scaling factor 𝑐 at each iteration. An efficient method for numerical application of pseudodifferential operators is therefore required; we use an algorithm introduced by Bao and Symes [17], which we review in Appendix A.

Having constructed 𝑐, we approximate the solution of normal equations by 𝑚_inv=𝑐Adj(𝑁)𝑚_mig. (13)

Equations (12) and (13) show that computation of the adjugate action is key to this approach, so we examine it more carefully. We will address the two parameter cases (𝑛=2) explicitly, both for inspiration and because it is this case for which we present examples in the next section. Writing 𝑁𝑁=¹¹𝑁₁₂𝑁₁₂𝑁₂₂, (14)Adj(𝑁) is (𝑁)=𝑁Adj²²−𝑁₁₂−𝑁₁₂𝑁₁₁,(𝑁)𝑁Adj𝑚=²²𝑚₁−𝑁₁₂𝑚₂−𝑁₁₂𝑚₁+𝑁₁₁𝑚₂. (15) It appears to be necessary to compute the various components of 𝑁 separately: for example, if the underlying theory is acoustic scattering, then 𝑁₂₂ might represent input and output density perturbations, and so on, so that each 𝑁_𝑖𝑗 would appear to require a complete modeling/migration cycle.

In fact, it is possible to compute Adj(𝑁)𝑚 using only applications of 𝑁, thus both reducing the cost of the application and eliminating the need to write special component-to-component migration codes. In the 2×2 case, only one application of 𝑁 is required! In fact, Adj(𝑁)=𝐽^𝑇𝑁𝐽, (16) here, 𝐽 is the so-called symplectic matrix, 𝑚𝐽=0−𝐼𝐼0,so𝐽¹𝑚₂=−𝑚₂𝑚₁. (17) Note that units are implicitly changed (as they are in the application of the adjugate operator) and that in practice it is presumed that the discrete representations of the two parameters are the same, so that, for example, density can be swapped for velocity and vice versa.

Equation (16) implies that the application of the adjugate on the migrated image requires one application of 𝑁, rather than the four that a naive implementation would suggest. Thus, approximating the inverse for 𝑛=2 requires two applications of the normal operator, plus a small amount of additional data manipulation.

The situation is more complicated for 𝑛=3, but a similar result holds: it is possible to compute the application of the adjugate using five applications of the normal operator, a considerable improvement over the twenty-four that the naive algorithm based on the definition (6) would appear to require; note that separate application of every component 𝑁_𝑖𝑗 by wave-equation methods is just as expensive as application of the entire operator 𝑁. See Nammour [38] for details. We conjecture that a similar drastic reduction is possible for general 𝑛.

As a final note, we can show explicitly for 𝑛=2 the sense in which equation (9) is an approximation, while equation (8) is exact. Denoting the commutator of two operators by [𝑎,𝑏]=𝑎𝑏−𝑏𝑎, Adj(𝑁)𝑁=𝑁₂₂𝑁₁₁−𝑁²₁₂+0𝑁1001₂₂,𝑁₁₂𝑁₁₁,𝑁₁₂𝑁₁₁,𝑁₂₂(𝑁)0𝑁=det𝐼+²²,𝑁₁₂𝑁₁₁,𝑁₁₂𝑁₁₁,𝑁₂₂. (18) The theory of pseudodifferential operators shows that the commutator of two pseudodifferential operators is of lower order than their composition, that is, relatively negligible for highly oscillatory input fields. Using ≈ to mean “differs by a lower order error”, Adj(𝑁)𝑁≈𝑁Adj(𝑁)≈det(𝑁). (19)

3. Numerical Examples: Layered Variable Density Acoustics

As a first application to the two parameter inversion, we construct a variable density acoustics model perturbation consisting of a thin oscillatory velocity layer and a thin oscillatory density layer in a different place (see Figures 1 and 2). The background model is homogeneous with 𝑣𝑝=2km/s and 𝑑𝑛=2000kg/m³. The two thin layers are treated as perturbations.

Figure 1 

vp, velocity perturbation.

Figure 2 

dn, density perturbation.

We simulated reflection data for this model using the IWAVE software developed by The Rice Inversion project in linearized (Born) modeling mode [45]. The model extends around 1.7 km in depth and 6.5 km horizontally. One source is put in the middle, and receivers are laid out to created an offset ranging [−2.7 km, 2.7 km], spaced 20 m apart and 40 m below the (absorbing) surface. The isotropic point source wavelet was a 2.5–5–15–20 Hz trapezoidal zero-phase bandpass filter. 3 seconds of data were recorded at 151 receiver positions. The boundary conditions were absorbing on all sides of the domain, in particular, free surface effects were not modeled. Since the model is layered, only a single source gather need be modeled.

All (reverse-time) migrations were also carried out with IWAVE [46]. The output is a layered model (i.e., two depth traces, one for each velocity and density perturbations), which we shall display as layered fields.

Migrating the Born data shows how migration mixes the effects of the two models in the two components of the migrated images (Figure 3). We shall refer to the migrated images as 𝑚_mig1 and 𝑚_mig2 to remain consistent with our notation where the vector of migrated images is 𝑚_mig. This example, albeit simple, stresses a new challenge of multiparameter inversion. For one parameter inversion, the events in the migrated image corresponded to events in the true model. In multiparameter inversion, events in the migrated images may correspond to an event in one or more of the components of the model. It is virtually impossible to tell that these migrated images correspond to a model with separate events for velocity and density without successful inversion.

(a) 𝑚mig1, velocity component of the migrated image

(b) 𝑚mig2, density component of the migrated image

Figure 3 

Migrated images mixing the contributions from density and velocity, and effecting a phase space scaling.

Applying the scheme outlined above, we form 𝐽^𝑇𝑁𝐽𝑚_mig≈det(𝑁)𝑚. (20) The result is shown in Figure 4 and shows how one application of the normal operator effectively separated the contributions of the velocity and density events. It remains to effect an amplitude correction by approximating an inverse to det(𝑁). For this end, we are required to form 𝑁det(𝑁)𝑚≈det(𝑁)𝑚_mig, shown in Figure 5.

(a) (𝐽𝑇𝑁𝐽𝑚mig)1≈det(𝑁)𝑚1

(b) (𝐽𝑇𝑁𝐽𝑚mig)2≈det(𝑁)𝑚2

Figure 4 

The application of the adjugate separates the velocity and density contributions. This result is a phase space scaling of the true model.

(a) det(𝑁)𝑚mig1

(b) det(𝑁)𝑚mig2

Figure 5 

Scaling of the migrated images by det(𝑁), used to undo the determinant.

The final step corrects the amplitudes of det(𝑁)𝑚 by undoing the effect of det(𝑁), which yields an approximate inverse. This final step complements the separation we obtained earlier with an amplitude correction Figure 6 shows that the approximate inverse compares favorably with true model. An interesting observation on this result is the fact that the velocity model is better recovered than the density model: traces of the velocity event in the density model are more apparent than those of the density event in the velocity model. This observation is in accordance with the theoretical fact that the recovery of velocity in variable density acoustics is better conditioned than the recovery of density. The inversion is successful in that the inverted model succeeds in fitting 70% of the target data (see Figure 7).

(a) invvp

(b) invdn

Figure 6 

The approximate inverse. The contributions from velocity and density are separated, and the amplitudes are corrected.

(a) Data misfit

(b) Target Data

Figure 7 

Data misfit versus target data. The inverted model fits 70% of the data.

4. Discussion and Conclusion

We have presented a method to approximate the solution of the multiparameter linearized inverse problem an extension of Cramer's rule for matrices of pseudodifferential operators. The method consists of two steps: (1)reducing the multiparameter problem to a one parameter problem, which yields an amplitude scaling of the solution; (2)correcting the amplitudes of the result from the previous step to approximate the solution.

The application of the work flow to a layered variable density acoustics example shows how the effects of different material parameters, indistinguishable in the migrated images, are separated in the first step. The amplitude correction step successfully yields an approximate solution to the linearized inverse problem.

The work flow presented above applies without modification to any model. Results of application to models more complex than the layered model described above will be presented elsewhere. We point out that other aspects of the physics of seismic data generation will need to be accommodated if this (or any other) algorithm is to extract accurate results from field data. For example, as shown already by Minkoff and Symes [39], even when the Born approximation is adequate, the radiation pattern of the source has a first-order effect on the variation of amplitude with angle and must therefore be included in the parameters to be estimated.

We have extended the Pseudodifferential Operator Algorithm (Appendix A) to three spatial dimensions (3D) which allows approximate inversion of 3D models. The symbol is represented by a truncated spherical harmonics expansion. The formulation of the method in this case remains intact, though the cost of modeling and migration increases. This extension of the method is described in detail, and examples are given in Nammour [38].

The explicit derivation analogous to 𝑛=2 for more parameters (namely 𝑛=3) is similar to the one presented here but requires tedious algebraic manipulations. The case 𝑛=3 becomes relevant in the linearized inverse problem for linear elasticity, for example. Six applications of the normal operator are required to approximate the solution of the linearized inverse problem. See Nammour [38] for a detailed discussion.

We end by reminding the reader that the approach to approximate inversion pursued here relies on the validity of the pseudodifferential characterization of the Hessian. As mentioned in the theory section, two of the three conditions (background model smooth on the wavelength scale, polarized data) are essential for any near-diagonal Hessian approximation, adaptive or not, to be effective. These limitations seem likely to be acceptable in some exploration settings, for example, inversion of marine data for supersalt sediment properties with soft sea floor, but may seriously obstruct this type of technique in other settings.

Appendices

A. The Pseudodifferential Operator Algorithm

Bao and Symes [17] presents an algorithm to numerically apply pseudodifferential operators, which we use to represent the scaling factor 𝑐 in (12). We will discuss this algorithm in this section as it is central to the feasibility of the method presented here.

This discussion is restricted to 2D, so we may write 𝐱=(𝑥,𝑧). Recall that a pseudodifferential operator is characterized by its symbol and defined by 𝑄_𝑚𝑢(𝑥,𝑧)=𝑞_𝑚𝑥,𝑧,𝑘_𝑥,𝑘_𝑧𝑢𝑘_𝑥,𝑘_𝑧𝑒^{𝑖(𝑥𝑘𝑥+𝑧𝑘𝑧)}𝑑𝑘_𝑥𝑑𝑘_𝑧, (A.1) where 𝑞_𝑚(𝑥,𝑧,𝑘_𝑥,𝑘_𝑧) is the principal symbol, homogeneous of degree 𝑚, and 𝑢=ℱ[𝑢] is the Fourier transform of 𝑢.

Thus, writing 𝑘_𝑥=𝜔cos𝜃, 𝑘_𝑧=𝜔sin𝜃 and using the homogeneity of 𝑞_𝑚, we have 𝑞_𝑚𝑥,𝑧,𝑘_𝑥,𝑘_𝑧=𝜔^𝑚𝑞_𝑚(𝑥,𝑧,𝜃). (A.2)

Notice that 𝑞_𝑚(𝑥,𝑧,𝜃)=𝑞_𝑚(𝑥,𝑧,cos𝜃,sin𝜃) is periodic and smooth in 𝜃, and hence it admits a rapidly converging Fourier expansion. We thus truncate the Fourier series, approximating the symbol by its first 𝐾+1 Fourier modes, 𝑞_𝑚(𝑥,𝑧,𝜃)≈^𝑙=𝐾/2_{𝑙=−𝐾/2}𝑐_𝑙(𝑥,𝑧)𝑒^𝑖𝑙𝜃=^𝑙=𝐾/2_{𝑙=−𝐾/2}𝜔^−𝑙𝑐_𝑙(𝑥,𝑧)𝑘_𝑥+𝑖𝑘_𝑧^𝑙. (A.3) Plugging (A.3) into (A.1), we obtain 𝑄_𝑚𝑢(𝑥,𝑧)≈^𝑙=𝐾/2_{𝑙=−𝐾/2}𝑐_𝑙(𝑥,𝑧)ℱ⁻¹𝜔^𝑚−𝑙𝑘_𝑥+𝑖𝑘_𝑧^𝑙𝑢𝑘_𝑥,𝑘_𝑧. (A.4)

Fourier transform theory identifies 𝜔^𝑚−𝑙 as the symbol of (−∇)^{(𝑚−𝑙)/2}, and 𝑘_𝑥 and 𝑘_𝑧 are, respectively, the symbols of 𝐷_𝑥=−𝑖𝜕_𝑥 and 𝐷_𝑧=−i𝜕_𝑧.

Sampling the field 𝑢(𝑥,𝑧) and the symbol 𝑞_𝑚(𝑥,𝑧,𝜃), 𝑈_𝑖𝑗=𝑢𝑥₀+(𝑖−1)Δ𝑥,𝑧₀+,𝑄𝑗−1Δ𝑧_𝑖𝑗𝑘=𝑞_𝑚𝑥₀+(𝑖−1)Δ𝑥,𝑧₀+,𝐾𝑗−1Δ𝑧,𝑘Δ𝜃𝑖=1,…,𝑀,𝑗=1,…,𝑁,𝑘=−2𝐾,…,2. (A.5)

Choosing Δ𝑘_𝑥=1/(𝑀−1)Δ𝑥 and Δ𝑘_𝑧=1/(𝑁−1)Δ𝑧 yields the unaliased discretizations of the symbols of the square root of the negative Laplacian, 𝐷_𝑥 and 𝐷_𝑧Ω_pr=2𝜋𝑝Δ𝑘_𝑥²+𝑟Δ𝑘_𝑧²,𝑋_pr=2𝜋𝑝Δ𝑘_𝑥,𝑍_pr=2𝜋𝑟Δ𝑘_𝑧,𝑀𝑝=−2𝑀,…,2𝑁,𝑟=−2𝑁,…,2. (A.6)

Equation (A.4) suggests the following algorithm to estimate 𝑄_𝑚𝑢 [17]. All Fourier transforms refer to a discrete Fourier transform. (1)Compute 𝑈_pr=ℱ[𝑈_𝑖𝑗]. (2)For each 𝑖∈[1,𝑀] and 𝑗∈[1,𝑁], compute 𝑄_𝑖𝑗={𝑄_𝑖𝑗𝑙}^𝐾/2_{𝑙=−𝐾/2} the discrete Fourier transform of 𝑄_𝑖𝑗={𝑄_i𝑗𝑘}^𝐾/2_{𝑘=−𝐾/2}.(3)Initialize (𝑄𝑈)_𝑖𝑗=0, for 𝑖∈[1,𝑀], 𝑗∈[1,𝑁], For 𝑙=−𝐾/2∶𝐾/2(a)compute {𝑅^𝑙_𝑖𝑗}^𝑀,𝑁_{𝑖=1,𝑗=1}=ℱ⁻¹[Ω^𝑚−𝑙_pr(𝑋_pr+𝑖𝑍_pr)^𝑙𝑈_pr] for 𝑝=−𝑀/2,…,𝑀/2 and 𝑟=−𝑁/2,…,𝑁/2,(b)accumulate(𝑄𝑈)_𝑖𝑗=(𝑄𝑈)_𝑖𝑗+𝑄_𝑖𝑗𝑙𝑅^𝑙_𝑖𝑗 End

A straightforward discretization of (A.1) has a computational complexity of 𝑂(𝑁⁴log(𝑁)). The algorithm described above uses FFT (Fast Fourier Transform) and thus exhibits a complexity of 𝑂(𝐾𝑁²(log(𝑁)+log(𝐾))). The appeal of this approach is that 𝐾 is independent of 𝑁. In fact, applications to reflection seismology require that the symbol be smooth and slowly varying in 𝜃, thus, may be captured accurately by a modest number of Fourier modes or, more explicitly, a small 𝐾.

The dependence on dip is captured in the angle variable 𝜃, and the method allows us to capture multiple dip events by increasing 𝐾>1.

B. The One-Parameter Case

In this appendix, we review the method developed in [11, 12] to construct an approximate inverse Hessian in the one-parameter case. The aim is to solve 𝑁𝑚=𝑚_mig, (B.1) where 𝑚_mig=𝐹^∗𝑑∈Range(𝑁). We seek a pseudo differential scaling factor and formulate its recovery as an optimization problem. Given the migrated image m_mig and the remigrated image 𝑁𝑚_mig, 𝑐=argmin_𝑐∈ΨDO𝑚_mig−𝑐𝑁𝑚_mig². (B.2) The scaling factor 𝑐 is chosen from a class of pseudodifferential operators represented numerically using the PsiDO algorithm (Appendix A). In this setting, the scaling factor approximates the action of the inverse of the normal operator on the migrated image 𝑚_mig. More specifically, 𝑚=𝑁⁻¹𝑚_mig≈𝑁⁻¹𝑐𝑁𝑚_mig≈𝑐𝑁⁻¹𝑁𝑚_mig=𝑐𝑚_mig∶=𝑚_inv. (B.3) The first of these equations expresses the true solution 𝑚, the second approximate equality follows from because pseudodifferential operators approximately commute. Defining 𝑚_inv∶=𝑐𝑚_mig thus yields an approximation to the true model 𝑚. Equation (B.3) shows that the scaling factor approximates the action of the inverse on the normal operator on the migrated image, and it is only in that sense that 𝑐 approximates 𝑁⁻¹.

The scaling factor 𝑐 is represented explicitly by the PsiDO algorithm, 𝑐𝑁𝑚_mig()≈𝑥,𝑧_{𝑙=𝐾/2𝑙=−𝐾/2}𝑐_𝑙(𝑥,𝑧)ℱ⁻¹𝜔^𝑚−𝑙𝑘_𝑥+𝑖𝑘_𝑧^𝑙𝑁𝑚_mig𝑘_𝑥,𝑘_𝑧. (B.4)

We enforce the continuity of 𝑐_𝑙(𝑥,𝑧) using a parsimonious basis technique. Let {𝜓_𝑗(𝑥,𝑧)}^𝐽_𝑗=1be a set of smooth shape functions (cubic b-splines, for example). Write 𝑐_𝑙(𝑥,𝑧)=^𝐽_𝑗=1𝑐^𝑗_𝑙𝜓_𝑗(𝑥,𝑧), (B.5) and plug into the objective function, 𝐽=𝑚_mig−𝑐𝑁𝑚_mig²=_{𝑥𝑖,𝑧𝑘}𝑚_mig𝑥_𝑖,𝑧_𝑘−𝑐𝑁𝑚_mig𝑥_𝑖,𝑧_𝑘². (B.6) We have denoted by 𝐽 the objective function in (B.2), and used the fact that the images 𝑚_mig and 𝑁𝑚_mig are discretized on a grid {𝑥_𝑖}^𝑁_𝑖=1, {𝑧_𝑘}^𝑀_𝑘=1 to write out the norm explicitly.

We use limited memory BFGS (lBFGS) to minimize the objective function (B.2), which requires the gradient to be supplied by the user. The gradient is easily derived as 𝜕𝐽𝜕𝑐^𝑙_𝑗=2_{𝑥𝑖,𝑧𝑘}𝑚_mig𝑥_𝑖,𝑧_𝑘−𝑐𝑁𝑚_mig𝑥_𝑖,𝑧_𝑘𝜓_𝑗𝑥_𝑖,𝑧_𝑘×ℱ⁻¹𝜔^𝑚−𝑙𝑘_𝑥+𝑖𝑘_𝑧^𝑙𝑁𝑚_mig𝑘_𝑥,𝑘_𝑧×𝑥_𝑖,𝑧_𝑘. (B.7) We have also explored enforcing the positivity of the scaling factor explicitly and other symmetries of the problem, which yields to slightly different representations of the scaling factor 𝑐 and in turn different objective function gradient. For more details on this point and a discussion about the ability of the method to resolve multiple dip events please consult [12].

Acknowledgments

This work was supported in part by the National Science Foundation under grant DMS 0620821, and by the sponsors of The Rice Inversion Project. The authors are grateful to an anonymous reviewer and to Guest Editor Sergey Fomel for thoughtful critiques.