Search

Full-waveform inversion (FWI) with extended sources first computes wavefields with data-driven source extensions, such that the simulated data in inaccurate velocity models match the observed counterpart well enough to prevent cycle skipping. Then, the source extensions are minimized to update the model parameters. This two-step workflow is iterated until both data and sources are matched. It was recently shown that the source extensions are the least-squares solutions of the recorded scattered data fitting problem. As a result, they are computed by propagating backward in time the deblurred FWI data residuals, where the deblurring operator is the inverse of the damped data-domain Hessian of the scattering-source estimation problem. Estimating the deblurred data residuals is the main computational bottleneck of time-domain extended-source FWI (ES-FWI). To mitigate this issue, we first estimate them when the inverse of the data-domain Hessians is approximated by matching filters in Fourier and short-time Fourier domains. Second, we refine them with conjugate-gradient iterations when necessary. Computing the matching filters and performing one conjugate-gradient iteration each require two simulations per source. Therefore, it is critical to design some workflows that minimize this computational burden. We implement time-domain ES-FWI with the augmented Lagrangian method. Moreover, we further extend its linear regime with a multiscale frequency continuation approach, which is combined with grid coarsening to mitigate the computational burden and regularize the inversion. Finally, we use total-variation regularization to deal with large-contrast reconstruction. We present synthetic cases where different inversion workflows carried out with data-domain Hessians of variable accuracy were assessed with the aim at converging toward accurate solutions while minimizing computational cost.

Full Waveform Inversion can be made immune to cycle skipping by matching the recorded data arbitrarily well from inaccurate subsurface models. To achieve this goal, the simulated wavefields can be computed in an extended search space as the solution of an overdetermined problem aiming at jointly satisfying the wave equation and fitting the data in a least-squares sense. Simply put, the wavefields are computed by solving the wave equation in the inaccurate background model with a feedback term to the data added to the physical source in the right-hand side. Then, the subsurface parameters are updated by canceling out these additional source terms, sometimes called unwisely wave-equation errors, to push the background model toward the true model in the left-hand side wave-equation operator. Although many studies were devoted to these approaches with promising numerical results, their governing physical principles and their relationships with classical FWI don't seem to be understood well yet. The goal of this tutorial is to review these principles in the theoretical framework of inverse scattering theory whose governing forward equation is the Lippmann-Schwinger equation. From this equation, we show how the data-assimilated wavefields embed an approximation of the scattered field generated by the sought model perturbation and how they modify the sensitivity kernel of classical FWI beyond the Born approximation. We also clarify how the approximation with which these wavefields approximate the unknown true wavefields is accounted for in the adjoint source of the parameter estimation problem. The theory is finally illustrated with numerical examples. Understanding the physical principles governing these methods is a necessary prerequisite to assessing their potential and limits and designing relevant heuristics to manage the latter.

Frequency-domain simulation of seismic waves plays an important role in seismic inversion, but it remains challenging in large models. The recently proposed physics-informed neural network (PINN), as an effective deep learning method, has achieved successful applications in solving a wide range of partial differential equations (PDEs), and there is still room for improvement on this front. For example, PINN can lead to inaccurate solutions when PDE coefficients are non-smooth and describe structurally-complex media. In this paper, we solve the acoustic and visco-acoustic scattered-field wave equation in the frequency domain with PINN instead of the wave equation to remove source singularity. We first illustrate that non-smooth velocity models lead to inaccurate wavefields when no boundary conditions are implemented in the loss function. Then, we add the perfectly matched layer (PML) conditions in the loss function of PINN and design a quadratic neural network to overcome the detrimental effects of non-smooth models in PINN. We show that PML and quadratic neurons improve the results as well as attenuation and discuss the reason for this improvement. We also illustrate that a network trained during a wavefield simulation can be used to pre-train the neural network of another wavefield simulation after PDE-coefficient alteration and improve the convergence speed accordingly. This pre-training strategy should find application in iterative full waveform inversion (FWI) and time-lag target-oriented imaging when the model perturbation between two consecutive iterations or two consecutive experiments can be small.

