Full Waveform Inversion Problem Research Articles

Efficient 3D elastic full-waveform inversion using wavefield reconstruction methods

In reverse time migration (RTM) or full-waveform inversion (FWI), forward and reverse time propagating wavefields are crosscorrelated in time to form either the image condition in RTM or the misfit gradient in FWI. The crosscorrelation condition requires both fields to be available at the same time instants. For large-scale 3D problems, it is not possible, in practice, to store snapshots of the wavefields during forward modeling due to extreme storage requirements. We have developed an approximate wavefield reconstruction method that uses particle velocity field recordings on the boundaries to reconstruct the forward wavefields during the computation of the reverse time wavefields. The method is computationally effective and requires less storage than similar methods. We have compared the reconstruction method to a boundary reconstruction method that uses particle velocity and stress fields at the boundaries and to the optimal checkpointing method. We have tested the methods on a 2D vertical transversely isotropic model and a large-scale 3D elastic FWI problem. Our results revealed that there are small differences in the results for the three methods.

A parallel evolution strategy for an earth imaging problem in geophysics

In this paper we propose a new way to compute a rough approximation solution, to be later used as a warm starting point in a more refined optimization process, for a challenging global optimization problem related to earth imaging in geophysics. The warm start consists of a velocity model that approximately solves a full-waveform inverse problem at low frequency. Our motivation arises from the availability of massively parallel computing platforms and the natural parallelization of evolution strategies as global optimization methods for continuous variables. Our first contribution consists of developing a new and efficient parametrization of the velocity models to significantly reduce the dimension of the original optimization space. Our second contribution is to adapt a class of evolution strategies to the specificity of the physical problem at hands where the objective function evaluation is known to be the most expensive computational part. A third contribution is the development of a parallel evolution strategy solver, taking advantage of a recently proposed modification of these class of evolutionary methods that ensures convergence and promotes better performance under moderate budgets. The numerical results presented demonstrate the effectiveness of the algorithm on a realistic 3D full-waveform inverse problem in geophysics. The developed numerical approach allows us to successfully solve an acoustic full-waveform inversion problem at low frequencies on a reasonable number of cores of a distributed memory computer.

Acoustic 2D full waveform inversion to solve gas cloud challenges

<p>The existing conventional inversion algorithm does not provide satisfactory results due to the complexity of propagated wavefield though the gas cloud. Acoustic full waveform inversion has been developed and applied to a realistic synthetic offshore shallow gas cloud feature with Student-t approach, with and without simultaneous sources encoding. As a modeling operator, we implemented the grid based finite-difference method in frequency domain using second order elastic wave equation. Jacobin operator and its adjoint provide a necessary platform for solving full waveform inversion problem in a reduced Hessian matrix. We invert gas cloud model in 5 frequency band selected from 1 to 12 Hz, each band contains 3 frequencies. The inversion results are highly sensitive to the misfit. The model allows better convergence and recovery of amplitude losses. This approach gives better resolution then the existing least-squares approach. In this paper, we implement the full waveform inversion for low frequency model with minimum number of iteration providing a better resolution of inversion results.</p>

Enhancement of near-surface elastic full waveform inversion results in regions of low sensitivities

Elastic full waveform inversion of high-resolution seismic data is a potentially very powerful option for imaging the shallow subsurface. Unfortunately, the success of traditional full waveform inversion applied to such problems is limited by a very uneven sensitivity distribution, which can be attributed to the uneven amplitudes of body and surface waves. As a result, very shallow structures are well resolved by fitting the large amplitude surface waves, but the imaging quality decreases rapidly with depth. To account for uneven sensitivity distributions, we present a novel scaling approach that enhances weak sensitivities in regions of interest. To this end, the column sums of the Jacobian matrix – each of them corresponding to one model parameter – are equalized prior to updating the model. The performance of this methodology is demonstrated by inverting two synthetic, but realistic, data sets. Both the P- and S-wave velocity images were improved significantly by applying the new scaling technique. Our results are particularly relevant for shallow elastic full waveform inversion problems, but we also see benefits of our technique for other surface-based geophysical methods, such as geoelectrics or electromagnetics.

