Abstract

Direct inverse methods solve the problem of interest; in addition, they communicate whether the problem of interest is the problem that we (the seismic industry) need to be interested in. When a direct solution does not result in an improved drill success rate, we know that the problem we have chosen to solve is not the right problem — because the solution is direct and cannot be the issue. On the other hand, with an indirect method, if the result is not an improved drill success rate, then the issue can be either the chosen problem, or the particular choice within the plethora of indirect solution methods, or both. The inverse scattering series (ISS) is the only direct inversion method for a multidimensional subsurface. Solving a forward problem in an inverse sense is not equivalent to a direct inverse solution. All current methods for parameter estimation, e.g., amplitude-variation-with-offset and full-waveform inversion, are solving a forward problem in an inverse sense and are indirect inversion methods. The direct ISS method for determining earth material properties defines the precise data required and the algorithms that directly output earth mechanical properties. For an elastic model of the subsurface, the required data are a matrix of multicomponent data, and a complete set of shot records, with only primaries. With indirect methods, any data can be matched: one trace, one or several shot records, one component, multicomponent, with primaries only or primaries and multiples. Added to that are the innumerable choices of cost functions, generalized inverses, and local and global search engines. Direct and indirect parameter inversion are compared. The direct ISS method has more rapid convergence and a broader region of convergence. The difference in effectiveness increases as subsurface circumstances become more realistic and complex, in particular with band-limited noisy data.

Introduction

Seismic processing is an inverse problem to determine the properties of a medium from measurements of a wavefield exterior to the medium. The ultimate inversion objective of seismic processing in seismic exploration is to use recorded reflection data to extract useful subsurface information that is relevant to the location and production of hydrocarbons. There is typically a coupled chain of intermediate steps and processing that takes place toward that objective, and I refer to each of those intermediate steps, stages, and tasks as objectives “associated with inversion” or inverse tasks toward the ultimate subsurface information extraction goal and objective. All seismic processing methods that are used to extract subsurface information make assumptions and have prerequisites.

A seismic method will be effective when those assumptions/conditions/requirements are satisfied. When those assumptions are not satisfied, the method can have difficulty and/or will fail. That failure can and will contribute to processing and interpretation difficulties with subsequent dry-hole exploration well drilling or drilling suboptimal appraisal and development wells.

Challenges in seismic processing and seismic exploration and production are derived from the violation of assumptions/requirements behind seismic processing methods. Advances in seismic processing effectiveness are measured in terms of whether the new capability results in/contributes to more successful plays and better informed decisions and an increased rate of successful drilling.

The purpose of seismic research is to identify and address seismic challenges and to thereby add more effective options to the seismic processing toolbox. These new options can be called upon when indicated, appropriate, and necessary as circumstances dictate.

No toolbox option is the appropriate choice under all circumstances. For example, the most effective method, from a technical perspective, might be more than is necessary and needed, under a given circumstance, and a less effective and often less costly option could be the appropriate and indicated choice. Under other more complex and daunting circumstances, the more effective and (perhaps) more costly option will be the only possible choice that is able to achieve the objective of that processing task and interpretation goal. The objective is to expand the number of options in the seismic toolbox to allow a capable response to a larger number of circumstances. As I will point out below, “identify the problem” is the first, the essential, and sometimes the most difficult (and often the most ignored and/or underappreciated) aspect of seismic research.

Identifying and delineating the violation of assumptions behind seismic processing methods is an absolutely essential first step in a strategy and plan for developing a response to prioritizing and pressing seismic exploration challenges. This paper provides a new insight, and advance for the first and critical step of addressing seismic processing challenges: problem identification.

I explain in detail and exemplify why only a direct inversion method can help us to decide whether the problem we (the seismic industry) are interested in addressing is, in fact, the problem we need to address.

Seismic processing methods can be classified as based on either statistical models and principles or wave-theory concepts and approaches. Wave-theory concepts used in seismic processing can be further catalogued as modeling and inversion.

In the next section, I describe these two wave theory approaches to seismic processing, that is, modeling and inversion, and I will further distinguish between direct and indirect inversion methods. That clarification represents a central theme and objective of this paper.

Modeling and inversion

Modeling, as a seismic processing tool, starts with a prescribed wavefield source mechanism and a model type (e.g., acoustic, elastic, anisotropic, or anelastic), and then properties are defined within the model type for a given medium (e.g., velocities, density, and attenuation Q). The modeling procedure then provides the seismic wavefield that the energy source produces at all points inside and outside the medium.

Inversion also starts with an assumed known and prescribed energy source outside the medium. In addition, the wavefield outside the medium is assumed to be recorded and known. The objective of seismic inversion is to use the latter source description and wavefield measurement information to make inferences about the subsurface medium that are relevant to the location and production of hydrocarbons.

Direct and indirect inversion

Inversion methods can be classified as direct or indirect. A direct inversion method solves an inverse problem (as its name suggests) directly. On the other hand, an indirect inversion method seeks to solve an inverse problem circuitously through indirect approaches that often call up assumed aligned objectives or conditions. There are times when the indirect approach will seek to satisfy necessary (but typically not sufficient) conditions, and properties, and it is often mistakenly considered and treated as though it was equivalent to a direct method and solution. Indirect methods come in many varieties; some are obvious, and others are more subtle and harder to identify as being indirect. Among indicators, identifiers, and examples of “indirect” inverse solutions (Weglein, 2015a) are (1) model matching, (2) objective/cost functions, (3) local and global-search algorithms, (4) iterative linear inversion, (5) methods corresponding to necessary but not sufficient conditions, e.g., common-image gather flatness as an indirect migration velocity analysis method, and (6) solving a forward problem in an inverse sense, e.g., amplitude-variation-with-offset (AVO) and full-waveform inversion (FWI). Regarding the last indirect indicator, item (6), I will show that solving a forward problem in an inverse sense is not equivalent to a direct inverse solution for those same objectives.

As a simple illustration, a quadratic equation  
ax2+bx+c=0
(1)
can be solved through a direct method as  
x=b±b24ac2a,
(2)
or it can be solved by an indirect method searching for x, such that, e.g., some functional of  
(ax2+bx+c)2
(3)
is a minimum.

In the next section, this example will be further discussed and examined as a way to introduce and develop fundamental concepts in a simple and transparent context. The lesson gleaned from that simple example will later (in this paper) be extended and applied to the more complicated and relevant seismic inverse formulations and methods. In Weglein (2013), there is an introduction to the subject of direct and indirect inverse solutions, which provides a useful background reference for this paper and contains several indirect inversion references (Blum, 1972; Keys and Weglein, 1983; Gauthier et al., 1986; Tarantola, 1986, 1987; Crase et al., 1990; Symes and Carazzone, 1991; Chavent and Jacewitz, 1995; Matson, 1997; Nolan and Symes, 1997; Weglein and Matson, 1998; Biondi and Sava, 1999; Brandsberg-Dahl et al., 1999; Pratt, 1999, Pratt and Shipp, 1999; Rickett and Sava, 2002; Weglein et al., 2002, 2011; Sava and Fomel, 2003; Biondi and Symes, 2004; Sava et al., 2005; Valenciano et al., 2006; Iledare and Kaiser, 2007; Ben-Hadj-ali et al., 2008; Symes, 2008; Vigh and Starr, 2008; Baumstein et al., 2009; Ben-Hadj-ali et al., 2009; Brossier et al., 2009; Hawthorn, 2009; Sirgue et al., 2009, 2010, 2012; Liang et al., 2010; Ferreira, 2011; Fichtner, 2011; Li et al., 2011, Luo et al., 2011; Anderson et al., 2012; Guasch et al., 2012, Kapoor et al., 2012; Weglein, 2012a, 2012b; Zhou et al., 2012; Zhang and Biondi, 2013).