Implementation of the standard full waveform inversion (FWI) poses difficulties as the initial model offsets from the true model. The wavefield reconstruction inversion (WRI) was proposed to mitigate these difficulties by relaxing the wave-equation constraint. In this abstract, working on the nonlinear term in the Hessian matrix of FWI, we develop a new approximate Hessian as an Augmented Gauss-Newton (AGN) Hessian including second-order derivative information. Moreover, we establish an intimate connection between an updating formula which results from approximate solve of the Newton's method with the AGN Hessian on the FWI problem and the WRI method. Our analysis opens new perspectives for developing efficient algorithms for FWI based on the Newton's method and highlights the importance of the nonlinear term in the Hessian matrix, which is ignored in most cases.

Full-waveform inversion (FWI) is a high-resolution and computationally intensive imaging technique to reconstruct unknown parameters in the computational model in which the waves propagate; however, an accurate model of only part of this medium is required for some applications. To decrease the computational burden of such problems, target-oriented FWI was proposed where the redatumed data on the part of the medium or localized solvers for the wave equation are used. On the other hand, the classical formulation of FWI suffers from non-linearity and ill-posedness, which makes FWI sensitive to the initial model, the low-frequency content of the data, and limited illumination. In this study, we propose a localized version of the alternating direction method of multipliers (ADMM)-based FWI method, which was proposed to solve these problems in classical FWI. In our localized FWI or LWI, the medium is decomposed into a few subdomains, where some of them are updated, and the others are kept fixed based on an adaptation of multi-block ADMM, which is a powerful algorithm for solving inverse problems with decomposition and block separability. Numerical tests on the Marmousi model for a time-lapse application confirm the computational efficiency and robustness against background velocity model errors.

Full-waveform inversion (FWI) is a seismic imaging method that provides quantitative inference about subsurface properties with a wavelength-scale resolution. Its frequency-domain formulation is computationally efficient when processing only a few discrete frequencies. However, classical FWI, which is formulated on the reduced-parameter space, requires starting the inversion with a sufficiently-accurate initial model and low frequency to prevent being stuck in local minima due to cycle skipping. FWI with extended search space has been proposed to mitigate this issue. It contains two main steps: first, data-assimilated (DA) wavefields are computed by allowing for wave-equation errors to match the data at receivers closely. Then, subsurface parameters are estimated from these wavefields by minimizing the wave-equation errors. The DA wavefields are the least-squares solution of an overdetermined system gathering the wave and observation equations. The numerical bandwidth of the resulting normal-equation system is two times that of the wave-equation system, which can be a limiting factor for 3D large-scale applications. Therefore, computing highly accurate DA wavefields at a reasonable computational cost is an issue in extended FWI. This issue is addressed here by rewriting the normal system such that its solution can be computed by solving the time-harmonic wave equation several times in sequence. Moreover, the computational burden of multi-right-hand side (RHS) simulations is mitigated with a sketching method. Finally, we solve the time-harmonic wave equation with the convergent Born series method, which conciliates accuracy and computational efficiency. Application of the new extended FWI algorithm on the salt benchmark shows that it reconstructs at a reasonable cost subsurface models that are similar to those obtained with the classical extended FWI.

The full-waveform inversion (FWI) addresses the computation and characterization of subsurface model parameters by matching predicted data to observed seismograms in the frame of nonlinear optimization. We formulate FWI as a nonlinearly constrained optimization problem, for which a regularization term is minimized subject to the nonlinear data matching constraint. Unlike FWI which is based on the penalty function, the method of multipliers solves the resulting optimization problems by using the augmented Lagrangian function; and leads to a two-step recursive algorithm. The primal step requires solving an unconstrained minimization problem like the traditional FWI with a difference that the data are replaced by the Lagrange multipliers. The dual step involves an update of the Lagrange multipliers. The overall performance of the algorithm is improved considering that this multiplier method does not require an exact solution of these primal-dual subproblems. In fact, convergence is attained when only one step of a gradient-based method is taken on both subproblems. The proposed algorithm greatly improves the overall performance of FWI such as convergence from inaccurate starting models and robustness with respect to the determination of the step length. Furthermore, it can be performed by the existing FWI engines with minimal change. We only have to replace the observed data at each iteration with the multipliers, thus all the nice properties of the traditional FWI algorithms are kept. Numerical experiments confirm that the multipliers waveform inversion can converge to a solution of the inverse problem in the absence of low-frequency data from an inaccurate initial model even with a constant step size.