An adjoint-based Jacobi-type iterative method for elastic full waveform inversion problem

Full waveform inversion (FWI) is a promising technique that is capable of creating high-resolution subsurface images of the earth from seismic data. However, it is computationally expensive, especially when 3D elastic wave equation is considered. In this paper an adjoint-based computational algorithm has been proposed to address the computational challenges. The important strategy in this work is to decouple the two Lamé parameters so that two half-sized subproblems are resulted and solved separately. Mathematically, the inverse problem is formulated as a PDE-constrained optimization problem in which the objective functional is defined as the misfit between observational and synthetic data. The parameters to be recovered are the spatially varying Lamé parameters which are of great interests to geophysicists for the purpose of hydrocarbon exploration and subsurface imaging. The gradient of the misfit functional with respect to the Lamé parameters is calculated by solving the adjoint elastic wave equation, whilst the Quasi-Newton method (L-BFGS) is used to minimize the misfit functional. Numerical experiments demonstrated that the new method is accurate, efficient and robust in recovering Lamé parameters from seismic data.

A frozen Gaussian approximation-based multi-level particle swarm optimization for seismic inversion

In this paper, we propose a frozen Gaussian approximation (FGA)-based multi-level particle swarm optimization (MLPSO) method for seismic inversion of high-frequency wave data. The method addresses two challenges in it: First, the optimization problem is highly non-convex, which makes hard for gradient-based methods to reach global minima. This is tackled by MLPSO which can escape from undesired local minima. Second, the character of high-frequency of seismic waves requires a large number of grid points in direct computational methods, and thus renders an extremely high computational demand on the simulation of each sample in MLPSO. We overcome this difficulty by three steps: First, we use FGA to compute high-frequency wave propagation based on asymptotic analysis on phase plane; Then we design a constrained full waveform inversion problem to prevent the optimization search getting into regions of velocity where FGA is not accurate; Last, we solve the constrained optimization problem by MLPSO that employs FGA solvers with different fidelity. The performance of the proposed method is demonstrated by a two-dimensional full-waveform inversion example of the smoothed Marmousi model.

An adaptive subspace trust-region method for frequency-domain seismic full waveform inversion

Full waveform inversion is currently considered as a promising seismic imaging method to obtain high-resolution and quantitative images of the subsurface. It is a nonlinear ill-posed inverse problem, the main difficulty of which that prevents the full waveform inversion from widespread applying to real data is the sensitivity to incorrect initial models and noisy data. Local optimization theories including Newton's method and gradient method always lead the convergence to local minima, while global optimization algorithms such as simulated annealing are computationally costly. To confront this issue, in this paper we investigate the possibility of applying the trust-region method to the full waveform inversion problem. Different from line search methods, trust-region methods force the new trial step within a certain neighborhood of the current iterate point. Theoretically, the trust-region methods are reliable and robust, and they have very strong convergence properties. The capability of this inversion technique is tested with the synthetic Marmousi velocity model and the SEG/EAGE Salt model. Numerical examples demonstrate that the adaptive subspace trust-region method can provide solutions closer to the global minima compared to the conventional Approximate Hessian approach and the L-BFGS method with a higher convergence rate. In addition, the match between the inverted model and the true model is still excellent even when the initial model deviates far from the true model. Inversion results with noisy data also exhibit the remarkable capability of the adaptive subspace trust-region method for low signal-to-noise data inversions. Promising numerical results suggest this adaptive subspace trust-region method is suitable for full waveform inversion, as it has stronger convergence and higher convergence rate.

A computational method for full waveform inversion of crosswell seismic data using automatic differentiation