The important quadratic equation example

The direct quadratic formula solution equation 2 explicitly and directly outputs the exact roots (for all values of a, b, and c) when the roots are real and distinct, a real double root, and imaginary and complex roots. The quadratic equation and quadratic solution provide a very simple and insightful example. How would a search algorithm know after a double root is found that it is the only root and not to keep looking and searching forever for a second, nonexistent, root? How would a search algorithm know to search for only real or for real and complex roots? How would a search algorithm accurately locate an irrational root such as 31.732, as x=(b±b24ac)/2a would directly and precisely and immediately produce? Indirect methods such as model matching and seeking and searching and determining roots as in equation 3 are ad hoc and do not derive from a firm framework and foundation and never provide the confidence that we (the seismic industry) are actually solving the problem of interest.

What is the point in discussing the quadratic formula? And what is the practical big deal about a direct solution?

How can this example and discussion of the quadratic equation possibly be relevant to exploration seismology? Please imagine for a moment that equation 1ax2+bx+c=0 was an equation whose inverse and solution for x given by equation 2x=(b±b24ac)/2a had seismic exploration drill location prediction consequence. And furthermore, suppose that this direct solution for x did not lead to successful and/or improved drilling decisions. Under the latter circumstance, we (the seismic industry) could not blame or question the method of solution of equation 1 because equation 2 is direct and unquestionably solving equation 1. If equation 2 was not producing useful and beneficial results, we know that our starting equation 1 is the issue, and we have identified the problem. The problem we thought we need to solve (equation 1) is not the problem we need to solve. In contrast with equation 3, an indirect method, any lack of drilling prediction improvement and added value or other negative exploration consequences could be due to either the equation you are seeking to invert and/or the boundless, unlimited selection, and the plethora of indirect methods using either partial or full recorded wavefields.

That lack of clarity and definitiveness within indirect methods obfuscates the underlying issue and makes identification of the problem (and what is behind a seismic challenge) considerably more difficult to identify and to define. Indirect methods with search engines, such as equation 3, lead to “workshops” for solving equation 1 and grasping at mega high-performance computing (HPC) straws (and capital expenditure investment for buildings full of HPC) that are required to search, seek, and locate “solutions.” The more HPC we invest in, and is required, the more we are literally “buying-in,” and as stake-holders, we become committed and therefore convinced of the unquestioned validity of the starting point and our indirect thinking and methodology.

Therefore, beyond the benefit of a direct method, such as equation 2 providing assurance that we are actually solving the problem of interest (equation 1), there is the unique problem location and identification benefit of a direct inverse when a seismic analysis, processing, and interpretation produces unsatisfactory E&P results.

To bring this (quadratic equation example) closer to the seismic experience, please imagine hypothetically that we are not satisfied (in terms of improved drill location and success rate) with a direct inverse of the elastic-isotropic equation for amplitude analysis. Because we were using a direct inversion solution, we know we need to go to a different starting point, perhaps with a more complete and realistic model of wave propagation because we can exclude the direct inverse solution method as the problem and issue. That is an example of determining that a problem of interest is not the same problem we need to be interested in.

How to distinguish between the “problem of interest” and the problem we need to be interested in

Direct inverse methods provide value for knowing that you have actually solved the problem of interest. Furthermore, with direct inverse solutions, there is the enormous additional value of determining whether our starting point, the problem of interest, is in fact the problem we need to be interested in.

Scattering theory and the forward and inverse scattering series: The basis of direct inversion theory and algorithms

Scattering theory is a form of perturbation theory. It provides a direct inversion method for all seismic processing objectives realized by a distinct isolated task subseries of the inverse scattering series (ISS) (Weglein et al., 2003). Each term in the ISS (and the distinct and specific collection of terms that achieve different specific inversion associated tasks) is computable (1) directly and (2) in terms of recorded reflection data and without any subsurface information known, estimated, or determined before, during, or after the task is performed and the specific processing objective is achieved.

For certain distinct tasks, and subseries, e.g., free-surface multiple elimination and internal multiple attenuation, the algorithms not only do not require subsurface information but in addition possess the absolutely remarkable property of being independent of the earth model type (Weglein et al., 2003). That is, the distinct ISS free-surface and internal multiple algorithms are unchanged, without a single line of code having the slightest change for acoustic, elastic, anisotropic, and anelastic earth models (Weglein et al., 2003; Wu and Weglein, 2014). For those who subscribe to indirect inversion methods as, e.g., the “be-all and end-all” of inversion with various model matching approaches, it would be a useful exercise for them to consider how they would formulate a model-type-independent model-matching scheme for free-surface and internal multiple removal. It is not conceivable, let alone realizable, to have a model-type-independent model matching scheme.

For the specific topic and focus of this paper, the inversion task of parameter estimation, there is an obvious need to specify the model type and what parameters are to be determined. Hence, it is for that parameter estimation/medium property objective, and that model-type-specific ISS subseries, that the difference between the problem of interest and the problem that we need to be interested in, is relevant, central, and significant. Only direct inversion methods for earth mechanical properties provide that assumed earth model-type evaluation, clarity, and distinction.

The basic operator identity that relates a change in a medium and the change in the wavefield

A direct inverse solution for parameter estimation can be derived from an operator identity that relates the change in a medium’s properties and the commensurate change in the wavefield. That operator identity is general and can accommodate any seismic model type, for example, acoustic, elastic, anisotropic, heterogeneous, and anelastic earth models. That operator identity can be the starting point and basis of (1) perturbative scattering-theory modeling methods and (2) a firm and solid math-physics foundation and framework for direct inverse methods.

Theory

Let us consider an energy source that generates a wave in a medium with prescribed properties. With the same energy source, let us consider a change in the medium and the resulting change in the wavefield inside and outside the medium. Scattering theory is a form of perturbation theory that relates a change (or perturbation) in a medium to a corresponding change (or perturbation) in the original wavefield. When the medium changes, the resulting wavefield changes. The direct inverse solution (Weglein et al., 2003; Zhang, 2006) for determining earth mechanical properties is derived from the operator identity that relates the change in a medium’s properties and the commensurate change in the wavefield within and exterior to the medium. Let L0, L, G0, and G be the differential operators and Green’s functions for the reference and actual media, respectively, that satisfy  
L0G0=δandLG=δ,
(4)
where δ is a Dirac delta function. I define the perturbation operator V and the scattered wavefield ψs as follows:  
VL0LandψsGG0.
(5)

The operator identity

The relationship (called the Lippmann-Schwinger or scattering theory equation)  
G=G0+G0VG
(6)
is an operator identity that follows from  
L1=L01+L01(L0L)L1,
(7)
and the definitions of L0, L, and V.

Direct forward series and direct inverse series

