Krylov-space-based quasi-newton preconditioner for full-wavefield inversion

Abstract

A method, including: storing, in a computer memory, seismic data acquired from a seismic survey of a subsurface region; and generating, with a computer, a final subsurface physical property model of the subsurface region by processing the seismic data with an iterative full wavefield inversion method, wherein the iterative full wavefield inversion method generates the final subsurface physical property model by iteratively applying a linear solver with a preconditioner that is generated from information from one or more previous iterations of the linear solver.

Description

FIELD OF THE INVENTION

Exemplary embodiments described herein pertain generally to the field of geophysical prospecting, and more particularly to geophysical data processing. An exemplary embodiment can increase the speed of convergence of full wavefield inversion (FWI).

BACKGROUND

This section is intended to introduce various aspects of the art, which may be associated with exemplary embodiments of the present technological advancement. This discussion is believed to assist in providing a framework to facilitate a better understanding of particular aspects of the present technological advancement. Accordingly, it should be understood that this section should be read in this light, and not necessarily as admissions of prior art.

Seismic inversion is a process of extracting subsurface information from the data measured at the surface of the earth acquired during a seismic survey. In a typical seismic survey, seismic waves are generated by a source positioned at desired locations. As the source generated wave propagates through the subsurface and travels back to the receiver locations where it is recorded.

Full waveform inversion (FWI) is a seismic processing method which can potentially exploit the full seismic record including events that are treated as “noise” by standard seismic processing algorithms. FWI iteratively minimizes an objective function based on a comparison of a simulated and measured seismic records. Even with today's available high-performance computing resources, one of the biggest challenges to FWI is still the computational cost. Nevertheless, the benefit of inferring a detailed representation of the subsurface using this method is expected to outweigh this cost, and development of new algorithms and workflows that lead to faster turnaround time is a key step towards making this technology feasible for field-scale applications, allowing users to solve larger scale problems faster. The computationally most intensive component of FWI is the simulations of the forward and adjoint wavefields. The number of total forward and adjoint simulations is proportional to the number of iterations, which is typically on the order of hundreds to thousands. Any method reducing number of FWI iterations will reduce to number of forward and adjoint simulation calls and the computational run time.

The crux of any FWI algorithm can be described as follows: using a given starting subsurface physical property model, synthetic seismic data are generated, i.e. modeled or simulated, by solving the wave equation using a numerical scheme (e.g., finite-difference, finite-element, spectral element, and etc.) which typically divides the model domain into a set of nonoverlapping cells (also referred as elements or blocks). The term velocity model or geophysical property model as used herein refers to an array of numbers, where each number, which may also be called a model parameter, is a value of velocity or another geophysical property in a cell. The synthetic seismic data are compared with the field seismic data and using a norm, an error or objective function is calculated. Using this objective function and an optimization algorithm, a modified subsurface model is generated which is used to simulate a new set of synthetic seismic data. This new set of synthetic seismic data is compared with the field data to generate a new objective function. This process is repeated until the optimization algorithm satisfactorily minimizes the objective function and the final subsurface model is generated. A global or local optimization method is used to minimize the objective function and to update the subsurface model. Further details regarding FWI can be found in U.S. Patent Publication 2011/0194379 to Lee et al., the entire contents of which are hereby incorporated by reference.

Common FWI methods iteratively minimize the objective function which is subject to the wavefield propagation—the physics of the problem. A (nonlinear) iteration i of FWI involves the following two steps: (1) compute a search direction for the current model m_i, d(m_i); and (2) search for an update to the current model which is a perturbation along the search direction and that reduces the objective function. The FWI processing starts from a given starting model m₀ provided by the user. FWI algorithms iteratively improve this starting model using an optimization technique,m_i+1=m_i+α_id_i,  (1)where α_i is a scalar parameter, d_i is the search direction and i is the nonlinear iteration number. The search direction is chosen along a globalization strategy [1,2]. For the second-order optimization methods, the search direction is obtained by solvingH_id_i=−g_i,  (2)where H_i can be Newton's Hessian or Gauss-Newton's Hessian. For the large scale optimization problems, the Hessian is both prohibitively large to store and compute explicitly. Instead, an approximate inverse Hessian H_i⁻¹ is used to calculate the search direction. There are several choices for this approximation, such as (i) quasi-Newton methods and (ii) truncated Newton's or Gauss-Newton methods (note that “(Gauss)-Newton” is used herein to refer to both Newton and Gauss-Newton methods).

SUMMARY