Full waveform inversion (FWI) is a model-based data-fitting technique that has been widely used to estimate model parameters in Geophysics. In this work, we propose an efficient computational approach to solve the FWI of crosswell seismic data. The FWI problem is mathematically formulated as a partial differential equation (PDE)-constrained optimization problem, which is numerically solved using a gradient-based optimization method. The efficiency and accuracy of FWI are mainly determined by the three main components: forward modeling, gradient calculation and model update which usually involves the gradient-based optimization algorithm. Given the large number of iterations needed by FWI, an accurate gradient is critical for the success of FWI, as it will not only speed up the convergence but also increase the accuracy of the solution. However computing the gradient still remains a challenging task even after the adjoint PDE has been derived. Automatic differentiation (AD) tools have been proved very effective in a variety of application areas including Geoscience. In this work we investigated the feasibility of integrating TAPENADE, a powerful AD tool into FWI, so that the FWI workflow is simplified to allow us to focus on the forward modeling and the model updating. In this paper we choose the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) method due to its robustness and fast convergence. Numerical experiments have been conducted to demonstrate the effectiveness, efficiency and robustness of the new computational approach for FWI.

Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures

ABSTRACTFull waveform inversion is a powerful tool for quantitative seismic imaging from wide‐azimuth seismic data. The method is based on the minimization of the misfit between observed and simulated data. This amounts to the solution of a large‐scale nonlinear minimization problem. The inverse Hessian operator plays a crucial role in this reconstruction process. Accounting accurately for the effect of this operator within the minimization scheme should correct for illumination deficits, restore the amplitude of the subsurface parameters, and help to remove artefacts generated by energetic multiple reflections. Conventional minimization methods (nonlinear conjugate gradient, quasi‐Newton methods) only roughly approximate the effect of this operator. In this study, we are interested in the truncated Newton minimization method. These methods are based on the computation of the model update through a matrix‐free conjugate gradient solution of the Newton linear system. We present a feasible implementation of this method for the full waveform inversion problem, based on a second‐order adjoint state formulation for the computation of Hessian‐vector products. We compare this method with conventional methods within the context of 2D acoustic frequency full waveform inversion for the reconstruction of P‐wave velocity models. Two test cases are investigated. The first is the synthetic BP 2004 model, representative of the Gulf of Mexico geology with high velocity contrasts associated with the presence of salt structures. The second is a 2D real data‐set from the Valhall oil field in North sea. Although, from a computational cost point of view, the truncated Newton method appears to be more expensive than conventional optimization algorithms, the results emphasize its increased robustness. A better reconstruction of the P‐wave velocity model is provided when energetic multiple reflections make it difficult to interpret the seismic data. A better trade‐off between regularization and resolution is obtained when noise contamination of the data requires one to regularize the solution of the inverse problem.

The variable projection method for waveform inversion with an unknown source function

This paper compares three alternative algorithms for simultaneously estimating a source wavelet at the same time as an earth model in full-waveform inversion: (i) simultaneous descent, (ii) alternating descent and (iii) descent with the variable projection method. The latter is a technique for solving separable least-squares problems that is well-known in the applied mathematics literature. When applied to full-waveform inversion, it involves making the source wavelet an implicit function of the earth model via a least-squares filter-estimation process. Since the source wavelet becomes purely a function of medium parameters, it no longer needs to be treated as a separate unknown in the inversion. Essentially, the predicted data are projected onto the measured data in a least-squares sense at every function evaluation, making use of the fact that the filter estimation problem is trivial when compared to the full-waveform inversion problem. Numerical tests on a simple 1D model indicate that the variable projection method gives the best result; actually producing results in quality that are very similar to control experiments with a known, correct wavelet.

Full Waveform Inversion and the Truncated Newton Method