The operator identity equation 6 (for a fixed-source function) is the exact relationship between changes in a medium and changes in the wavefield; it is a relationship between those quantities and not a solution. However, the operator identity equation 6 can be solved for G as  
G=(1G0V)1G0,
(8)
and expanded as  
G=G0+G0VG0+G0VG0VG0+.
(9)
The forward modeling of the wavefield G from equation 9 for a medium described by L is given in terms of the two parts of L, that is, L0 and V. The differential operator L0 enters through G0, and V enters as V itself. Equation 9 communicates that modeling using scattering theory requires a complete and detailed knowledge of the earth model type and medium properties within the model type. Equation 9 communicates that any change in medium properties, V, will lead to a change in the wavefield, GG0 that is always nonlinearly related to the medium property change, V. Equation 9 is called the Born or Neumann series in the scattering theory literature (see, e.g., Taylor, 1972). Equation 9 has the form of a generalized geometric series  
GG0=S=ar+ar2+=ar1rfor  |r|<1,
(10)
where I identify a=G0 and r=VG0 in equation 9, and  
S=S1+S2+S3+,
(11)
where the portion of S that is linear, quadratic, in r is  
S1=ar,S2=ar2,
(12)
and the sum is  
S=ar1r,for  |r|<1.
(13)
Solving equation 13 for r, in terms of S/a produces the inverse geometric series  
r=S/a1+S/a=S/a(S/a)2+(S/a)3+=r1+r2+r3+,when  |S/a|<1,
(14)
where rn is the portion of r that is nth order in S/a. When S is a geometric power series in r, then r is a geometric power series in S. The former is the forward series, and the latter is the inverse series. That is exactly what the inverse series represents: the inverse geometric series of the forward series equation 9. This is the simplest prototype of an inverse series for r, i.e., the inverse of the forward geometric series for S.
For the seismic inverse problem, I associate S with the measured data (see, e.g., Weglein et al., 2003)  
S=(GG0)ms=Data,
(15)
and the forward and inverse series follow from treating the forward solution as S in terms of V, and the inverse solution as V in terms of S (where S corresponds to the measured values of GG0). The inverse series is the analog of equation 14, where r1,r2, are replaced with V1,V2,:  
V=V1+V2+V3+,
(16)
where Vn is the portion of V that is nth order in the measured data D. Equation 9 is the forward-scattering series, and equation 16 is the ISS. The identity (equation 6) provides a generalized geometric forward series, a very special case of a Taylor series. A Taylor series of a function S(r) 
S(r)=S(0)+S(0)r+S(0)r22+ands(r)=S(r)S(0)=S(0)r+S(0)r22+,
(17)
whereas the geometric series is  
S(r)S(0)a=ar+ar2+.
(18)
The Taylor series equation 17 reduces to the special case of a geometric series equation 18 if  
S(0)=S(0)=S(0)2==a.
(19)

The geometric series equation 18 has an inverse series, whereas the Taylor series equation 17 does not. In general, a Taylor series does not have an inverse series. That is the reason that inversionists committed to a Taylor series starting point adopt the indirect linear updating approach, where a linear approximate Taylor series is inverted. They attempt through updating to make the linear form an ever more accurate approximate — and its premise and justification is entirely indirect and hence ad hoc — in the sense that some sort of iterative linear updating of a reference medium and model matching seek to satisfy a property that a solution might “reasonably” satisfy.

The relationship 9 provides a geometric forward series that honors equation 6 in contrast to a truncated Taylor series that does not.

All conventional current mainstream parameter estimation inversion, including iterative linear inversion, AVO, and FWI, are based on a forward Taylor series description of given data (where the chosen data can often be fundamentally and intrinsically inadequate from a direct inversion perspective), that do not honor and remain consistent with the identity equation 9.

Solving a forward problem in an inverse sense is not the same as solving an inverse problem directly

I will show that, in general, solving a forward problem in an inverse sense is not the same as solving an inverse problem directly. The exception is when the exact direct inverse is linear, as for example, in the theory of wave-equation migration (see, e.g., Claerbout, 1971; Stolt, 1978; Stolt and Weglein, 2012; Weglein et al., 2016). For wave-equation migration, given a velocity model, the migration and structure map output is a linear function of the input recorded reflection data.

To explain the latter statement, if I assume S=ar (i.e., that there is an exact linear forward relationship between S and r), then r=S/a is solving the inverse problem directly. In that case, solving the forward problem in an inverse sense is the same as solving the inverse problem directly; i.e., it provides a direct inverse solution.

However, if the forward exact relationship is nonlinear, for example,  
Sn=ar+ar2++arn,Snarar2arn=0,
(20)
and solving the forward problem 20 in an inverse sense for r will have n roots, r1,r2,,rn. As n, the number of roots . However, from the direct nonlinear forward problem S=ar/(1r), I found that the direct inverse solution r=S/(a+S) has one real root.

This discussion above provides an extremely simple, transparent, and compelling illustration of how solving a forward problem in an inverse sense is not the same as solving the inverse problem directly when there is a nonlinear forward and nonlinear inverse problem. The difference between solving a forward problem in an inverse sense (e.g., using equation 9 to solve for V) and solving an inverse problem directly (e.g., equations 2123) is much more serious, substantive, and practically significant the further I move away from a scalar single component acoustic framework. For example, it is hard to overstate the differences when examining the direct and indirect inversion of the elastic heterogeneous wave equation for earth mechanical properties and the consequences for structural and amplitude analysis and interpretation. This is a central flaw in many inverse approaches, including AVO and FWI (see Weglein, 2013).

The expansion of V in equation 16, in terms of G0 and D=(GG0)ms, the ISS (Weglein et al., 2003) can be obtained as  
G0V1G0=D,
(21)
 
G0V2G0=G0V1G0V1G0,
(22)
 
G0V3G0=G0V1G0V1G0V1G0G0V1G0V2G0G0V2G0V1G0,
(23)
To illustrate how to solve equations 2123, for V1, V2, and V3, consider the marine case with L0 corresponding to a homogeneous reference medium of water. Here, G0 is the Green’s function for propagation in water; D is the data measured, for example, with towed streamer acquisition; G is the total field that the hydrophone receiver records on the measurement surface; and G0 is the field that the reference wave (due to L0) would record at the receiver. The differential operator V then represents the difference between earth properties L and water properties L0. The solution for V is found using  
V=V1+V2+V3+,
(24)
where Vn is the portion of V that is nth order in the data D. Substituting equation 24 into the forward series equation 9, then evaluating equation 9 on the measurement surface and setting terms that are equal order in the data equal, I find equations 2123. Solving equation 21 for V1 involves the data D and G0 (water-speed propagator) and solving for V1 is analytic, and corresponds to a prestack water-speed Stolt f-k migration of the data D.

Hence, solving for V1 involves an analytic water-speed f-k migration of data D. Solving for V2 from equation 22 involves the same water-speed analytic Stolt f-k migration of G0V1G0V1G0, a quantity that depends on V1 and G0, where V1 depends on data and water speed and G0 is the water-speed Green’s function. Each term in the series produces Vn as an analytic Stolt f-k migration of a new “effective data,” where the effective data, the right side of equations 2123, are multiplicative combinations of factors that only depend on the data D and G0. Hence, every term in the ISS is directly computed in terms of data and water speed. That is the direct nonlinear inverse solution.

There are closed-form inverse solutions for a 1D earth and a normal incident plane wave (see, e.g., Ware and Aki, 1969), but the ISS is the only direct inverse method for a multidimensional subsurface.

The ISS provides a direct method for obtaining the subsurface properties contained within the differential operator L, by inverting the series order-by-order to solve for the perturbation operator V, using only the measured data D and a reference Green’s function G0, for any assumed earth model type. Equations 2123 provide V in terms of V1,V2,, and each of the Vi is computable directly in terms of D and G0. There is one equation (equation 21) that exactly produces V1, and V1 is the exact portion of V that is linear in the measured data D. The inverse operation to determine V1,V2,V3, is analytic, and it never is updated with band-limited data D. The band-limited nature of D never enters an updating process as occurs in iterative linear inversion, nonlinear AVO, and FWI.

The ISS and isolated task subseries

I can imagine that a set of tasks needs to be achieved to determine the subsurface properties V from recorded seismic data D. These tasks are achieved within equations 2123. The inverse tasks (and processing objectives) that are within a direct inverse solution are (1) free-surface multiple removal, (2) internal multiple removal, (3) depth imaging, (4) Q compensation without Q, and (5) nonlinear direct parameter estimation. Each of these five tasks has its own task-specific subseries from the ISS for V1,V2,, and each of those tasks is achievable directly and without subsurface information (see, e.g., Weglein et al., 2003, 2012; Innanen and Lira, 2010). In Appendix A, I review the details of equations 2123 for a 2D heterogeneous isotropic elastic medium.