Efficient frequency-domain Full Waveform Inversion (FWI) of long-offset/wide-azimuth node data can be designed with a few discrete frequencies. However, 3D frequency-domain seismic modeling remains challenging since it requires solving a large and sparse linear indefinite system per frequency. When such systems are solved with direct methods or hybrid direct/iterative solvers, based upon domain decomposition preconditioner, finite-difference stencils on regular Cartesian grids should be designed to conciliate compactness and accuracy, the former being necessary to mitigate the fill-in induced by the Lower-Upper (LU) factorization. Compactness is classically implemented by combining several second-order accurate stencils covering the eight cells surrounding the collocation point, leading to the so-called 27-point stencil. Accuracy is obtained by applying optimal weights on the different stiffness and consistent mass matrices such that numerical dispersion is jointly minimized for several number of grid points per wavelength ($G$). However, with this approach, the same weights are used at each collocation point, leading to suboptimal accuracy in heterogeneous media. In this study, we propose a straightforward recipe to improve the accuracy of the 27-point stencil. First, we finely tabulate the values of $G$ covering the range of wavelengths spanned by the subsurface model and the frequency. Then, we estimate with a classical dispersion analysis in homogeneous media the corresponding table of optimal weights that minimize dispersion for each $G$ treated separately. We however apply a Tikhonov regularization to guarantee smooth variation of the weights with $G$. Finally, we build the impedance matrix by selecting the optimal weights at each collocation point according to the local wavelength, hence leading to a wavelength-adaptive stencil.

Partial differential equation (PDE) constrained optimization problems such as seismic full waveform inversion (FWI) frequently arise in the geoscience and related fields. For such problems, many observations are usually gathered by multiple sources, which form the right-hand-sides of the PDE constraint. Solving the inverse problem with such massive data sets is computationally demanding, in particular when dealing with large number of model parameters. This paper proposes a novel randomized source sketching method for the efficient resolution of multisource PDE constrained optimization problems. We first formulate the different source-encoding strategies used in seismic imaging into a unified framework based on a randomized sketching. To this end, the source dimension of the problem is projected in a smaller domain by a suitably defined projection matrix that gathers the physical sources in super-sources through a weighted summation. This reduction in the number of physical sources decreases significantly the number of PDE solves while suitable sparsity-promoting regularization can efficiently mitigate the footprint of the cross-talk noise to maintain the convergence speed of the algorithm. We implement the randomized sketching method in an extended search-space formulation of frequency-domain FWI, which relies on the alternating direction method of multipliers (ADMM). Numerical examples carried out with a series of well documented 2D benchmarks demonstrate that the randomized sketching algorithm reduces the cost of large-scale problems by at least one order of magnitude compared to the original deterministic algorithm.

The augmented Lagrangian (AL) method provides a flexible and efficient framework for solving extended-space full-waveform inversion (FWI), a constrained nonlinear optimization problem whereby we seek model parameters and wavefields that minimize the data residuals and satisfy the wave equation constraint. The AL-based wavefield reconstruction inversion, also known as iteratively refined wavefield reconstruction inversion, extends the search space of FWI in the source dimension and decreases sensitivity of the inversion to the initial model accuracy. Furthermore, it benefits from the advantages of the alternating direction method of multipliers (ADMM), such as generality and decomposability for dealing with non-differentiable regularizers, e.g., total variation regularization, and large scale problems, respectively. In practice any extension of the method aiming at improving its convergence and decreasing the number of wave-equation solves would have a great importance. To achieve this goal, we recast the method as a general fixed-point iteration problem, which enables us to apply sophisticated acceleration strategies like Anderson acceleration. The accelerated algorithm stores a predefined number of previous iterates and uses their linear combination together with the current iteration to predict the next iteration. We investigate the performance of the proposed accelerated algorithm on a simple checkerboard model and the benchmark Marmousi II and 2004 BP salt models through numerical examples. These numerical results confirm the effectiveness of the proposed algorithm in terms of convergence rate and the quality of the final estimated model.