Full waveform inversion (FWI) is a powerful method for reconstructing subsurface parameters from local measurements of the seismic wavefield. This method consists in minimizing the distance between predicted and recorded data. The predicted data are computed as the solution of a wave-propagation problem. Conventional numerical methods for the resolution of FWI problems are gradient-based methods, such as the preconditioned steepest descent, or more recently the $l$-BFGS quasi-Newton algorithm. In this study, we investigate the desirability of applying a truncated Newton method to FWI. The inverse Hessian operator plays a crucial role in the parameter reconstruction. The truncated Newton method allows one to better account for this operator. This method is based on the computation of the Newton descent direction by solving the corresponding linear system through an iterative procedure such as the conjugate gradient method. The large-scale nature of FWI problems requires us, however, to carefully implement this method to avoid prohibitive computational costs. First, this requires working in a matrix-free formalism and the capability of computing efficiently Hessian-vector products. To this purpose, we propose general second-order adjoint state formulas. Second, special attention must be paid to define the stopping criterion for the inner linear iterations associated with the computation of the Newton descent direction. We propose several possibilities and establish a theoretical link between the Steihaug--Toint method, based on trust regions, and the Eisenstat and Walker stopping criterion, designed for a method globalized by linesearch. We investigate the application of the truncated Newton method to two test cases: The first is a standard test case in seismic imaging based on the Marmousi II model. The second one is inspired by a near-surface imaging problem for the reconstruction of high-velocity structures. In the latter case, we demonstrate that the presence of large amplitude multiscattered waves prevents standard methods from converging while the truncated Newton method provides more reliable results.

Quasi-Newton full-waveform inversion with a projected Hessian matrix

We present a method, in realistic-size full-waveform inversion (FWI), to explicitly construct a projected Hessian matrix and its inverse matrix, which we subsequently used to solve FWI with a quasi-Newton method. Newton’s method is practically unfeasible in solving realistic-size FWI problems because of the prohibitive cost (computing time and memory consumption) of calculating the Hessian matrix and the inverse Hessian. Therefore, the Gauss-Newton method and various quasi-Newton methods are proposed to approximate a Hessian matrix. Particularly, current quasi-Newton FWI (QNFWI) commonly uses the limited-memory BFGS (L-BFGS) method, which, however, only implicitly approximates an inverse Hessian. We repose FWI as a sparse optimization problem in a sparse model space, which contains substantially fewer model parameters that are constrained by structures of the model. With respect to fewer parameters in the sparse model, we can avoid the “limited-memory” approximation and are able to explicitly compute and store a projected Hessian matrix that saves the computational time and required memory. We constructed such a projected Hessian matrix by adapting the classic BFGS method to a projected BFGS (P-BFGS) method in the sparse space. Using the projected Hessian matrix and its inverse, we can apply the P-BFGS method to solve FWI with a quasi-Newton method. In QNFWI with P-BFGS because we invert for a sparse model with much fewer parameters, the memory required to compute the projected Hessian is negligible compared to either forward modeling or gradient calculation. QNFWI with P-BFGS converges in fewer iterations than conjugate-gradient based methods and QNFWI with L-BFGS.

Compressed implicit Jacobian scheme for elastic full-waveform inversion

SUMMARY We present a regularized Gauss–Newton (GN) inversion method for solving the elastic fullwaveform inversion problem in the frequency domain. The main bottleneck of this method is the Jacobian matrix storage and the computational cost of calculating the GN step (the inner-loop calculation). In this work, we managed to reduce the memory usage by calculating the Jacobian matrix on the fly in each inner-loop iteration. By doing so the computational cost of calculating the GN step increases; however, this overhead is mitigated by compressing the fieldmatricesusingtheadaptivecrossapproximationscheme.Forsomecases,thiscompressed implicit Jacobian scheme may even speed-up the GN step calculation and further regularizes the GN method. As examples, we present inversion results of cross-well seismic and surface seismic data.

Monte Carlo full-waveform inversion of crosshole GPR data using multiple-point geostatistical a priori information