Direct inverse and indirect inverse

Because iterative linear inversion is the concept and thinking behind many inverse approaches, I determined to make explicit the difference between that approach and a direct inverse method. The direct 2D elastic isotropic inverse solution described in Appendix A is not iterative linear inversion. Iterative linear inversion starts with equation 21. In that approach, I solve for V1 and then change the reference medium iteratively. The new differential operator L0 and the new reference medium G0 satisfy  
L0=L0V1andL0G0=δ.
(25)
In the indirect iterative linear approach, all steps basically relate to the linear relationship equation 21 with a new reference background medium, with differential operator L0 and a new reference Green’s function G0, where in terms of the new updated reference L0 equation 21 becomes  
G0V1G0=D=(GG0)ms,
(26)
where V1 is the portion of V linear in data (GG0)ms. We can continue to update L0 and G0, and we hope that the indirect procedure is solving for the perturbation operator V. In contrast, the direct inverse solution (equations 16 and A-6) calls for a single unchanged reference medium for computing V1,V2,. For a homogeneous reference medium, V1,V2, are each obtained by a single unchanged analytic inverse. We remind ourselves that the inverse to find V1 from data is the same exact unchanged analytic inverse operation to find V2,V3, from equations 21, 22, , which is completely distinct and different from equations 25 and 26 and higher iterates.

For ISS direct inversion, there are no numerical inverses, no generalized inverses, no inverses of matrices that are computed from and contain noisy band-limited data. The latter issue is terribly troublesome and difficult and is a serious practical problem, which does not exist or occur with direct ISS methods. The inverse of operators that contain and depend on band-limited noisy data is a central and intrinsic characteristic and practical pitfall of indirect methods, model matching, updating, and iterative linear inverse approaches (e.g., AVO and FWI).

Are there any circumstances in which the indirect iterative linear inversion and the direct ISS parameter estimation would be equivalent?

Are there any circumstances in which the ISS direct parameter inversion subseries would be equivalent to and correspond to the indirect iterative linear approach? Let us consider the simplest acoustic single-reflector model and a normal incident plane-wave reflection data experiment with ideal full band-width perfect data. Let the upper half-space have velocity c0 and the lower half-space have velocity c1 and then analyze these two methods (direct ISS parameter estimation and indirect iterative linear inversion) to use the reflected data event to determine the velocity of the lower half-space, c1. Yang and Weglein (2015) examine and analyze this problem and compare the results of the direct ISS method and the indirect iterative linear inversion. They show that the direct ISS inversion to estimate c1 converged to c1 under all circumstances and all values of c0 and c1. In contrast, the indirect linear iterative inversion had a limited range of values of c0 and c1 where it converged to c1, and in that range, it converged much slower than the direct ISS parameter estimation for c1. The iterative linear inverse simply shut down and failed when the reflection coefficient R was greater than 1/4 (see Appendix B and Yang, 2014).

The direct ISS parameter estimation method converged to c1 for any value of the reflection coefficient R. Hence, under the simplest possible circumstance, and providing the iterative linear method with an analytic Fréchet derivative, as a courtesy from and a gift delivered to the linear iterative from the ISS direct inversion method, the ranges of usefulness, validity, and relative effectiveness were never equivalent or comparable. With band-limited data and more complex earth models (e.g., elastic multiparameter), this gap in the range of validity, usefulness, and effectiveness will necessarily widen (see Zhang, 2006; Weglein, 2013). The indirect iterative linear inversion and the direct ISS parameter-estimation method are never equivalent, and there are absolutely no simple or complicated circumstances in which they are equally effective. The distinct ISS free-surface-multiple elimination subseries and internal-multiple attenuation subseries are not only not dependent on subsurface properties, but they are precisely the same unchanged algorithms for any earth model type.

There was an earlier time when free-surface multiples were modeled and subtracted. Multiple-removal methods have moved on. Parameter-estimation methods continue to be firmly connected to model matching and subtraction. That stark and immense difference between iterative linear updating model matching and the direct inversion inverse scattering methods is an essential point to consider and comprehend for those interested in understanding these methodologies and their seismic processing and interpretation consequences and value. It is not conceivable to even formulate an iterative linear model matching method that is not dependent on a specified model type — let alone to compare it with ISS model-type-independent algorithms.

Direct ISS parameter inversion: A time-lapse application

The direct inverse ISS elastic parameter estimation method (equation A-6) was successfully applied (Zhang et al., 2006) in a time-lapse sense to discriminate between pressure and fluid saturation changes. Traditional time-lapse estimation methods were unable to predict and match that direct inversion ISS discrimination.

Further substantive differences between iterative linear model matching inversion and direct inversion from the Lippmann-Schwinger equation and the ISS

The difference between iterative linear and the direct inverse of equation A-6 is much more substantive and serious than merely a different way to solve G0V1G0=D (equation 21), for V1. If equation 21 is someone’s entire basic theory, you can mistakenly think that  
D^PP=G^0PV^1PPG^0P
(27)
is sufficient to update (generalizing equations 25 and 26)  
D^PP=G^0PV^1PPG^0P.
(28)
Please note that ^ indicates that variables are transformed to PS space. This step loses contact with and violates the basic operator identity G=G0+G0VG for the elastic wave equation. The fundamental identity G=G0+G0VG for the elastic wave equation is a nonlinear multiplicative matrix relationship. For the forward and inverse series, the input and output variables are matrices. The inverse solution for a change in an earth mechanical property has a nonlinear coupled dependence on all the data components  
(DPPDPSDSPDSS),
(29)
in 2D and the P, SH, SV 3×3 generalization in 3D (Stolt and Weglein, 2012, chapter 7).
A unique expansion of VG0 in orders of measurement values of (GG0) is  
VG0=(VG0)1+(VG0)2+.
(30)

The scattering-theory equation allows that forward series form the opportunity to find a direct inverse solution. Substituting equation 30 into equation 9 and setting the terms of equal order in the data to be equal, I have D=G0V1G0, where the higher order terms are V2,V3,, as given in Weglein et al. (2003, p. R33, equations 714).

For the elastic equation, V is a matrix and the relationship between the data and V1 is  
(DPPDPSDSPDSS)=(G0P00G0S)(V1PPV1PSV1SPV1SS)(G0P00G0S),
(31)
 
V1=(V1PPV1PSV1SPV1SS),
(32)
 
V=(VPPVPSVSPVSS),
(33)
 
V=V1+V2+,
(34)
where V1,V2 are linear, quadratic contributions to V in terms of the data  
D=(DPPDPSDSPDSS).
(35)
The changes in elastic properties and density are contained in  
V=(VPPVPSVSPVSS),
(36)
and that leads to direct and explicit solutions for the changes in mechanical properties in orders of the data  
D=(DPPDPSDSPDSS),
(37)
 
Δγγ=(Δγγ)1+(Δγγ)2+,
(38)
 
Δμμ=(Δμμ)1+(Δμμ)2+,
(39)
 
Δρρ=(Δρρ)1+(Δρρ)2+,
(40)
where γ, μ, and ρ are the bulk modulus, shear modulus, and density, respectively.

The ability of the forward series to have a direct inverse series derives from (1) the identity among G, G0, and V provided by the scattering equation and then (2) the recognition that the forward solution can be viewed as a geometric series for the data D, in terms of VG0. The latter derives the direct inverse series for VG0 in terms of the data.