Full waveform inversion (FWI) is beginning to be used to characterize weak seismic events at different scales, an example of which is microseismic event (MSE) characterization. However, FWI with unknown sources is a severely underdetermined optimization problem, and hence requires strong prior information about the sources and/or the velocity model. The frequency-domain wavefield reconstruction inversion method (WRI) has shown promising results to mitigate the nonlinearity of the FWI objective function that is generated by cycle-skipping. WRI relies on the reconstruction of data-assimilated wavefields, which approach the true wavefields near the receivers, a helpful feature when the source is added as an additional optimization variable. We present an adaptation of a recently proposed version of WRI based on the alternating direction method of multipliers (ADMM) that first finds the location of the MSEs and then reconstructs the wavefields and the source signatures jointly. Finally, the subsurface model is updated to focus the MSEs at their true location. The method does not require prior knowledge of the number of MSEs. The inversion is stabilized by sparsifying regularizations separately tailored to the source location and velocity model subproblems. The method is tested on the Marmousi model using one MSE and two clusters of MSEs with two different initial velocity models, an accurate one and a rough one, as well as with added noise. In all cases, the method accurately locates the MSEs and recovers their source signatures.

The emergence of long-offset sparse stationary-recording surveys carried out with ocean bottom nodes (OBN) makes frequency-domain full waveform inversion (FWI) attractive to manage compact volume of data and perform attenuation imaging. One challenge of frequency-domain FWI is the forward problem, which requires the solution of large and sparse linear systems with multiple right-hand sides. While direct methods are suitable for dense acquisitions and problems involving less than 100 million unknowns, iterative solver are more suitable for large computational domains covered by sparse OBN surveys. Here, we solve these linear systems with a Krylov subspace method preconditioned with the two-level Optimized Restricted Additive Schwarz (ORAS) domain decomposition preconditioner, the prefix optimized referring to the use of absorbing conditions at the subdomain interfaces. We implement this method with finite differences on uniform grid and finite elements on unstructured tetrahedral meshes. A simulation in a model where the velocity linearly increases with depth allows us to validate the accuracy of the two schemes against an analytical solution while highlighting how their relative cost varies with the band of propagated wavelengths. A simulation in the overthrust model involving up to 2 billions of parameters allows us to tune the method and highlights its scalability., Comment: arXiv admin note: substantial text overlap with arXiv:2004.07930, arXiv:2004.06309

Wave propagation problems are notoriously difficult to solve. Time-harmonic problems are especially challenging in mid and high frequency regimes. The main reason is the oscillatory nature of solutions, meaning that the number of degrees of freedom after discretisation increases drastically with the wave number, giving rise to large complex-valued problems to solve. Additional difficulties occur when the problem is defined in a highly heterogeneous medium, as is often the case in realistic physical applications. For time-discretised problems of Maxwell type, the main challenge remains the significant kernel in curl-conforming spaces, an issue that impacts on the design of robust preconditioners. This has already been addressed theoretically for a homogeneous medium but not yet in the presence of heterogeneities. In this review we provide a big-picture view of the main difficulties encountered when solving wave propagation problems, from the first step of their discretisation through to their parallel solution using two-level methods, by showing their limitations on a few realistic examples. We also propose a new preconditioner inspired by the idea of subspace decomposition, but based on spectral coarse spaces, for curl-conforming discretisations of Maxwell's equations in heterogeneous media.