We present a general Monte Carlo full-waveform inversion strategy that integrates a priori information described by geostatistical algorithms with Bayesian inverse problem theory. The extended Metropolis algorithm can be used to sample the a posteriori probability density of highly nonlinear inverse problems, such as full-waveform inversion. Sequential Gibbs sampling is a method that allows efficient sampling of a priori probability densities described by geostatistical algorithms based on either two-point (e.g., Gaussian) or multiple-point statistics. We outline the theoretical framework for a full-waveform inversion strategy that integrates the extended Metropolis algorithm with sequential Gibbs sampling such that arbitrary complex geostatistically defined a priori information can be included. At the same time we show how temporally and/or spatiallycorrelated data uncertainties can be taken into account during the inversion. The suggested inversion strategy is tested on synthetic tomographic crosshole ground-penetrating radar full-waveform data using multiple-point-based a priori information. This is, to our knowledge, the first example of obtaining a posteriori realizations of a full-waveform inverse problem. Benefits of the proposed methodology compared with deterministic inversion approaches include: (1) The a posteriori model variability reflects the states of information provided by the data uncertainties and a priori information, which provides a means of obtaining resolution analysis. (2) Based on a posteriori realizations, complicated statistical questions can be answered, such as the probability of connectivity across a layer. (3) Complex a priori information can be included through geostatistical algorithms. These benefits, however, require more computing resources than traditional methods do. Moreover, an adequate knowledge of data uncertainties and a priori information is required to obtain meaningful uncertainty estimates. The latter may be a key challenge when considering field experiments, which will not be addressed here.

Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data

We have applied the finite-difference contrast-source inversion (FDCSI) method to seismic full-waveform inversion problems. The FDCSI method is an iterative nonlinear inversion algorithm. However, unlike the nonlinear conjugate gradient method and the Gauss-Newton method, FDCSI does not solve any full forward problem explicitly in each iterative step of the inversion process. This feature makes the method very efficient in solving large-scale computational problems. It is shown that FDCSI, with a significant lower computation cost, can produce inversion results comparable in quality to those produced by the Gauss-Newton method and better than those produced by the nonlinear conjugate gradient method. Another attractive feature of the FDCSI method is that it is capable of employing an inhomogeneous background medium without any extra or special effort. This feature is useful when dealing with time-lapse inversion problems where the objective is to reconstruct changes between the baseline and the monitor model. By using the baseline model as the background medium in crosswell seismic monitoring problems, high quality time-lapse inversion results are obtained.

Full Waveform Inversion for Seismic Velocity and Anelastic Losses in Heterogeneous Structures

We present a least-squares optimization method for solving the nonlinear full waveform inverse problem of determining the crustal velocity and intrinsic at- tenuation properties of sedimentary valleys in earthquake-prone regions. Given a known earthquake source and a set of seismograms generated by the source, the in- verse problem is to reconstruct the anelastic properties of a heterogeneous medium with possibly discontinuous wave velocities. The inverse problem is formulated as a constrained optimization problem, where the constraints are the partial and ordinary differential equations governing the anelastic wave propagation from the source to the receivers in the time domain. This leads to a variational formulation in terms of the material model plus the state variables and their adjoints. We employ a wave propaga- tion model in which the intrinsic energy-dissipating nature of the soil medium is mod- eled by a set of standard linear solids. The least-squares optimization approach to inverse wave propagation presents the well-known difficulties of ill posedness and multiple minima. To overcome ill posedness, we include a total variation regulariza- tion functional in the objective function, which annihilates highly oscillatory material property components while preserving discontinuities in the medium. To treat multi- ple minima, we use a multilevel algorithm that solves a sequence of subproblems on increasingly finer grids with increasingly higher frequency source components to re- main within the basin of attraction of the global minimum. We illustrate the metho- dology with high-resolution inversions for two-dimensional sedimentary models of the San Fernando Valley, under SH-wave excitation. We perform inversions for both the seismic velocity and the intrinsic attenuation using synthetic waveforms at the observer locations as pseudoobserved data.