Viewing the forward problem and series as the Taylor series  
D(m)=D(m0)+D(m0)Δm+D(m0)2Δm2+,
(41)
in which the derivatives are Fréchet derivatives, in terms of Δm, does not offer a direct inverse series, and hence there is no choice but to solve the forward series in an inverse sense. It is that fact that results in all current AVO and FWI methods being modeling methods that are solved in an inverse sense. Among references that solve a forward problem in an inverse sense in P-wave AVO are Clayton and Stolt (1981), Shuey (1985), Stolt and Weglein (1985), Boyse and Keller (1986), Stolt (1989), Beylkin and Burridge (1990), Castagna and Smith (1994), Goodway et al. (1997), Burridge et al. (1998), Smith and Gidlow (2000), Foster et al. (2010), and Goodway (2010). The intervention of the explicit relationship among G, G0, and V (the scattering equation) in a Taylor series-like form produces a geometric series and a direct inverse solution.
The linear equations are  
(D^PPD^PSD^SPD^SS)=(G^0P00G^0S)(V^1PPV^1PSV^1SPV^1SS)(G^0P00G^0S),
(42)
 
D^PP=G^0PV^1PPG^0P,
(43)
 
D^PS=G^0PV^1PSG^0S,
(44)
 
D^SP=G^0SV^1SPG^0P,
(45)
 
D^SS=G^0SV^1SSG^0S,
(46)
 
D˜PP(kg,0;kg,0;ω)=14(1kg2νg2)a˜ρ(1)(2νg)14(1+kg2νg2)a˜γ(1)(2νg)+2kg2β02(νg2+kg2)α02a˜μ(1)(2νg),
(47)
 
D˜PS(νg,ηg)=14(kgνg+kgηg)a˜ρ(1)(νgηg)β022ω2kg(νg+ηg)(1kg2νgηg)a˜μ(1)(νgηg),
(48)
 
D˜SP(νg,ηg)=14(kgνg+kgηg)a˜ρ(1)(νgηg)+β022ω2kg(νg+ηg)(1kg2νgηg)a˜μ(1)(νgηg),and
(49)
 
D˜SS(kg,ηg)=14(1kg2ηg2)a˜ρ(1)(2ηg)[ηg2+kg24ηg22kg2ηg2+kg2]a˜μ(1)(2ηg),
(50)
where aγ(1), aμ(1), and aρ(1) are the linear estimates of the changes in bulk modulus, shear modulus, and density, respectively. Here, kg is the Fourier conjugate to the receiver position xg and νg and ηg are the vertical wavenumbers for the P- and S-reference waves, respectively, where  
νg2+kg2=ω2α02,
(51)
 
ηg2+kg2=ω2β02,
(52)
and α0 and β0 are the P- and S-velocities in the reference medium, respectively. The direct quadratic nonlinear equations are  
(G^0P00G^0S)(V^2PPV^2PSV^2SPV^2SS)(G^0P00G^0S)=(G^0P00G^0S)(V^1PPV^1PSV^1SPV^1SS)(G^0P00G^0S)(V^1PPV^1PSV^1SPV^1SS)(G^0P00G^0S),
(53)
 
G^0PV^2PPG^0P=G^0PV^1PPG^0PV^1PPG^0PG^0PV^1PSG^0SV^1SPG^0P,
(54)
 
G^0PV^2PSG^0S=G^0PV^1PPG^0PV^1PSG^0SG^0PV^1PSG^0SV^1SSG^0S,
(55)
 
G^0SV^2SPG^0P=G^0SV^1SPG^0PV^1PPG^0PG^0SV^1SSG^0SV^1SPG^0P,
(56)
 
G^0SV^2SSG^0S=G^0SV^1SPG^0PV^1PSG^0SG^0SV^1SSG^0SV^1SSG^0S.
(57)

Because V^1PP relates to D^PP, V^1PS relates to D^PS, and so on, the four components of the data will be coupled in the nonlinear elastic inversion. I cannot perform the direct nonlinear inversion without knowing all components of the data. Thus, the direct nonlinear solution determines the data needed for a direct inverse. That, in turn, defines what a linear estimate means. That is, a linear estimate of a parameter is an estimate of a parameter that is linear in data that can directly invert for that parameter. Because DPP, DPS, DSP, and DSS are needed to determine aγ, aμ, and aρ directly, a linear estimate for any one of these quantities requires simultaneously solving equations 4750 (for further details, see, e.g., Weglein et al., 2009).

Those direct nonlinear formulas are like the direct solution for the quadratic equation mentioned above and solve directly and nonlinearly for changes in the velocities, α, β, and the density ρ in a 1D elastic earth. Stolt and Weglein (2012) present the linear equations for a 3D earth that generalize equations 4750. Those formulas prescribe precisely what data you need as input, and they dictate how to compute those sought-after mechanical properties, given the necessary data. There is no search or cost function, and the unambiguous and unequivocal data needed are full-multicomponent data — PP, PS, SP, and SS — for all traces in each of the P- and S-shot records. The direct algorithm determines first the data needed and then the appropriate algorithms for using those data to directly compute the sought-after changes in the earth’s mechanical properties. Hence, any method that calls itself inversion (let alone full-wave inversion) for determining changes in elastic properties, and in particular, the P-wave velocity α, and that inputs only P-data, is more off base, misguided, and lost than the methods that sought two or more functions of depth from a single trace. You can model-match P-data ad nauseum, which takes a lot of computational effort and people with advanced degrees in math and physics computing Fréchet derivatives, and it requires sophisticated LP-norm cost functions and local or global search engines, so it must be reasonable, scientific, and worthwhile. Why can I not use just PP-data to invert for changes in VP, VS, and density because Zoeppritz says that I can model PP from those quantities and because I have, using PP-data with angle variation, enough dimension? As stated above, data dimension is good, but it is not good enough for a direct inversion of those elastic properties.

Adopting equations 27 and 28 as in AVO and FWI, there is a violation of the fundamental relationship between changes in a medium and changes in a wavefield, G=G0+G0VG, which is as serious as considering problems involving a right triangle and violating the Pythagorean theorem. That is, iteratively updating PP data with an elastic model violates the basic relationship between changes in a medium V and changes in the wavefield GG0 for the simplest elastic earth model.

This direct inverse method for parameter estimation provides a platform for amplitude analysis and a solid framework and direct methodology for the goals and objectives of indirect methods such as AVO and FWI. A direct method for the purposes of amplitude analysis provides a method that derives from, respects, and honors the fundamental identity and relationship G=G0+G0VG. Iteratively inverting multicomponent data has the correct data, but it does not correspond to a direct inverse algorithm. To honor G=G0+G0VG, you need the data and the algorithm that the direct inverse prescribes. Not recognizing the message that an operator identity and the elastic wave equation unequivocally communicate is a fundamental and significant contribution to the gap in effectiveness in current AVO and FWI methods and application (equation A-6). This analysis generalizes to 3D with P, SH, and SV data.

The role of direct and indirect methods

There is a role for direct and indirect methods in practical real-world applications. In our view, indirect methods are to be called upon for recognizing that the world is more complicated than the physics that we assume in our models and methods. For the part of the world that you are capturing in your model and physics, nothing compares to direct methods for clarity and effectiveness. An optimal indirect method would seek to satisfy a cost function that derives from a property of the direct method. In that way, the indirect and direct methods would be aligned, consistent, and cooperative for accommodating the part of the world described by your physical model (with a direct inverse method) and the part that is outside (with an indirect method).

The indirect method of model matching primaries and multiples (so-called FWI)

All model matching inverse approaches are indirect methods. Iterative linear inversion model matching is an indirect search methodology, which is ad hoc and without a firm and solid foundation and theoretical and conceptual framework. Nevertheless, we can imagine and understand that model matching primaries and multiples, rather than only primaries, could improve upon matching only primaries. However, model matching primaries and multiples remains ad hoc and indirect and is always on much shakier footing than direct inversion for the same inversion goals and objectives. Direct ISS inversion for parameter estimation only requires and inputs primaries.