In an exemplary embodiment, a method can include: storing, in a computer memory, seismic data acquired from a seismic survey of a subsurface region; and generating, with a computer, a final subsurface physical property model of the subsurface region by processing the seismic data with an iterative full wavefield inversion method, wherein the iterative full wavefield inversion method generates the final subsurface physical property model by iteratively applying a linear solver with a preconditioner that is generated from information from one or more previous iterations of the linear solver.

In an exemplary embodiment, the linear solver is a Krylov-space method.

In an exemplary embodiment, the linear solver is a conjugate gradient method.

In an exemplary embodiment, the method can further include generating the preconditioner with a limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) method.

In an exemplary embodiment, the method can further include generating the preconditioner with a quasi-Newton method.

In an exemplary embodiment, the method can further include: storing, in a computer memory, a change in an optimization parameter of the full wavefield inversion method and a change in a gradient of a cost function used in the full wavefield inversion method for each of a plurality of iterations of the linear solver; generating the preconditioner based on the change in the optimization parameter and the change in the gradient of the cost function from each of the plurality of iterations of the linear solver; and applying the preconditioner to a subsequent iteration of the linear solver, relative to the plurality of iterations.

In an exemplary embodiment, the preconditioner is a fixed preconditioner, and the preconditioner does not change when solving a linear system and it is only based on changes in the optimization parameter and changes in the gradient of the cost function from previous iterations of the linear solver.

In an exemplary embodiment, the preconditioner is a variable preconditioner, and the preconditioner can change when solving a linear system, and is based on changes in the optimization parameter and changes in the gradient of the cost function from previous iterations of the linear solver and a current linear iteration of the linear solver.

In an exemplary embodiment, the linear solver is non-flexible.

In an exemplary embodiment, the linear solver is flexible.

In an exemplary embodiment, the change in the optimization parameter for a standard quasi-Newton algorithm is replaced with a change in a search direction of a linear system.

In an exemplary embodiment, the change in the gradient for the standard quasi-Newton algorithm is replaced with a change in a residual of the linear solver.

In an exemplary embodiment, the optimization parameter is the search direction and the gradient is a residual of a linear system.

In an exemplary embodiment, the linear solver is a Krylov-space method.

In an exemplary embodiment, the method can further include managing hydrocarbons based on the final subsurface physical property model of the subsurface region.

In an exemplary embodiment, the method can further include: creating, with a processor, an image of the subsurface region from the final subsurface physical property model.

In an exemplary embodiment, the method can further include: using the final subsurface physical property model in interpreting a subsurface region for hydrocarbon exploration or production.

In an exemplary embodiment, the method can further include drilling for hydrocarbons at a location determined using the final subsurface physical property model of the subsurface region.

In an exemplary embodiment, the linear solver is a generalized minimal residual method.

BRIEF DESCRIPTION OF THE DRAWINGS

While the present disclosure is susceptible to various modifications and alternative forms, specific example embodiments thereof have been shown in the drawings and are herein described in detail. It should be understood, however, that the description herein of specific example embodiments is not intended to limit the disclosure to the particular forms disclosed herein, but on the contrary, this disclosure is to cover all modifications and equivalents as defined by the appended claims. It should also be understood that the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating principles of exemplary embodiments of the present invention. Moreover, certain dimensions may be exaggerated to help visually convey such principles.

FIG. 1 illustrates the linear and nonlinear iterative methods that are included in the present technological advancement.

FIG. 2 illustrates an exemplary method embodying the present technological advancement.

FIG. 3 describes an exemplary application of the present technological advancement.

DESCRIPTION OF THE INVENTION

While the present disclosure is susceptible to various modifications and alternative forms, specific example embodiments thereof have been shown in the drawings and are herein described in detail. It should be understood, however, that the description herein of specific example embodiments is not intended to limit the disclosure to the particular forms disclosed herein, but on the contrary, this disclosure is to cover all modifications and equivalents as defined by the appended claims. It should also be understood that the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating principles of exemplary embodiments of the present invention. Moreover, certain dimensions may be exaggerated to help visually convey such principles.

The present technological advancement can increase the speed for convergence of FWI by several factors when second order methods are used as the optimization technique. The present technological advancement can uniquely combine two known optimization techniques: the quasi-Newton methods (such as L-BFGS) (first method) and the truncated Newton-or Gauss-Newton)-method (second method). In the present technological advancement, the second method can be used as an optimization algorithm and the first method can be used to speed up the convergence of the second method as a preconditioner.

Quasi-Newton Methods