Full waveform inversion (FWI) requires an accurate estimation of source signatures. Due to the coupling between the source signatures and the subsurface model, small errors in the former can translate into large errors in the latter. When direct methods are used to solve the forward problem, classical frequency-domain FWI efficiently processes multiple sources for source signature and wavefield estimations once a single Lower-Upper (LU) decomposition of the wave-equation operator has been performed. However, this efficient FWI formulation is based on the exact solution of the wave equation and hence is highly sensitive to the inaccuracy of the velocity model due to the cycle skipping pathology. Recent extended-space FWI variants tackle this sensitivity issue through a relaxation of the wave equation combined with data assimilation, allowing the wavefields to closely match the data from the first inversion iteration. Then, the subsurface parameters are updated by minimizing the wave-equation violations. When the wavefields and the source signatures are jointly estimated with this approach, the extended wave equation operator becomes source dependent, hence making direct methods ineffective. In this paper, we propose a simple method to bypass this issue and estimate source signatures efficiently during extended FWI. The proposed method replaces each source with a blended source during each data-assimilated wavefield reconstruction to make the extended wave equation operator source independent. Besides computational efficiency, the additional degrees of freedom introduced by spatially distributing the sources allows for a better signature estimation at the physical location when the velocity model is rough. Numerical tests on the Marmousi II and 2004 BP salt synthetic models confirm the efficiency and the robustness of the proposed method.

Extended full-waveform inversion (FWI) has shown promising results for accurate estimation of subsurface parameters when the initial models are not sufficiently accurate. Frequency-domain applications have shown that the augmented Lagrangian (AL) method solves the inverse problem accurately with a minimal effect of the penalty parameter choice. Applying this method in the time domain, however, is limited by two main factors: (1) The challenge of data-assimilated wavefield reconstruction due to the lack of an explicit time-stepping and (2) The need to store the Lagrange multipliers, which is not feasible for the field-scale problems. We show that these wavefields are efficiently determined from the associated data (projection of the wavefields onto the receivers space) by using explicit time stepping. Accordingly, based on the augmented Lagrangian, a new algorithm is proposed which performs in "data space" (a lower dimensional subspace of the full space) in which the wavefield reconstruction step is replaced by reconstruction of the associated data, thus requiring optimization in a lower dimensional space (convenient for handling the Lagrange multipliers). We show that this new algorithm can be implemented efficiently in the time domain with existing solvers for the FWI and at a cost comparable to that of the FWI while benefiting from the robustness of the extended FWI formulation. The results obtained by numerical examples show high-performance of the proposed method for large scale time-domain FWI.

Full waveform inversion (FWI) is a nonlinear PDE constrained optimization problem, which seeks to estimate constitutive parameters of a medium such as phase velocity, density, and anisotropy, by fitting waveforms. Attenuation is an additional parameter that needs to be taken into account in viscous media to exploit the full potential of FWI. Attenuation is more easily implemented in the frequency domain by using complex-valued velocities in the time-harmonic wave equation. These complex velocities are frequency-dependent to guarantee causality and account for dispersion. Since estimating a complex frequency-dependent velocity at each grid point in space is not realistic, the optimization is generally performed in the real domain by processing the phase velocity (or slowness) at a reference frequency and attenuation (or quality factor) as separate real parameters. This real parametrization requires an a priori empirical relation (such as the nonlinear Kolsky-Futterman (KF) or standard linear solid (SLS) attenuation models) between the complex velocity and the two real quantities, which is prone to generate modeling errors if it does not represent accurately the attenuation behavior of the subsurface. Moreover, it leads to a multivariate inverse problem, which is twice larger than the actual size of the medium and ill-posed due to the cross-talk between the two classes of real parameters. To alleviate these issues, we present a mono-variate algorithm that solves directly the optimization problem in the complex domain by processing in sequence narrow bands of frequencies under the assumption of band-wise frequency dependence of the sought complex velocities.