For all multidimensional seismic applications, the only direct inverse solution is provided by the operator identity equation 6 and is in the form of a series of equations 2123, the ISS (Weglein et al., 2003). It can achieve all processing objectives within a single framework and a single set of equations 2123 without requiring any subsurface information. There are distinct isolated-task inverse scattering subseries derived from the ISS, which can perform free-surface multiple removal (Carvalho et al., 1992; Weglein et al., 1997), internal multiple removal (Araújo et al., 1994; Weglein et al., 2003), depth imaging (e.g., Shaw, 2005; Liu, 2006; Weglein et al., 2012), parameter estimation (Zhang, 2006; Li, 2011; Liang, 2013; Yang and Weglein, 2015), and Q compensation without needing, estimating, or determining Q (Innanen and Weglein, 2007; Lira, 2009; Innanen and Lira, 2010), and each achieves its objective directly and without subsurface information. The direct inverse solution (e.g., Weglein et al., 2003, 2009) provides a framework and a firm math-physics foundation that unambiguously defines the data requirements and the distinct algorithms to perform each and every associated task within the inverse problem, directly and without subsurface information.

Having an ad hoc, indirect method as the starting point places a cloud over issue identification when less-than-satisfactory results arise with field data. In addition, we saw that direct inversion parameter estimation has a significantly lower dependence on the low-frequency data components in comparison with indirect methods such as nonlinear AVO and FWI.

Only a direct solution can provide algorithmic clarity, confidence, and effectiveness. The current industry-standard AVO and FWI, using variants of model-matching and iterative linear inverse, are indirect methods, and iteratively linearly updating P data or multicomponent data (with or without multiples) does not correspond to, and will not produce, a direct solution.

All direct inverse methods for structural determination and amplitude analysis require only primaries

In Weglein (2016), the role of primaries and multiples in imaging is examined and analyzed. The most capable and interpretable migration method derives from predicting a source and receiver experiment at depth. For data consisting of primaries and multiples, a discontinuous velocity model is needed to achieve that predicted experiment at depth. With that discontinuous velocity model, free-surface and internal multiples play no role in the migration and the exact same image results with or without multiples (see Weglein, 2016). For a smooth velocity model, multiples will result in false and misleading images and must be removed before the migration and migration-inversion of primaries.

In Weglein et al. (2003), ISS direct depth imaging (without a velocity model or subsurface information) removes free-surface and internal multiples prior to the distinct subseries that input primaries and perform depth imaging and amplitude analysis, respectively, each directly and without subsurface information and only using and requiring primaries.

Hence, all direct inversion methods, those with and those without subsurface/velocity information, require only primaries for complete structural determination and amplitude analysis. Methods that seek to use multiples to address issues from less than a complete acquisition of primaries are seeking an appropriate image of an unrecorded primary.

Indirect methods are ad hoc without a clear or firm math-physics foundation and framework, and they start without knowing whether “the indirect solution” is in fact a solution. A more complete or fuller data set being matched between model data and field data, each with primaries and multiples, could at times improve upon matching only primaries, but the entire approach is indirect and ad hoc with or without multiples, and it lacks the benefits of a direct method. With indirect methods, there is no framework and theory to rely on, and no confidence that a solution is forthcoming under any circumstances.

If I seek the parameters of an elastic heterogeneous isotropic subsurface, then the differential operator in the operator identity is the differential operator that occurs in the elastic, heterogeneous, isotropic wave equation. From 40 years of AVO and amplitude analysis application in the petroleum industry, the elastic isotropic model is the baseline minimally realistic and acceptable earth model type for amplitude analysis, for example, for AVO and FWI. Then, taking the operator identity (called the Lippmann-Schwinger, or scattering theory, equation) for the elastic-wave equation, I can obtain a direct inverse solution for the changes in the elastic properties and density. The direct inverse solution specifies the data required and the algorithm to achieve a direct parameter estimation solution. In this paper, I explain how this methodology differs from all current AVO and FWI methods, which are, in fact, forms of model matching. Multicomponent data consisting of only primaries are needed for a direct inverse solution for subsurface properties. This paper focuses on one specific inverse task, parameter estimation, within the overall and broader set of inversion objectives and tasks. Furthermore, the impact of band-limited data and noise are discussed and compared for the direct ISS parameter estimation and indirect (AVO and FWI) inversion methods.

In this paper, I focused on analyzing and examining the direct inverse solution that the ISS inversion subseries provides for parameter estimation. The distinct issues of (1) data requirements, (2) model type, and (3) inversion algorithm for the direct inverse are all important (Weglein, 2015b). For an elastic heterogeneous medium, I show that the direct inverse requires multicomponent/PS (P- and S-component) data and prescribes how that data are used for a direct parameter estimation solution (Zhang and Weglein, 2006).

Conclusion

In this paper, I describe, illustrate, and analyze the considerable conceptual, substantive, and practical benefit and added value that a direct parameter inversion from the ISS provides in comparison with all current indirect inverse methods (e.g., AVO and FWI) for amplitude analysis goals and objectives. A direct method provides (1) a solution that we (the seismic industry) can have confidence that it is in fact solving the defined problem of interest and (2) in addition, when the method does not improve the drilling decisions, then we know that the issue is that the problem of interest is not the problem that we need to be interested in. On the other hand, indirect methods such as AVO and FWI have a plethora of approaches and paths, and when less-than-satisfactory results occur, we do not know whether the issue is the chosen problem of interest or the choice among innumerable indirect solutions, or both.

All scientific methods make assumptions — and seismic processing and interpretation methods are no exception. When the assumptions behind seismic methods are satisfied, the methods are useful and effective and can support successful drill decisions. When the assumptions are not satisfied, the methods can have difficulty or can fail. The latter breakdown can contribute to unsuccessful ill-informed drill decisions, dry-hole drilling, or suboptimal appraisal and development wells.

The objective of seismic research is to provide new and effective toolbox capability for processing and interpretation that will improve the drill success rate and reduce dry-hole and suboptimal drilling decisions. Toward that end, the starting point in seismic research is to identify the outstanding prioritized problems and challenges that need to be addressed and solved.

The ability to clearly and unambiguously define the origin and root cause behind seismic issues, problems, breakdown, and challenges is an essential and critically important step in designing and executing a strategy to provide new and more capable methods to the seismic processing and interpretation toolbox.

Direct inversion methods can provide that problem definition and clarity. They are also unique in providing the confidence that the problem of interest is actually being addressed. For ISS parameter estimation, although the recorded data are of course band limited, the band-limited data are never used to compute the updated inverse operator for the next iterated linear step because the inverse operator is fixed and analytic for every term in the ISS. That is one of several important and substantive differences pointed out in this paper between the direct inverse ISS parameter estimation method and all indirect inversion methods, e.g., AVO and FWI. I provide an explicit analytic example and comparison between direct ISS parameter estimation and the indirect linear updating model matching concepts behind AVO and FWI.