Quasi-Newton methods replace H_i⁻¹ with its approximation in equation (2) when solving for d. These methods approximate the inverse Hessian operator (H) using gradient and model parameter changes throughout nonlinear iteration (this can be contrasted with (Gauss) Newton's method, which does not approximate the Hessian, but (approximately) solves Equation (2)).

Not all quasi-Newton methods are directly applicable to FWI problems due to its large-scale nature. However, quasi-Newton methods are modified and extended in several ways to make them suitable for large-scale optimization problems. The members of the quasi-Newton methods suitable for large scale optimization problems are so called limited-memory quasi-Newton methods. The limited-memory BFGS (Broyden-Fletcher-Goldfarb-Shanno) (L-BFGS) algorithm is the most common member suitable for FWI, as it is robust and computationally inexpensive and easy to implement [2]. All the preconditioner approaches introduced in the rest of the discussion will be based on the L-BFGS algorithm. However, the present technological advancement is not limited to the use of the L-BFGS algorithm; and other quasi-Newton algorithms such as limited memory SR1 method could also be used [2].

The inverse Hessian approximation H_i+1⁻¹ will be dense, so the cost of storing and manipulating it is computationally prohibitive for an FWI problem. To circumvent this problem, in limited memory BFGS approximation, {s_i,y_i} pairs are used to approximate the action of H_i+1⁻¹ to a given vector. The vector pairs {s_i,y_i} are defined ass_i=m_i+1−m_i, and  (3)y_i=g_i+1−g_i.  (4)

In other words, vector s_i is the change in optimization parameter and vector y_i is the change in the gradient g_i, at nonlinear iteration i.

The resulting recursive algorithm computes the application of the approximate inverse Hessian on a vector q by the L-BFGS approach using m pairs of {s_i,y_i} is given in Algorithm (1) below. Note that given vector q and m vector pairs {s_i,y_i}, the algorithm returns the vector p, which is multiplication of the approximate inverse Hessian with vector q.

Algorithm 1: Two-loop L-BFGS algorithm [2]
FOR (k = m, m − 1, ... ,1) (m is the total number of pairs and k is the
iterator for the loop)
αk=skT⁢qykT⁢sk
q = q − αkyk
END FOR
p = (Hi0)−1q
FOR (k = 1,2, ... , m)
β=ykT⁢pykT⁢sk
p = p + (αk − β)sk
END FOR
RETURN p

To complete the given L-BFGS algorithm (1), an initial estimate (H_i⁰)⁻¹ needs to be provided to the algorithm. A method for choosing (H_i⁰)⁻¹ that has proven to be effective in practice is to set

$\begin{array}{cc}{\left(H_i^0\right)}^{-1}=\frac{s_{i-1}^T{y}_{i-1}}{y_{i-1}^T{y}_{i-1}}I,& \left(5\right)\end{array}$ where the multiplier in front of identity matrix (I) is the scaling factor that attempts to estimate the size of the true Hessian matrix along the most recent search direction. The choice of scaling factor is important to ensure that the initial Hessian approximation is accurately scaled. Alternatively, the initial estimate (H_i⁰)⁻¹ can be set to an available geophysics based preconditioner (such as U.S. Patent Publication No. 2015/0073755, the entirety of which is incorporated by reference). This mechanism enables different types of preconditioners to be combined with the present technological advancement.

Truncated (Gauss)-Newton Method

Another approach for solving the (Gauss)-Newton system (2) is to use an iterative method. In contrast to the Quasi-Newton method, this approach uses the Hessian operator (or its approximation as in Gauss-Newton method) directly. However, due to difficulty of storing and explicitly computing the Hessian operator, one can only employ an iterative method to solve (2) since these methods do not explicitly require the Hessian operator in (2) but rather they require application of the Hessian operator to a given vector. A preferred approach is to use one of the so-called Krylov space methods. A Krylov space method for solving a linear system Ax=b is an iterative method starting from some initial approximation x₀ and the corresponding residual r₀=b−Ax₀, and iterating until the process possibly finds the exact solution or a stopping criteria is satisfied. These methods only require application of the Hessian operator to a given vector. For an FWI, application of Hessian on a vector can require computing at least one forward and one adjoint simulations of the wavefields. The linear iterations usually terminate using inner convergence criteria to improve speed of the nonlinear convergence of FWI[1, 2, 6]. The following discussion will utilize the conjugate gradient method as a non-limiting example of a Krylov space method.

After the search direction is computed by approximately solving (2) with an iterative linear Krylov-space solver (such as the conjugate gradient method), the FWI model is updated using a line search strategy. This procedure is repeated till convergence (see FIG. 1). The linear iterations for solving the system (2) with the conjugate gradient method are referred to as inner iterations 101, and the nonlinear iterations for updating the optimization parameter m_i are referred to as outer iterations 103, and the relationship between the two are depicted in FIG. 1.

The accompanying Appendix provides additional information regarding implementation of the conjugate gradient method.

Preconditioning

When using the truncated (Gauss)-Newton method, the performance of linear solvers used to solve (2) can be improved by using a preconditioner. In this case, instead of solving equation (2), the following system is solved with the Krylov space methodsB_i⁻¹H_id_i=−B_i⁻¹g_i,  (6)where B_i⁻¹ is the preconditioner to equation system (2) which is typically an approximation to the numerical inverse of the Hessian for optimization problems. One of the roles of the preconditioner is to reduce the condition number of the Hessian, so that equation system (6) can be solved more efficiently with less linear iterations. Note that, in (6) the preconditioner is applied from the left side of the operator H_i. There are other alternatives applications of the preconditioner, and the present technological advancement is not limited to the example provided here. (See Appendix, [2,6]).

Exemplary Embodiments

In following non-limiting exemplary embodiments of the present technological advancement, three methods are combined in a unique way: (i) Quasi-Newton method, (ii) truncated (Gauss)-Newton method and (iii) preconditioning. These exemplary embodiments use the truncated (Gauss)-Newton method as the optimization algorithm. In addition, the quasi-Newton approximation of the inverse Hessian is used as a preconditioner to the (Gauss)-Newton system. In other words, Algorithm 1 is used as a preconditioner when solving the system (2) using a Krylov-space method. To create the preconditioner of Algorithm 1, either information from the outer nonlinear iterations 103 or inner linear iterations 101 can be used to construct the preconditioner [1]. The main difference in these approaches (and one example of where the present technological advancement differentiates itself from the state-of-the art [1]) is the way the vector pairs {s_i,y_i} are created and used in the application of Algorithm 1 by the present technological advancement.

The state-of-the-art quasi-Newton preconditioning approach essentially approximates the inverse Hessian using the information captured from the outer iterations. The present technological advancement, on the contrary, introduces a quasi-Newton preconditioner which approximates the inverse Hessian using the history of inner iterations 101. The present technological advancement may include the use of any additional preconditioner to Algorithm 1 as a starting initial estimate of equation (5) for the iterative process. The present technological advancement significantly improves the convergence speed of FWI relative to the state-of-the-art preconditioning methods with an additional negligible computational cost.

First, it is observed that the solution of equation (2) via an iterative Krylov space method is equivalent to minimization of the following unconstrained quadratic optimization problem of the form

$\begin{array}{cc}\underset{d_i}{\min }f\left(d_i\right)=\frac{1}{2}{d}_i^T{H}_i{d}_i+d_i^T{g}_i..& \left(7\right)\end{array}$

It is noted that solving (2) with a Krylov subspace method, such as the conjugate gradient method, is equivalent to minimizing an objective function in the form of (7). This idea is combined with the quasi-Newton method to create a preconditioner for the truncated Newton (or truncated Gauss-Newton) method.

In solving (2) with a Krylov-subspace method, the optimization parameter is the search direction d_i and the gradient is equivalent to residual of the linear system (2). At a given inner linear iteration l, the gradient of the objective function (7) (also the residual of linear system (2)) isr^l=H_id_i^l+g_i.  (8).

In a given linear iteration of the conjugate gradient method, the solution of (2) is updated withd_i^l+1=d_i^l+μ^lγ^l.  (9).

And the residual (gradient) is updated withr^l+1=r^l+μ^lH_iγ^l,  (10)where μ^l is the step length and γ^l is the search direction for the linear system (see [2, Appendix] for details).

The present technological advancement constructs the quasi-Newton preconditioner with the information captured from this inner optimization, {s_i,y_i} pairs in equation (7). Using equation (3) and (4) along with (8) and (9) we gety^l=r^l+1−r^l=μ^lH_iγ^l, and  (11)s^l=d_i^l+1−d_i^l=μ^lγ^l.  (12)In contrast, the quasi-Newton method uses outer nonlinear iterations to create {s_l,y_l} pairs.

Noting that scaling factor μ^l in (11) and (12) can be omitted because it cancels out in the application of the preconditioner, leading tos^l=γ^l and y^l=H_iγ^l.  (13)

To generate the preconditioner, the vectors in (13) are substituted into Algorithm 1, wherein the output is the multiplication of the approximate inverse Hessian with vector q, wherein this vector q can be the residual vector r in the preconditioned conjugate gradient method (see Appendix).

These choices of vectors (13) and their use, for example, are ways in which the present technological advancement distinguishes over the conventional techniques. In other words, on using inner iterations 101 of Krylov-space method (conjugate gradient method) as a linear solver, it is recognized that: “the change in the residual of inner iterations” plays the same role as “change in the gradient of outer iterations”, and “the change in the solution of the linear system” plays the same role as “the change in the solution of the outer optimization system”.

FIG. 1 is a depiction of the outer and inner iterations that are part of the truncated (Gauss)-Newton method, and how they relate to generating an updated physical properties subsurface model. The linear or inner iterations 101_i are repeated for equations (7)-(12) until a predetermined stopping criteria is reached, at which point d_i is determined. As part of these iterations, L {s_l,y_l} pairs are generated in linear iterations 101_i, stored, and used to create the preconditioner for use in subsequent iterations. d_i, the output of the first linear iteration 101_i in FIG. 1, is used to update the model m at the outer iteration 103_i, which is repeated for the subsequent outer iterations. The outer 103 and inner 101 iterations are repeated until the updated model converges or some other predetermined stopping criteria is satisfied.

FIG. 2 is a flow chart for an exemplary method of implementing the present technological advancement. The exemplary method of FIG. 2 builds on the truncated (Gauss)-Newton method of FIG. 1 by combining it with quasi-Newton methods and preconditioning.

In step 201, in the first non-linear iteration 103i, equation (2) is solved using the conjugate gradient method, without preconditioning, by performing several linear iterations 101_i of the conjugate gradient method. Since (s_k, y_k) pairs are not available for the first iteration, the preconditioner of the present technological advancement is not used in the first nonlinear iteration. However, it is also possible to use another preconditioner when performing step 201.

In step 203, in the first nonlinear iteration 103_i, when solving equation (2), vectors γ^l, and its product with the Hessian operator H_iγ^l are stored in computer memory. Effectively, the right most terms in equations (9) and (10) can be stored for each iteration 1. As noted above, the scaling factor will cancel out and effectively it does not matter whether the scaling factor is stored along with vectors γ^l and its product with the Hessian operator H_iγ^l.

In step 205, the model is updated in outer iteration 103_i. The output of the inner iterations 101_i is d_i, which is used to update the model in the outer iteration 103_i.

In step 207, in the second nonlinear iteration 103_i+1, the stored pairs of {γ^l, H_iγ^l} from the previous nonlinear iteration 103_i are used to construct quasi-Newton “inverse Hessian mat-vec (matrix-free) operator” (Algorithm 1). This operator is used for preconditioning the conjugate gradient algorithm as given in equation (6). Only the pairs stored in the previous outer iteration(s) are used in the preconditioner application when a non-flexible conjugate gradient method is used. However, flexible version of the conjugate gradient method can also be used with the present technological advancement. The flexible version allows for a variable preconditioner within inner iterations 101. In this case, a slight modification in the application of the preconditioner is possible. Accordingly, the flexible version of the conjugate gradient method can use all pairs of {γ^l, H_iγ^l} vectors (i.e., those from the previous nonlinear iteration 103_i and those generated during the current nonlinear iteration 103_i+1). Further details of this flexible preconditioner approach are included in the Appendix.

The L-BFGS preconditioner requires a starting inverse Hessian, which is above approximated by scaled identity (5). The correct scaling can be crucial for the performance of the preconditioner. When Algorithm 1 is used as part of the present technological advancement to generate the preconditioner, the state-of-the-art preconditioner (5) constructed from the information at the outer iterations can be used as a starting inverse Hessian. Thus, the present technological advancement can combine information obtained from both the outer and inner iterations.

Provided that H_i in (2) is positive definite for a Gauss-Newton system, it can be shown that the resultant algorithm produces a positive definite operator, and guarantees a decent direction and robust algorithm. If one uses a truncated Newton's method, additional care must be taken to preserve positive definiteness [2].

In step 209, in the second linear iterations 101_i+1, the process continues to store the pairs of {γ^l, H_iγ^l} vectors generated from the current iteration (i+1) during the solving of equation (6) using the conjugate gradient method.

In step 211, the model is updated in outer iteration 103_i+1. The output of the inner iterations 101_i+1 is d_i+1, which is used to update the model in the outer iteration 103_i+1.

In step 213, it is determined whether convergence criteria or other predetermined stopping criteria has been reached for the updated physical property model. If not, the process returns to step 207 for another iteration. Iterations can continue until the convergence or stopping criteria is satisfied.

When the convergence or stopping criteria is satisfied, the process proceeds to step 215, in which a final physical property subsurface model is generated.

In step 217, the final physical property subsurface model can be used to manage hydrocarbon exploration. As used herein, hydrocarbon management includes hydrocarbon extraction, hydrocarbon production, hydrocarbon exploration, identifying potential hydrocarbon resources, identifying well locations, determining well injection and/or extraction rates, identifying reservoir connectivity, acquiring, disposing of and/or abandoning hydrocarbon resources, reviewing prior hydrocarbon management decisions, and any other hydrocarbon-related acts or activities.

The preconditioner of the present technological advancement can be used along with alternative iterative methods. The above examples use the conjugate gradient method as the linear solver. The conjugate gradient method can be replaced with another iterative method, such as the generalized minimal residual method (GMRES). To be applicable to FWI, such alternative methods need to be iterative like the conjugate gradient method.

Construction of the preconditioner of the present technological advancement is not limited to Algorithm 1. Other alternative iterative methods such as limited memory SR1 can be used. The update vectors of Algorithm 1 are the change in the residual at each iteration and the change of the iteration of the linear solver.

FIG. 3 illustrates a comparison between the present technological advancement and conventional technology. The comparison is based on the data synthetically generated using a geologically appropriate model, and uses a state-of-the-art truncated-Gauss-Newton method as the baseline algorithm. FIG. 3 displays the convergence speed-up with the present technological advancement when it is used for a sequential-source [5] FWI problem constructing-the model.

In all practical applications, the present technological advancement must be used in conjunction with a computer, programmed in accordance with the disclosures herein. Preferably, in order to efficiently perform FWI, the computer is a high performance computer (HPC), known as to those skilled in the art. Such high performance computers typically involve clusters of nodes, each node having multiple CPU's and computer memory that allow parallel computation. The models may be visualized and edited using any interactive visualization programs and associated hardware, such as monitors and projectors. The architecture of system may vary and may be composed of any number of suitable hardware structures capable of executing logical operations and displaying the output according to the present technological advancement. Those of ordinary skill in the art are aware of suitable supercomputers available from Cray or IBM.

The present techniques may be susceptible to various modifications and alternative forms, and the examples discussed above have been shown only by way of example. However, the present techniques are not intended to be limited to the particular examples disclosed herein. Indeed, the present techniques include all alternatives, modifications, and equivalents falling within the spirit and scope of the appended claims.

REFERENCES

The following references are each incorporated by the reference in their entirety: [1] L. Metivier, R. Brossier, J. Virieux, and S. Operto, “Full Waveform Inversion and Truncated Newton Method”, SIAM J. Sci. Comput., 35(2), B401-B437; [2] J. Nocedal and J. Wright, “Numerical Optimization”, 2nd Edition, Springer; [3] V. Akcelik, G. Biros, O. Ghattas, J. Hill, D. Keyes, and B. van Bloemen Waanders, “Parallel algorithms for PDE constrained optimization, in Parallel Processing for Scientific Computing”, SIAM, 2006; [4] J. L. Morales and J. Nocedal, “Automatic Preconditioning by Limited Memory Quasi-Newton Updating”, SIAM J. Optim, 10(4), 1079-1096; [5] J. R. Krebs, J. E. Anderson, D. Hinkley, R. Neelamani, S. Lee, A. Baumstein and M. D. Lacasse, “Fast Full-wavefield Seismic Inversion Using Encoded Sources”, Geophysics, 74; and [6] D. A. Knoll and D. E. Keyes, “Jacobian-free Newton-Krylov Methods: A Survey of Approaches and Applications”, SIAM J. Sci. Comp. 24:183-200, 2002.

APPENDIX

The conjugate gradient method is an algorithm for finding the numerical solution of symmteric and positive-definite systems of linear equations. The conjugate gradient method is often implemented as an iterative algorithm, applicable to sparse systems that are too large to be handled by a direct implementation or other direct methods such as the Cholesky decomposition. Large sparse systems often arise when numerically solving partial differential equations or optimization problems. The conjugate gradient method can also be used to solve unconstrained optimization problems such as least squares minimization problem.

Description of the Method

Suppose solving the following system of linear equationsAx=b,  (A1)for the vector x, where b is the known vector, and A is the known n-by-n matrix which is symmetric (i.e., A^T=A), positive definite (x^TAx>0 for all non-zero vectors x in ), and real. The unique solution of this system is given by x*.

The Conjugate Gradient Method as a Direct Method

Two non-zero vectors u and v are conjugate with respect to A ifu^TAv=0.  (A2)Since A is symmetric and positive definite, the left-hand side of equation (A2) defines an inner productu, v:=Au, v=u, A^Tv=u, Av=u^TAv,  (A3)where , is the inner product operator of two vectors. The two vectors are conjugate if and only if they are orthogonal with respect to this inner product operator. Being conjugate is a symmetric relation: if u is conjugate to v, then v is conjugate to u.

Suppose thatP={p_k: ∀i≠k, k∈[1, n], and p_i, p_k=0}  (A4)is a set of n mutually conjugate directions. Then P is a basis of , so the solution x, of Ax=b can be expanded within P space asx*=Σ_i=1ⁿα_ip_i,  (A5)which leads tob=Ax,=Σ_i=1ⁿα_iAp_i,  (A6)For any p_k ∈ P,p_k^Tb=p_k^TAx*=Σ_i=1ⁿα_ip_k^TAp_i=α_kp_k^TAp_k,  (A7)because ∀i≠ k,p_i and p_k are mutually conjugate.

$\begin{array}{cc}{\alpha }_k=\frac{p_k^{T}b}{p_k^T{\mathrm{Ap}}_k}=\frac{\langle p_k,b\rangle }{{\langle p_k,p_k\rangle }_A}.& \left(\mathrm{A8}\right)\end{array}$ This gives the following method for solving the equation Ax=b find a sequence of n conjugate directions, and then compute the coefficients α_k.

The Conjugate Gradient Method as an Iterative Method

An iterative conjugate gradient method allows to approximately solve systems of linear equations where n is so large that the direct method is computationally intractable. Suppose an initial guess for x* by x₀ and assume without loss of generality that x₀=0. Starting with x₀, while searching for the solution, in each iteration a metric to determine if the current iterate is closer to x* (that is unknown). This metric comes from the fact that the solution x* is the unique minimizer of the following quadratic function

$\begin{array}{cc}f\left(x\right)=\frac{1}{2}{x}^T\mathrm{Ax}-x^{T}b,x\in {ℝ}^n,& \left(\mathrm{A9}\right)\end{array}$ and as this function ƒ becomes smaller, solution x gets closer to x*. The search (descent) direction for function ƒ in (A9) equals to the negative gradient b−Ax. Starting from a guessed solution x₀ (x₀=0 in case of no guessed solution) at the k^th step, this descent direction isr_k=b−Ax_k  (A10)The conjugation constraint described previously is an orthonormal-type constraint and hence the algorithm bears resemblance to Gram-Schmidt orthonormalization. This gives the following expression for conjugate of r_k

$\begin{array}{cc}{p}_k=r_k-\sum _{i<k}\frac{p_i^T{\mathrm{Ar}}_k}{p_i^T{\mathrm{Ap}}_i}{p}_i& \left(\mathrm{A11}\right)\end{array}$ Following this direction, the next optimal location is

$\begin{array}{cc}{x}_{k+1}=x_k+{\alpha }_k{p}_k\mathrm{where},& \left(\mathrm{A12}\right)\\ {\alpha }_k=\frac{p_k^{T}b}{p_k^T{\mathrm{Ap}}_k}=\frac{p_k^T{r}_{k-1}}{p_k^T{\mathrm{Ap}}_k},& \left(A13\right)\end{array}$ where the last equality holds because p_k and x_k−1 are conjugate.

Conjugate Gradient Algorithm

The above algorithm gives the straightforward explanation of the conjugate gradient method. Seemingly, the algorithm as stated requires storage of all previous searching directions and residue vectors, as well as many matrix-vector multiplications, and thus can be computationally expensive. However, a closer analysis of the algorithm shows that r_k+1 is conjugate to p_i for all i<k, and therefore only r_k, p_k, and x_k are needed to construct r_k+1, p_k+1, and x_k+1. Furthermore, only one matrix-vector multiplication is needed in each iteration.

A modified algorithm is detailed below for solving Ax=b where A is a real, symmetric, positive-definite matrix, with an input vector x₀ (a guessed solution otherwise 0).

Algorithm A1: A conjugate gradient algorithm.
r0 = b − Ax0
p0 = r0
k = 0
REPEAT
αk=rkT⁢rkpkT⁢Apk
xk+1 = xk + αkpk
rk+1 = rk − αkApk
IF rk+1 is sufficiently small, EXIT REPEAT
βk=rk+1T⁢rk+1rkT⁢rk
pk+1 = rk+1 + βkpk
k = k + 1
END REPEAT
RETURN xk+1

Preconditioned Conjugate Gradient Method

Preconditioning speeds up convergence of the conjugate gradient method. A preconditioned conjugate gradient algorithm is given in Algorithm A2, which requires an application of preconditioner operator B⁻¹ on a given vector in addition to the steps in Algorithm A1.

Algorithm A2: A preconditioned conjugate gradient algorithm.
r0 = b − Ax0
z0 = B−1r0
p0 = z0
k = 0
REPEAT
αk=rkT⁢rkpkT⁢Apk
xk+1 = xk + αkpk
rk+1 = rk − αkApk
IF rk+1 is sufficiently small, EXIT REPEAT
zk+1 = B−1rk+1
βk=rk+1⁢zk+1Trk⁢zkT
pk+1 = rk+1 + βkpk
k = k +1
END REPEAT
RETURN xk+1

The preconditioner matrix B has to be symmetric positive-definite and fixed, i.e., cannot change from iteration to iteration. If any of these assumptions on the preconditioner is violated, the behavior of Algorithm A2 becomes unpredictable and its convergence can not be guaranteed.

Flexible Preconditioned Conjugate Gradient Method

For some numerically challenging applications, Algorithm A2 can be modified to accept variable preconditioners, changing between iterations, in order to improve the convergence performance of Algorithm A2. For instance, the Polak-Ribière formula

$\begin{array}{cc}{\beta }_k=\frac{z_{k+1}^T\left(r_{k+1}-r_k\right)}{z_k^T{r}_k},& \left(\mathrm{A14}\right)\end{array}$ instead of the Fletcher-Reeves formula used in Algorithm A2,

$\begin{array}{cc}{\beta }_k=\frac{z_{k+1}^T{r}_{k+1}}{z_k^T{r}_k},& \left(\mathrm{A15}\right)\end{array}$ may dramatically improve the convergence of the preconditioned conjugate-gradient method. This version of the preconditioned conjugate gradient method can be called flexible, as it allows for variable preconditioning. The implementation of the flexible version requires storing an extra vector. For a fixed preconditioner, z_k+1^Tt_k=0 both Polak-Ribière and Fletcher-Reeves formulas are equivalent. The mathematical explanation of the better convergence behavior of the method with the Polak-Ribière formula is that the method is locally optimal in this case, in particular, it does not converge slower than the locally optimal steepest descent method.

Claims

1. A method, comprising: storing, in a computer memory, seismic data acquired from a seismic survey of a subsurface region;generating, with a computer, a final subsurface physical property model of the subsurface region by processing the seismic data with an iterative full wavefield inversion method, wherein the iterative full wavefield inversion method includes a non-linear outer iteration process comprising a plurality of outer iterations, each of which updates physical property values m of the physical property model and a nested linear inner iteration process, wherein: (1) each iteration i of the plurality of non-linear outer iterations includes a nested linear inner iteration process for determining a search direction di to be used in updating the physical property values mi of said non-linear outer iteration; and(2) each of the second and subsequent outer iterations of the plurality of non-linear outer iterations further includes applying a preconditioner to the nested linear inner iteration process of such respective second and subsequent non-linear outer iteration, wherein the preconditioner at such non-linear outer iteration i+1 is generated based at least in part upon the nested linear inner iteration process of the immediately previous outer iteration i and wherein the preconditioner is generated using vectors {γl, Hiγl} solved for in the nested inner iteration process, where γl is the search direction for the linear system and where Hiγl is the product of the vector γl with a Hessian operator Hi; andcreating, with a process, an image of the subsurface region from the final subsurface physical property model.

2. The method of claim 1, comprising using a Krylov-space method as a linear solver for the nested inner iteration processes.

3. The method of claim 1, comprising using a conjugate gradient method as a linear solver for the nested inner iteration processes.

4. The method of claim 2, wherein each preconditioner is generated based on the nested inner iteration process of the immediately previous outer iteration using a limited-memory Broyden-Fletcher-Goldfarb-Shanno (BFGS) method.

5. The method of claim 2, wherein each preconditioner is generated based on the nested inner iteration process of the immediately previous outer iteration using a quasi-Newton method.

6. The method of claim 1, wherein the preconditioner of each of the second and subsequent outer iterations is a variable preconditioner such that the preconditioner can change when solving a linear system, and the preconditioner is generated based on the nested iteration process of the then-current outer iteration in addition to the nested inner iteration process of the immediately previous outer iteration.

7. The method of claim 1, further comprising managing hydrocarbons based on the final subsurface physical property model of the subsurface region.

8. The method of claim 1, further comprising: using the final subsurface physical property model in interpreting a subsurface region for hydrocarbon exploration or production.

9. The method of claim 1, further comprising drilling for hydrocarbons at a location determined using the final subsurface physical property model of the subsurface region.

10. The method of claim 1, wherein a generalized minimal residual method is used as a linear solver for the nested inner iteration processes.