Frequency-domain full-waveform inversion (FWI) is suitable for long-offset stationary-recording acquisition, since reliable subsurface models can be reconstructed with a few frequencies and attenuation is easily implemented without computational overhead. In the frequency domain, wave modeling is a Helmholtz-type boundary-value problem which requires to solve a large and sparse system of linear equations per frequency with multiple right-hand sides (sources). This system can be solved with direct or iterative methods. While the former are suitable for FWI application on 3D dense OBC acquisitions covering spatial domains of moderate size, the later should be the approach of choice for sparse node acquisitions covering large domains (more than 50 millions of unknowns). Fast convergence of iterative solvers for Helmholtz problems remains however challenging in high frequency regime due to the non definiteness of the Helmholtz operator, on one side and on the discretization constraints in order to minimize the dispersion error for a given frequency, on the other side, hence requiring efficient preconditioners. In this study, we use the Krylov subspace GMRES iterative solver combined with a two-level domain-decomposition preconditioner. Discretization relies on continuous Lagrange finite elements of order 3 on unstructured tetrahedral meshes to comply with complex geometries and adapt the size of the elements to the local wavelength ($h$-adaptivity). We assess the accuracy, the convergence and the scalability of our method with the acoustic 3D SEG/EAGE Overthrust model up to a frequency of 20~Hz., Comment: arXiv admin note: substantial text overlap with arXiv:2004.06309

Frequency-domain full-waveform inversion (FWI) is suitable for long-offset stationary-recording acquisition, since reliable subsurface models can be reconstructed with a few frequencies and attenuation is easily implemented without computational overhead. In the frequency domain, wave modelling is a Helmholtz-type boundary-value problem which requires to solve a large and sparse system of linear equations per frequency with multiple right-hand sides (sources). This system can be solved with direct or iterative methods. While the former are suitable for FWI application on 3D dense OBC acquisitions covering spatial domains of moderate size, the later should be the approach of choice for sparse node acquisitions covering large domains (more than 50 millions of unknowns). Fast convergence of iterative solvers for Helmholtz problems remains however challenging due to the non definiteness of the Helmholtz operator, hence requiring efficient preconditioners. In this study, we use the Krylov subspace GMRES iterative solver combined with a multi-level domain-decomposition preconditioner. Discretization relies on continuous finite elements on unstructured tetrahedral meshes to comply with complex geometries and adapt the size of the elements to the local wavelength ($h$-adaptivity). We assess the convergence and the scalability of our method with the acoustic 3D SEG/EAGE Overthrust model up to a frequency of 20~Hz and discuss its efficiency for multi right-hand side processing.

Regularization is necessary for solving nonlinear ill-posed inverse problems arising in different fields of geosciences. The base of a suitable regularization is the prior expressed by the regularizer, which can be non-adaptive or adaptive (data-driven). In this paper, we propose general black-box regularization algorithms for solving nonlinear inverse problems such as full-waveform inversion (FWI), which admit empirical priors that are determined adaptively by sophisticated denoising algorithms. The nonlinear inverse problem is solved by a proximal Newton method, which generalizes the traditional Newton step in such a way to involve the gradients/subgradients of a (possibly non-differentiable) regularization function through operator splitting and proximal mappings. Furthermore, it requires to account for the Hessian matrix in the regularized least-squares optimization problem. We propose two different splitting algorithms for this task. In the first, we compute the Newton search direction with an iterative method based upon the first-order generalized iterative shrinkage-thresholding algorithm (ISTA), and hence Newton-ISTA (NISTA). The iterations require only Hessian-vector products to compute the gradient step of the quadratic approximation of the nonlinear objective function. The second relies on the alternating direction method of multipliers (ADMM), and hence Newton-ADMM (NADMM), where the least-square optimization subproblem and the regularization subproblem in the composite are decoupled through auxiliary variable and solved in an alternating mode. We compare NISTA and NADMM numerically by solving full-waveform inversion with BM3D regularizations. The tests show promising results obtained by both algorithms. However, NADMM shows a faster convergence rate than Newton-ISTA when using L-BFGS to solve the Newton system.