All seismic processing methods depend on the amplitude and phase of seismic data. Different processing methods that seek to achieve a certain specific processing goal can have different relative sensitivities to noise and bandwidth. Amplitude analysis for determining earth mechanical property changes is one of the most sensitive. Methods that achieve seismic goals as a sequence of separate intermediate steps have a natural advantage over methods that seek to combine goals. Achieving an intermediate, easier goal that is less demanding can significantly enhance the ability to achieve the subsequent more demanding seismic processing objectives. The indirect methods that seek to locate structure and identify changes in earth mechanical properties at once have a terrible dependence on missing low-frequency data. However, if I first locate a structure by wave-equation migration (a process that is insensitive to missing low frequency data), then in principle, I can determine the earth mechanical property changes with a single frequency within the bandwidth. The ISS direct amplitude analysis method described, exemplified, tested, and compared in this paper assumes that a set of less-daunting seismic processing tasks, using an ISS task specific subseries, has been achieved (e.g., multiple removal, depth imaging) before this task is undertaken. To have a fair comparison, the indirect model matching method is tested with a data with a well-located single reflector, and hence there are no imaging issues or multiples in the problem. That allows a pristine, clear, and definitive comparison of the amplitude analysis — parameter estimation function of the prototype direct ISS method and the corresponding indirect model-matching iterative updating approach. There are important issues of resolution and illumination, which will impact the results of this paper, with advances in migration theory and algorithms that avoid all high-frequency approximations in the imaging principles and wave-propagation models that can improve resolution and illumination.

Direct and indirect methods can play an important role and function in seismic processing, in which the former accommodates and addresses the assumed physics and the latter provides a channel for real-world phenomena beyond the assumed physics. Both are called for within a comprehensive and effective seismic processing and interpretation strategy.

Acknowledgments

The author is grateful to all M-OSRP sponsors for their support and encouragement of this research. J. Mayhan is thanked for his invaluable assist in preparing this paper. The author is deeply grateful and appreciative to H. Bui for her encouragement to write and submit this paper. The author thanks B. Goodway and J. Zhang for their detailed, thoughtful, and constructive reviews. All of their suggestions enhanced and benefited this paper. The author is deeply grateful to H. Zhang, H. Liang, X. Li, and S. Jiang for their historic contributions to direct ISS elastic parameter estimation.

The operator identity and direct inverse solution for a 2D heterogeneous isotropic elastic medium

I describe the forward and direct inverse method for a 2D elastic heterogeneous earth (see Zhang, 2006).

The 2D elastic wave equation for a heterogeneous isotropic medium (Zhang, 2006) is  
Lu=(fxfz)andL^(ϕPϕS)=(FPFS),
(A-1)
where u, fx, and fz are the displacement and forces in displacement coordinates and ϕP, ϕS, and FP, FS are the P and S-waves and the force components in P- and S-coordinates, respectively. The operators L and L0 in the actual and reference elastic media are  
L=[ρω2(1001)+(xγx+zμzx(γ2μ)z+zμxz(γ2μ)x+xμzzγz+xμx)],
(A-2)
 
L0=[ρω2(1001)+(γ0x2+μ0z2(γ0μ0)xz(γ0μ0)xzμ0x2+γ0z2)],
(A-3)
and the perturbation V is  
VL0L=[aρω2+α02xaγx+β02zaμzz(α02aγ2β02aμ)x+β02xaμzx(α02aγ2β02aμ)z+β02zaμxaρω2+α02zaγz+β02xaμx],
(A-4)
where the quantities aρρ/ρ01, aγγ/γ01, and aμμ/μ01 are defined in terms of the bulk modulus, shear modulus, and density (γ0, μ0, ρ0, γ, μ, ρ) in the reference and actual media, respectively.
The forward problem is found from the identity equation 9 and the elastic wave equation A-1 in PS-coordinates as  
G^G^0=G^0V^G^0+G^0V^G^0V^G^0+,(D^PPD^PSD^SPD^SS)=(G^0P00G^0S)(V^PPV^PSV^SPV^SS)(G^0P00G^0S)+(G^0P00G^0S)(V^PPV^PSV^SPV^SS)(G^0P00G^0S)(V^PPV^PSV^SPV^SS)(G^0P00G^0S)+,
(A-5)
and the inverse solution, equations 2123, for the elastic equation A-1 is  
(D^PPD^PSD^SPD^SS)=(G^0P00G^0S)(V^1PPV^1PSV^1SPV^1SS)(G^0P00G^0S),(G^0P00G^0S)(V^2PPV^2PSV^2SPV^2SS)(G^0P00G^0S)=(G^0P00G^0S)(V^1PPV^1PSV^1SPV^1SS)(G^0P00G^0S)(V^1PPV^1PSV^1SPV^1SS)(G^0P00G^0S),
(A-6)
where V^PP=V^1PP+V^2PP+V^3PP+ and any one of the four matrix elements of V requires the four components of the data  
(D^PPD^PSD^SPD^SS).
(A-7)
The 3D heterogeneous isotropic elastic generalization of the above 2D forward and direct inverse elastic isotropic method begins with the linear 3D form found in Stolt and Weglein (2012, p. 159).

In summary, from equation A-5, D^PP can be determined in terms of the four elements of V. The four components V^PP, V^PS, V^SP, and V^SS require the four components of D. That is what the general relationship G=G0+G0VG requires; i.e., a direct nonlinear inverse solution is a solution order-by-order in the four matrix elements of D (in 2D). The generalization of the forward series equation A-5 and the inverse series equation A-6 for a direct inversion of an elastic isotropic heterogeneous medium in 3D involves the 3×3 data, D, and V matrices in terms of P, SH, and SV data and start with the linear G0V1G0=D (Stolt and Weglein, 2012, p. 179).

Numerical examples for a 1D normal incident wave on an acoustic medium

Numerical examples for a 1D normal incident wave on an acoustic medium are shown in this section. First, I examine and compare the convergence of the ISS direct inversion and iterative inversion. Second, the rate of convergence of the ISS inversion subseries is examined and studied using an analytic example, where the ISS method converges and the iterative linear method does not and where both methods converge.

The operator identity for a 1D acoustic medium

For a normal incidence plane wave on a 1D acoustic medium (where only the velocity is assumed to vary), the model I consider here consists of two half-spaces with acoustic velocities c0 and c1 and an interface located at z=a as shown in Figure B-1. If I put the source and receiver on the surface, z=0, the pressure wave  
D(t)=Rδ(t2a/c0)
(B-1)
will be recorded, where the reflection coefficient R=(c1c0)/(c1+c0). For this example, D(t) is the only input to the direct ISS inverse and the iterative inversion methods. Because I will assume knowledge of the velocity in the upper half-space, c0, the location of the reflector at z=a is not an issue. I will focus on only determining the change of velocity across the reflector at z=a. The operators L0 and L in the reference and actual acoustic media are  
L0=d2dz2+ω2c02andL=d2dz2+ω2c2(z),
(B-2)
and I characterize the velocity perturbation as  
α(z)1c02c2(z).
(B-3)
The perturbation V (Weglein et al., 2003) can be expressed as  
V(z)=L0L=ω2c02ω2c2(z)=k02α(z),
(B-4)
where ω is the angular frequency and k0=ω/c0. The functions c0 and c(z) are the reference and local acoustic velocity, respectively. Therefore, the inverse series of V (equation 16) becomes  
α(z)=α1(z)+α2(z)+α3(z)+.
(B-5)
That is,  
V1=k02α1,V2=k02α2,.
(B-6)
From the ISS (equations 2123), Shaw and Weglein (2004) isolate the leading order imaging subseries and the direct nonlinear inversion subseries.
In this section, I will focus on studying the convergence properties of the ISS inversion subseries. The inversion-only terms isolated from the ISS (Zhang, 2006; Li, 2011) are  
α(z)=α1(z)12α12(z)+316α13(z)+.
(B-7)
For a 1D normal incidence case, the linear equation 21 solves for α1 in terms of the single trace data D(t) (Shaw and Weglein, 2004) as  
α1(z)=4zD(z)dz,
(B-8)
where z=c0t/2. For a single reflector, inserting data D (equation B-1) gives  
α1=4RH(za),
(B-9)
where R is the reflection coefficient R=(c1c0)/(c1+c0) and H is the Heaviside function. When z>a, substituting α1 into equation B-7, the ISS direct nonlinear inversion subseries in terms of R can be written as (where α is the magnitude of α(z) for z>a)  
α=4R8R2+12R3+=4Rn=0(n+1)(R)n.
(B-10)
After solving for α, the inverted velocity c(z) can be obtained through c1=c0(1α)1/2 (equation B-4).
Considering the convergence property of the series for α or the inversion subseries, I can calculate the ratio test  
|αn+1αn|=|(n+2)(R)n+1(n+1)(R)n|=|n+2n+1R|.
(B-11)
If limn|(αn+1/αn)|<1, this subseries converges absolutely. That is,  
|R|<limnn+1n+2=1.
(B-12)
Therefore, the ISS direct nonlinear inversion subseries converges when the reflection coefficient |R| is less than one, which is always true. Hence, for this example, the ISS inversion subseries will converge under any velocity contrasts between the two media.
For the iterative linear inversion, I use the first linear estimate of α=α11 to compute the first estimate of c1=c11. Then, I choose the first estimate of c1=c0(1α11)1/2c11 as the new reference velocity, c01=c0(1α11)1/2, where α11=4R1 and R1=(c1c0)/(c1+c0). Repeating the linear process with a new reflection coefficient R2 (again exploiting the analytic inverse generously provided by ISS to benefit the iterative linear inverse approach) gives  
R2=c1c01c1+c01,  α12=4R2  and  c12=c01(1α12)1/2=c02,
(B-13)
 
