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 the 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 (iteratively refined WRI [IR-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 is 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 indicate the resilience of IR-WRI to cycle skipping and noise, as well as its ability to reconstruct with high-fidelity, large-contrast salt bodies and subsalt structures starting the inversion from crude initial models and a 3 Hz starting frequency.

Full-waveform inversion (FWI) is an iterative nonlinear waveform matching procedure, which seeks to reconstruct unknown model parameters from partial waveform measurements. The nonlinear and ill-posed nature of FWI requires sophisticated regularization techniques to solve it. In most applications, the model parameters may be described by physical properties (e.g., wave speeds, density, attenuation, and anisotropy) that are piecewise functions of space. Compound regularizations are, thus, beneficial to capture these different functions 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 mitigates the so-called cycle skipping pathology. Our hybrid regularizations rely on the 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 parameters are 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, and piecewise linear). We implement these compound regularizations in WRI using the alternating direction method of multipliers (ADMM). Then, we assess our regularized WRI for seismic imaging applications. Using a wide range of subsurface models, we conclude that the 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.

The search space of full-waveform inversion (FWI) can be extended via a relaxation of the wave equation to increase the linear regime of the inversion. This wave-equation relaxation is implemented by solving jointly (in a least-squares sense) the wave equation weighted by a penalty parameter and the observation equation such that the reconstructed wavefields closely match the data, hence preventing cycle skipping at receivers. Then, the subsurface parameters are updated by minimizing the temporal and spatial source extension generated by the wave-equation relaxation to push back the data-assimilated wavefields toward the physics. This extended formulation of FWI has been efficiently implemented in the frequency domain with the augmented Lagrangian method in which the overdetermined systems of the data-assimilated wavefields can be solved separately for each frequency with linear-algebra methods and the sensitivity of the optimization to the penalty parameter is mitigated through the action of the Lagrange multipliers. Applying this method in the time domain is, however, hampered by two main issues: the computation of data-assimilated wavefields with explicit time-stepping schemes and the storage of the Lagrange multipliers capturing the history of the source residuals in the state space. These two issues are solved by recognizing that the source residuals on the right side of the extended wave equation, when formulated in a form suitable for explicit time stepping, are related to the extended data residuals through an adjoint equation. This relationship first allows us to relate the extended data residuals to the reduced data residuals through a normal equation in the data space. Once the extended data residuals have been estimated by solving (exactly or approximately) this normal equation, the data-assimilated wavefields are computed with explicit time-stepping schemes by cascading an adjoint and a forward simulation. Second, the relationship between the source and the data residuals allows us to project the Lagrange multipliers tied to the wave-equation constraint in the data space at the expense of additional adjoint simulations. This leads to a formulation of the extended FWI that is very similar to that of the classic reduced-space FWI. The computational cost of the extended FWI is two times that of the classic reduced-space FWI, and the storage requirement is the same. Applications to different synthetic benchmarks involving large contrasts confirm the resilience of the method to cycle skipping even when the quadratic problem for the extended data residuals is solved approximately with one inner iteration of a steepest-descent method.

Full-waveform inversion (FWI) can be made immune to cycle skipping by matching the recorded data with traveltime errors smaller than one-half period 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 aimed at jointly satisfying the wave equation and fitting the data in a least-squares sense. This leads to data-assimilated wavefields that are computed by solving the wave equation in the inaccurate background model with a data-dependent source extension added to the source term. Then, the subsurface parameters are updated by canceling out these additional source terms, sometimes inaccurately called wave equation errors, to push the background model toward the true model in the left-side wave equation operator. Although many studies are devoted to these approaches with promising numerical results, their governing physical principles and their relationships with classical FWI do not seem to be understood well yet. The goal of this tutorial is to review these principles in the framework of inverse scattering theory whose governing forward equation is the Lippmann-Schwinger equation. From this equation, we find 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 and in the full Newton Hessian 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 assess their potential and limits and design relevant heuristics to manage the latter.

SUMMARY: Regularization is necessary for solving non-linear 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), smooth or non-smooth, variational-based or not. Nevertheless, tailoring a suitable and easy-to-implement prior for describing geophysical models is a non-trivial task. In this paper, we propose two generic optimization algorithms to implement arbitrary regularization in non-linear inverse problems such as full-waveform inversion (FWI), where the regularization task is recast as a denoising problem. We assess these optimization algorithms with the plug-and-play block matching (BM3D) regularization algorithm, which determines empirical priors adaptively without any optimization formulation. The non-linear inverse problem is solved with 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 non-linear objective function. The second relies on the alternating direction method of multipliers (ADMM), and hence Newton-ADMM (NADMM), where the least-squares optimization subproblem and the regularization subproblem in the composite objective function are decoupled through auxiliary variable and solved in an alternating mode. The least-squares subproblem can be solved with exact, inexact, or quasi-Newton methods. We compare NISTA and NADMM numerically by solving FWI with BM3D regularization. The tests show promising results obtained by both algorithms. However, NADMM shows a faster convergence rate than NISTA when using L-BFGS to solve the Newton system.

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 nonsmooth and describe structurally complex media. Thus, we solve the acoustic and visco-acoustic scattered-field (Lippmann-Schwinger) wave equation in the frequency domain with PINN instead of the wave equation to remove the source singularity. We first illustrate that nonsmooth 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 to better couple the real and imaginary parts of the wavefield. Moreover, we design new neurons by replacing the classical affine function with a quadratic function in the argument of the activation function to better capture the nonsmooth features of the wavefields. We find that the PML condition and the quadratic functions improve the results including handling attenuation and discuss the reason for this improvement. We also illustrate that a network trained to predict a wavefield for a specific medium can be used as an initial model of the neural network for predicting other wavefields corresponding to PDE-coefficient alterations and improving the convergence speed accordingly. This pretraining strategy should find applications in iterative full-waveform inversion and time-lag target-oriented imaging when the model perturbation between two consecutive iterations or two consecutive experiments is small.