Wavefield reconstruction inversion (WRI) extends the search space of Full Waveform Inversion (FWI) by allowing for wave equation errors during wavefield reconstruction to match the data from the first iteration. Then, the wavespeeds are updated from the wavefields by minimizing the source residuals. Performing these two tasks in alternating mode breaks down the nonlinear FWI as a sequence of two linear subproblems, relaying on the bilinearity of the wave equation. We solve this biconvex optimization with the alternating-direction method of multipliers (ADMM) to cancel out efficiently the data and source residuals in iterations and stabilize the parameter estimation with appropriate regularizations. Here, we extend WRI to viscoacoustic media for attenuation imaging. Attenuation reconstruction is challenging because of the small imprint of attenuation in the data and the cross-talks with velocities. To address these issues, we recast the multivariate viscoacoustic WRI as a triconvex optimization and update wavefields, squared slowness, and attenuation factor in alternating mode at each WRI iteration. This requires to linearize the attenuation-estimation subproblem via an approximated trilinear viscoacoustic wave equation. The iterative defect correction embedded in ADMM corrects the errors generated by this linearization, while the operator splitting allows us to tailor $\ell{1}$ regularization to each parameter class. A toy numerical example shows that these strategies mitigate cross-talk artifacts and noise from the attenuation reconstruction. A more realistic synthetic example representative of the North Sea validates the method.

Extended formulation of Full Waveform Inversion (FWI), called Wavefield Reconstruction Inversion (WRI), offers potential benefits of decreasing the nonlinearity of the inverse problem by replacing the explicit inverse of the ill-conditioned wave-equation operator of classical FWI (the oscillating Green functions) with a suitably defined data-driven regularized inverse. This regularization relaxes the wave-equation constraint to reconstruct wavefields that match the data, hence mitigating the risk of cycle skipping. The subsurface model parameters are then updated in a direction that reduces these constraint violations. However, in the case of a rough initial model, the phase errors in the reconstructed wavefields may trap the waveform inversion in a local minimum leading to inaccurate subsurface models. In this paper, in order to avoid matching such incorrect phase information during the early WRI iterations, we design a new cost function based upon phase retrieval, namely a process which seeks to reconstruct a signal from the amplitude of linear measurements. This new formulation, called Wavefield Inversion with Phase Retrieval (WIPR), further improves the robustness of the parameter estimation subproblem by a suitable phase correction. We implement the resulting WIPR problem with an alternating-direction approach, which combines the Majorization-Minimization (MM) algorithm to linearise the phase-retrieval term and a variable splitting technique based upon the alternating direction method of multipliers (ADMM). This new workflow equipped with Tikhonov-total variation (TT) regularization, which is the combination of second-order Tikhonov and total variation regularizations and bound constraints, successfully reconstructs the 2004 BP salt model from a sparse fixed-spread acquisition using a 3~Hz starting frequency and a homogeneous initial velocity model.

Full waveform inversion (FWI) is a nonlinear waveform matching procedure, which suffers from cycle skipping when the initial model is not kinematically-accurate enough. To mitigate cycle skipping, wavefield reconstruction inversion (WRI) extends the inversion search space by computing wavefields with a relaxation of the wave equation in order to fit the data from the first iteration. Then, the subsurface parameters are updated by minimizing the source residuals the relaxation generated. Capitalizing on the wave-equation bilinearity, performing wavefield reconstruction and parameter estimation in alternating mode decomposes WRI into two linear subproblems, which can solved efficiently with the alternating-direction method of multiplier (ADMM), leading to the so-called IR-WRI. Moreover, ADMM provides a suitable framework to implement bound constraints and different types of regularizations and their mixture in IR-WRI. Here, IR-WRI is extended to multiparameter reconstruction for VTI acoustic media. To achieve this goal, we first propose different forms of bilinear VTI acoustic wave equation. We develop more specifically IR-WRI for the one that relies on a parametrisation involving vertical wavespeed and Thomsen's parameters delta and epsilon. With a toy numerical example, we first show that the radiation patterns of the virtual sources generate similar wavenumber filtering and parameter cross-talks in classical FWI and IR-WRI. Bound constraints and TV regularization in IR-WRI fully remove these undesired effects for an idealized piecewise constant target. We show with a more realistic long-offset case study representative of the North Sea that anisotropic IR-WRI successfully reconstruct the vertical wavespeed starting from a laterally homogeneous model and update the long-wavelengths of the starting epsilon model, while a smooth delta model is used as a passive background model.