Rn+1=c1c0nc1+c0n,  α1n+1=4Rn+1  and  c1n=c0n1(1α1n)1/2=c0n,
(B-14)
where α1n=nth estimate of α1 and c1n=nth estimate of c1. The questions are (1) under what conditions does c1n approach c1, and (2) when it converges, what is its rate of convergence?

From the above analysis, I can see that the ISS method for α always converges and the resulting α can be used to find c1. For the iterative linear inverse, there are values of α1, such that you cannot compute a real c11. When α11>1 and 4R>1, R>1/4 and you cannot compute an updated reference velocity and the method simply shuts down and fails. The ISS never computes a new reference and does not suffer that problem, with the series for α always converging and then outputting c1, the correct unknown velocity below the reflector.

The convergence of the ISS direct inversion and iterative inversion

In this section, I will examine and compare the convergence property of the ISS inversion (equation B-10) and the iterative linear inversion for different velocity contrasts in the 1D acoustic case. In the 1D normal incident acoustic model (Figure B-1), only one parameter (velocity) varies and a plane wave propagates into the medium. There is only a single reflector, and I assume the velocity is known above the reflector and unknown below the reflector. I will compare the convergence of the perturbation α and the inversion results by using the ISS direct nonlinear method and the iterative linear method.

With the reference velocity c0=1500m/s, two analytic examples with different velocity contrasts for c1=2000 and 3000m/s are examined. Figure B-2 shows the estimated α by the ISS method (green line) for c1=2000m/s. The red line represents the actual α that is calculated from the model. The horizontal axis represents the order of the ISS inversion subseries. The vertical axis shows the value of α. The updated estimation of α using the iterative inversion method (blue line) is shown in Figure B-3. The horizontal axis represents the iteration numbers in the iterative inversion method. From Figures B-2 and B-3, I can see that at the small velocity contrast, the estimated α by ISS method becomes the actual α after about five orders of calculation and the updated estimation of α by the iterative inversion method goes to zero as expected because after several iterations, the updated model is close to and approaching to the actual model. Figure B-4 represents the velocity estimation. The green and blue lines represent the estimated velocity by using the ISS inversion method and the iterative inversion method, respectively. We can see that at the small velocity contrast, both methods converge and produce correct velocity after five orders of iterations and the ISS inversion method converges faster than the iterative inversion method.

Figure B-5 shows the estimated α by the ISS method (green line) for c1=3000m/s. When the velocity contrast is larger, i.e., R>0.25, the iterative inversion method cannot be computable, but the ISS inversion method always converges (see the green line in Figure B-5) after the summation of more orders in computing α.

As we know, the reflection coefficient R is almost always less than 0.2 in practice, so that the ISS method and the iterative method converge, but the ISS method converges faster than the iterative method. Moreover, for more complicated circumstances (e.g., the elastic nonnormal incidence case), the difference between the ISS method and the iterative method is much greater, not just on the algorithms, but also on data requirements and on how the band-limited noisy nature of the seismic data impacts the inverse operators in the iterative method but not in the ISS method.

The rate of convergence of the ISS inversion subseries

The rate of convergence of the estimated α for the ISS inversion subseries (equation B-10) is analytically examined and studied. Because α is always convergent when R<1, the summation of this subseries (Zhang, 2006) is  
α=4Rn=0(n+1)(R)n=4R1(1+R)2.
(B-15)
If the error between the estimated and the actual α is monotonically decreasing, it means that the subseries is a term-by-term added-value improvement toward determining the actual medium properties. If this error is increasing before decreasing, it means that the estimate of α becomes worse before it gets better. The error for the first order and the error for the second order have the relation  
|αα1α2|>|αα1|,
(B-16)
i.e.,  
|4R3R2+2R3(1+R)2|>|4RR22R(1+R)2|.
(B-17)
After simplification, it gives  
R2+R1>0.
(B-18)
I can solve it and obtain the reflection coefficient R<[(15)/2]=1.618 or R>[(1+5)/2]=0.618. Therefore, when R>0.618, the error increases first. Similarly, if the error for the third order is greater than that for the second order, I get R>0.667. If the error for the fourth order is greater than that for the third order, I obtain R>0.721. In summary, when R>0.618, the error increases and the estimated α gets worse before getting better. The sum of terms in the direct inverse ISS solution (for very large contrasts) requires certain partial sums to be temporarily worse in order for the entire series to produce the correct velocity. The dashed green line in Figure B-6 shows that when the reflection coefficient R is equal to 0.618, the error for the first order is equal to the error for the second order.

As the analytic calculation, when the reflection coefficient R is smaller than 0.618, this inversion subseries gives a monotonically term-by-term added-value improvement toward determining c1. When the reflection coefficient is larger than 0.618, the ISS inversion series still converges, but the estimation of α will become worse before it gets better. Each term in the series works toward the final goal. Sometimes when more terms in the series are included, the estimation looks temporarily worse, but once it starts to improve the estimation at a specific order, the approximations never become worse again, and every single term after that order will produce an improved estimation. The locally worse partial sum behavior is, in fact, purposeful and essential for convergence to and for computing the exact velocity. The direct inverse solution fulfills its commitment to always predict c1 and not necessarily to having order-by-order improvement. The ISS direct inversion always converges in contrast to the iterative linear inverse method. This property has also been indicated by Carvalho (1992) in the free-surface-multiple elimination subseries; e.g., what appears to make a second-order free-surface multiple larger with a first-order free-surface algorithm is actually helpful and necessary for preparing the second-order multiple to be removed by the higher order terms.

Arthur Weglein holds the Hugh Roy and Lillie Cranz Cullen Distinguished University Chair in Physics, is Professor in the Department of Physics and the Department of Earth and Atmospheric Sciences, and is Director of the Mission-Oriented Seismic Research Program at the University of Houston. He received his Ph.D. from the City University of New York (CUNY/CCNY), was a research scientist at Cities Service Company, SOHIO Petroleum Company, ARCO Exploration and Production Technology, and Schlumberger Cambridge Research, and was a visiting professor at the Universidade Federal da Bahia. He served as SEG Distinguished Lecturer in 2003, he received the Townsend Harris Medal from CUNY/CCNY in 2008, the Reginald Fessenden Medal from SEG in 2010, and in 2016, he received the highest SEG honor and recognition, the Maurice Ewing Medal.

Freely available online through the SEG open-access option.