The nonlinear and ill-posed nature of full waveform inversion (FWI) requires us to use sophisticated regularization techniques to solve it. In most applications, the model parameters may be described by physical properties (e.g., wave speeds, density, attenuation, anisotropic parameters) which are piecewise functions of space. Compound regularizations are thus necessary to reconstruct properly such parameters by FWI. We consider different implementations of compound regularizations in the wavefield reconstruction inversion (WRI) method, a formulation of FWI that extends its search space and prevent the so-called cycle skipping pathology. Our hybrid regularizations rely on Tikhonov and total variation (TV) functionals, from which we build two classes of hybrid regularizers: the first class is simply obtained by a convex combination (CC) of the two functionals, while the second relies on their infimal convolution (IC). In the former class, the model of parameters is required to simultaneously satisfy different priors, while in the latter the model is broken into its basic components, each satisfying a distinct prior (e.g. smooth, piecewise constant, piecewise linear). We implement these types of compound regularizations in the WRI optimization problem using the alternating direction method of multipliers (ADMM). Then, we assess our regularized WRI in the framework of seismic imaging applications. Using a wide range of subsurface models, we conclude that compound regularizer based on IC leads to the lowest error in the parameter reconstruction compared to that obtained with the CC counterpart and the Tikhonov and TV regularizers when used independently.

Full waveform inversion (FWI) is a waveform matching procedure, which can provide a subsurface model with a wavelength-scale resolution. However, this high resolution makes FWI prone to cycle skipping, which drives the inversion to a local minimum when the initial model is not accurate enough. Other sources of nonlinearities and ill-posedness are noise, uneven illumination, approximate wave physics and parameter cross-talks. All these sources of error require robust and versatile regularized optimization approaches to mitigate their imprint on FWI while preserving its intrinsic resolution power. To achieve this goal, we implement bound constraints and total variation (TV) regularization in the so-called frequency-domain wavefield-reconstruction inversion (WRI) with the alternating direction method of multipliers (ADMM). In the ADMM framework, WRI relies on an augmented Lagrangian function, a combination of penalty and Lagrangian functions, to extend the FWI search space by relaxing the wave-equation constraint during early iterations. Moreover, ADMM breaks down the joint wavefield reconstruction plus parameter estimation problem into a sequence of two linear subproblems, whose solutions are coordinated to provide the solution of the global problem. The decomposability of ADMM is further exploited to interface in a straightforward way bound constraints and TV regularization with WRI via variable splitting and proximal operators. The resilience of our regularized WRI formulation to cycle skipping and noise as well as its resolution power are illustrated with two targets of the large-contrast BP salt model. Starting from a 3Hz frequency and a crude initial model, the extended search space allows for the reconstruction of the salt and subsalt structures with a high fidelity.

Full waveform inversion (FWI) is an iterative nonlinear waveform matching procedure subject to wave-equation constraint. FWI is highly nonlinear when the wave-equation constraint is enforced at each iteration. To mitigate nonlinearity, wavefield-reconstruction inversion (WRI) expands the search space by relaxing the wave-equation constraint with a penalty method. The pitfall of this approach resides in the tuning of the penalty parameter because increasing values should be used to foster data fitting during early iterations while progressively enforcing the wave-equation constraint during late iterations. However, large values of penalty parameter lead to ill-conditioned problems. Here, this tuning issue is solved by replacing the penalty method by an augmented Lagrangian method equipped with operator splitting (IR-WRI as iteratively-refined WRI). It is shown that IR-WRI is similar to a penalty method in which data and sources are updated at each iteration by the running sum of the data and source residuals of previous iterations. Moreover, the alternating direction strategy exploits the bilinearity of the wave equation constraint to linearize the subsurface model estimation around the reconstructed wavefield. Accordingly, the original nonlinear FWI is decomposed into a sequence of two linear subproblems, the optimization variable of one subproblem being passed as a passive variable for the next subproblem. The convergence of WRI and IR-WRI are first compared with a simple transmission experiment, which lies in the linear regime of FWI. Under the same conditions, IR-WRI converges to a more accurate minimizer with a smaller number of iterations than WRI. More realistic case studies performed with the Marmousi II and the BP salt models show the resilience of IR-WRI to cycle skipping and noise.