<figure><img src="image_1.png" /></figure>

HAL
open science

# Full Waveform Inversion beyond the Born approximation: A tutorial review

S. Operto

<figure><img src="image_2.png" /></figure>

## To cite this version:

S. Operto. Full Waveform Inversion beyond the Born approximation: A tutorial review. 2023. [hal-04010372](https://hal.science/hal-04010372)

HAL Id: hal-04010372

https://hal.science/hal-04010372

Preprint submitted on 1 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# FULL WAVEFORM INVERSION BEYOND THE BORN APPROXIMATION: A TUTORIAL REVIEW

A PREPRINT

## Stéphane Operto

University Cote d'Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. operto@geoazur.unice.fr

## Ali Gholami

Institute of Geophysics, University of Tehran, Tehran, Iran. agholami@ut.ac.ir

## Hossein S. Aghamiry

University Cote d'Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. aghamiry@geoazur.unice.fr

## Gaoshan Guo

University Cote d'Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. guo@geoazur.unice.fr

## Frichnel Mamfoumbi

University Cote d'Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. mamfoumbi@geoazur.unice.fr

## Stephen Beller

University Cote d'Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. beller@geoazur.unice.fr

December 20, 2022

## ABSTRACT

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, sometime called unwisely wave-equation errors, to push the background model towards 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 assess their potential and limits and design relevant heuristics to manage the latter.
# 1 Introduction

Full Waveform Inversion (FWI) (Gauthier et al., 1986; Mora, 1987; Pratt et al., 1998; Tarantola, 1984) has become the baseline seismic imaging method to build high-resolution multi-parameter subsurface models from full-azimuth long-offset data since a decade (Sirgue et al., 2010). Despite its popularity, FWI remains a highly nonlinear problem due to the oscillatory nature of seismic signals, which makes the classical least-squares data-fitting misfit function highly multimodal. In this context, a central issue of FWI is cycle skipping, which traps FWI in spurious minima as soon as the simulated data don't predict the recorded counterpart with an error lower than half a period (Virieux & Operto, 2009). Remembering that FWI is an inverse scattering problem solved with gradient-based methods, this stringent condition results from the weak-scattering Born approximation with which the nonlinear data fitting problem is linearized around a background model to recast it as a sequence of surrogate linear problems. On the acquisition side, sources with low-frequency content ($\ge$ 1.5 Hz) have been designed to mitigate cycle skipping (e.g. Baeten et al., 2013; Brenders et al., 2022). On the imaging side, new FWI formulations match the data arbitrarily well from inaccurate velocity models, hence preventing cycle skipping. To achieve this goal, a first optimization subproblem is solved to build a more convex distance in place of the usual least-squares data misfit distance before minimizing this new distance for parameter estimation. Several approaches can be viewed to match the data with inaccurate models such as optimal transport (e.g., Engquist et al., 2016; Métivier et al., 2018; Yang et al., 2018), trace-by-trace matching filters (e.g., Warner & Guasch, 2016), dynamic warping based approaches (e.g., Ma & Hale, 2013) and error in constitutive equation approaches (e.g., Banerjee et al., 2013). The later extend the search space of classical FWI by adding some degrees of freedom in the forward-problem equation to match the data with the desired accuracy. These approaches include the so-called Extended-Born modeling approaches (e.g Barnier et al., 2012; Biondi & Almomin, 2014; Symes, 2008), the contrast-source method (e.g., Abubakar et al., 2009; van den Berg & Kleinman, 1997), the wavefield reconstruction inversion method (e.g., van Leeuwen & Herrmann, 2016, 2013) and the extended-source FWI (e.g., Huang et al., 2018a; Wang et al., 2016). Extended-Born modeling approaches extends the search space in the model domain with subsurface offsets or time lags, while the last three methods, which mainly differ in the parametrization of the optimization variables, extend the search space in the source domain.

<figure><img src="image_3.png" /><figcaption>Figure 1: Sketch of new formulations of FWI aiming at minimizing cycle skipping. A first optimization subproblem is solved to build a more convex distance than the least-squares data misfit function by matching the observables well enough. This can be performed by optimal transport, dynamic warping, matching filters, and approximated wave equation relaxation. In the later case, the extended search space can be implemented in the model domain (Extended-Born modeling methods) or in the source domain (contrast-source, wavefield reconstruction, extended-source methods). These approaches modify the adjoint source of the parameter-estimation subproblem. Some approaches also modify the sensitivity kernel of classical FWI.</figcaption></figure>

All these approaches modify the source of the adjoint-state equation in the descent direction of the model updating and hence can be recast as a generalized FWI (Tarantola & Valette, 1982), where a data-domain covariance matrix in the least-squares data-fitting misfit function transforms the data residuals such that the weighted misfit function is more convex. Some of them also modify the sensitivity kernel of FWI when the covariance matrix depends on the model
parameters as we will show in this study.

This study focuses on the error in constitutive equation approaches where the search space extension is implemented in the source domain, which will be called extended-space (ES) FWI. After the seminal paper of van Leeuwen & Herrmann (2013), ES-FWI have been investigated in several studies (Table 1) with however specific angles of approach of mathematical, algorithmic or applied natures but none of them provided a comprehensive overview of the physical principles governing these methods and their underlying approximations. The goal of this tutorial is to fill this gap by gathering the lessons learnt from these studies into a self-contained review of these principles. By clarifying these physical principles, we will clearly draw the similarities and differences between classical FWI and ES-FWI. That is, we will show that ES-FWI relies on the same diffraction-based physical principles as those governing classical FWI, while mitigating the detrimental effects of the Born approximation with more accurate approximation of the wavefields. Many inverse scattering problems revolve around the Lippmann-Schwinger equation giving the scattered field by a model perturbation. Accordingly, this equation will be the cornerstone of this tutorial. Let's add that this preserved connection between ES-FWI and the physics of diffraction somehow contrasts with other approaches that are more inspired by signal or image processing like matching filters and dynamic warping based approaches or by optimization theory like optimal transport. Finally, we would like to summarize the governing idea of ES-FWI by borrowing two sentences from the preface of Chavent & Sabatier (1996) introducing the approximate inverse scattering method of Fiddy & Pommet (1996). In a similar manner to ES-FWI, this method aims at improving the ability of inverse scattering theory based on the Born approximation to reconstruct strong scattering objects from the Lippmann-Schwinger equation: The method is better than Born one but it keeps the physical features that appear in this approximate method. It goes also as close as possible to exact methods without losing physical features.

In the theory section, we first review the bilinearity of the wave equation to remind that ES-FWI aims at reconstructing wavefields that best approximate the true ones. Second, we review the Lippmann-Schwinger equation and the Born approximation before reminding their role in classical FWI. Then, we explain how improved wavefields can be reconstructed from the Lippmann-Schwinger equation and the recorded scattered data before explaining how they modify the sensitivity kernel of FWI. Finally, we explain how the approximation with which these wavefields are reconstructed is accounted for in the adjoint source of ES-FWI. We briefly review some important numerical aspects of ES-FWI to design workable algorithms before illustrating the theory with two numerical examples.

# 2 Theory

## 2.1 On the bilinearity of the wave equation

To introduce the governing idea of ES-FWI, we may begin with reminding the bilinearity of the wave equation. In this study, we will restrict ourselves to the Helmholtz equation written in matrix form as:

$$A(\mathbf{m})\mathbf{u} = \mathbf{b}, \qquad (2.1)$$

where $A(\mathbf{m}) = \omega^2\mathrm{diag}(\mathbf{m}) + \nabla^2$ is the so-called impedance matrix, $\omega$ is the angular frequency, $\mathrm{diag}(\bullet)$ is a diagonal matrix formed with the coefficients of vector $\bullet$, model parameters $\mathbf{m}$ are squared slowness, $\nabla^2$ is the Laplace operator, $\mathbf{u}$ is the monochromatic wavefield and $\mathbf{b}$ is the source. Moreover, $\bullet_i$ denotes the $i^{th}$ entry of the vector $\bullet$ through the rest of the tutorial.

Equation 2.1 is linear in $\mathbf{u}$ for known $\mathbf{m}$. Permuting $\mathbf{m}$ and $\mathbf{u}$ in the mass term, $\mathrm{diag}(\mathbf{m})\mathbf{u} = \mathrm{diag}(\mathbf{u})\mathbf{m}$, recasts equation 2.1 as

$$L(\mathbf{u})\mathbf{m} = \mathbf{y}(\mathbf{u}), \qquad (2.2)$$

where $L(\mathbf{u}) = \omega^2\mathrm{diag}(\mathbf{u})$ and $\mathbf{y}(\mathbf{u}) = \mathbf{b} - \nabla^2\mathbf{u}$. Equation 2.2 is affine in $\mathbf{m}$ for known $\mathbf{u}$ and $\mathbf{m}$ can be inferred from $\mathbf{u}$ by a pointwise division since $L(\mathbf{u})$ is a full-rank diagonal matrix:

$$\mathbf{m} = \frac{\mathbf{y}(\mathbf{u})}{L(\mathbf{u})}. \qquad (2.3)$$

This implies that, if we could record a monochromatic wavefield triggered by a monochromatic source everywhere in the subsurface, then we could reconstruct $\mathbf{m}$ exactly (Figure 2). This condition is indeed never satisfied since data are typically recorded near the surface, which makes FWI an underdetermined problem. Accordingly, fitting the data at receivers with an error smaller than half the period is a necessary condition to avoid cycle skipping but it is not sufficient to guarantee convergence toward the global minimizer since the wavefield should be matched not only at receivers but also everywhere in the targeted domain.
<table><caption>Table 1: Main references about ES-FWI since the first publications on the contrasted source method.</caption><tr><td>Chavent & Sabatier (1996)</td><td>Lecture notes on Inverse problems of wave propagation and diffraction.</td></tr><tr><td>van den Berg & Kleinman (1997)</td><td>Contrast source inversion method for strongly scattering objects.</td></tr><tr><td>Abubakar et al. (2008, 2009, 2011)</td><td>Application of contrast-source method to seismic waveform inversion.</td></tr><tr><td>van Leeuwen & Herrmann (2013)</td><td>Wavefield Reconstruction Inversion (WRI) with alternating directions.</td></tr><tr><td>van Leeuwen et al. (2014)</td><td>Tentative analogy between source extension and adjoint wavefield.</td></tr><tr><td>van Leeuwen & Herrmann (2016)</td><td>Reduced form of WRI (variable projection) and mathematical aspects.</td></tr><tr><td>Wang et al. (2016, 2017)</td><td>Time-domain WRI with approximated extended source.</td></tr><tr><td>Fu & Symes (2017)</td><td>Discrepancy-based penalty method for adaptive penalty parameter.</td></tr><tr><td>Huang et al. (2018a)</td><td>Source-signature independent extended-source FWI.</td></tr><tr><td>Huang et al. (2018b)</td><td>Volume extended-source FWI (first temporal sample of extended source).</td></tr><tr><td>Fang et al. (2018)</td><td>Source signature estimation in WRI.</td></tr><tr><td>Aghamiry et al. (2019c)</td><td>Frequency-domain WRI with Augmented Lagrangian (AL) (IR-WRI).</td></tr><tr><td>Aghamiry et al. (2019a)</td><td>Extension of IR-WRI to VTI media.</td></tr><tr><td>Aghamiry et al. (2019a, 2020b)</td><td>TV & compound regularizations, and bound constraints in IR-WRI.</td></tr><tr><td>Peters & Herrmann (2019)</td><td>Frequency-domain wavefield reconstruction in 3D large-scale models.</td></tr><tr><td>van Leeuwen (2019)</td><td>Recognize WRI as a generalized form of FWI.</td></tr><tr><td>Aghamiry et al. (2020d)</td><td>IR-WRI with phase retrieval.</td></tr><tr><td>Aghamiry et al. (2020c, 2021a)</td><td>IR-WRI for attenuation imaging.</td></tr><tr><td>Symes (2020)</td><td>Mathematical analysis of WRI convexity.</td></tr><tr><td>Lee & Pyun (2020)</td><td>Weighted scattering source minimization in WRI.</td></tr><tr><td>Aghamiry et al. (2021d)</td><td>Source signature estimation in IR-WRI.</td></tr><tr><td>Aghamiry et al. (2021c)</td><td>On the analogy between the contrast source method and WRI.</td></tr><tr><td>Rizzuti et al. (2021)</td><td>Dual formulation of WRI in the data domain.</td></tr><tr><td>Gholami et al. (2022)</td><td>Efficient data-domain algorithm for time-domain IR-WRI.</td></tr><tr><td>Aghazade et al. (2022b)</td><td>IR-WRI with Anderson acceleration.</td></tr><tr><td>Aghazade et al. (2022a)</td><td>IR-WRI with randomized source sketching.</td></tr><tr><td>Hajjaj et al. (2022)</td><td>Wavefield reconstruction in the framework of Marchenko theory.</td></tr><tr><td>Lin et al. (2022)</td><td>Data-domain Hessian with point-spread function.</td></tr><tr><td>Guo et al. (2022)</td><td>Data-domain Hessian with matching filter & truncated iterative method.</td></tr></table>

Accordingly, the objective of ES-FWI is to estimate wavefields that best match the true wavefields by solving a first optimization problem before updating the subsurface parameters from these improved wavefields. Since the true wavefields cannot be estimated in one iteration due to the underdetermined nature of the problem, this two-step cycle is typically iterated until convergence.

This objective raises the underlying questions:

1. How to estimate the best-fitting wavefields from parsimonious data and inaccurate $\mathbf{m}$?

2. How accurately do they match the true wavefields?

3. How do they modify the sensitivity kernel of classical FWI?

4. How does the (limited) accuracy with which they match the true wavefields accounted for in the so-called adjoint source of the FWI gradient?

Answering these questions will draw the similarities and differences between FWI and ES-FWI.

## 2.2 FWI as a constrained optimization problem: full space versus reduced space methods

Before answering these questions, we review how the search-space extension is implemented in ES-FWI from the optimization viewpoint.

Broadly defined, FWI is a partial differential equation constrained optimization problem

$$ \min_{\mathbf{u}_s, \mathbf{m}} \sum_{s=1}^{n_s} \|\mathbf{P}\mathbf{u}_s - \mathbf{d}_s^*\|_2^2 \quad \text{subject to} \quad A(\mathbf{m})\mathbf{u}_s - \mathbf{b}_s = \mathbf{0}, s = 1, \dots, n_s, \qquad (2.4) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Bilinearity of the wave equation. (a) 2004 BP salt model. (b) Three-Hertz monochromatic wavefield triggered by a source in the middle of the model. (c) Reconstructed model from equation 2.3 (Aghamiry et al., 2019c).</figcaption></figure>

where $n_s$ denotes the number of sources, $d_s^*$ the recorded data for source $s$, and the optimization variables gather the wavefields $\mathbf{u}_s$ and the subsurface parameters $\mathbf{m}$. We assume a stationary-recording acquisition (P doesn't depend on $s$) for sake of compact notation but the method can be applied to other acquisition geometries.

In the full space approach, this constrained optimization problem can be solved iteratively with the Lagrange multiplier method

$$ \min_{\mathbf{u}_s, \mathbf{m}} \max_{\mathbf{v}_s} \sum_{s=1}^{n_s} \|\mathbf{P}\mathbf{u}_s - \mathbf{d}_s^*\|_2^2 + \sum_{s=1}^{n_s} \langle \mathbf{v}_s, \mathbf{A}(\mathbf{m})\mathbf{u}_s - \mathbf{b}_s \rangle_{\mathcal{U}}, \quad (2.5) $$

where $\mathbf{v}_s$ denotes the Lagrange multipliers, $\langle \cdot, \cdot \rangle$ the inner product, and $\mathcal{U}$ the state space. The solution of this multivariate problem is found at the saddle point of the above Lagrangian function. When this constrained optimization problem is solved with Newton algorithms, $\mathbf{u}_s$, $\mathbf{m}$ and $\mathbf{v}_s$ are updated jointly at each iteration by solving a normal-equation system, the so-called Karush-Kuhn-Tucker (KKT) system, whose unknowns gather the multivariate descent direction, the right-hand side is the multivariate gradient and the left operator is the multivariate Hessian (e.g., Epanomeritakis et al., 2008). This linear system is prohibitively expensive to solve due to its size.

In classical FWI, the wave-equation constraint is strictly satisfied at each iteration, namely the full space spanned by the three classes of variables is projected onto the parameter space. Accordingly, the wavefields $\mathbf{u}_s$ are eliminated from the optimization variables in equation 2.4 by enforcing their closed-form expression as a function of $\mathbf{m}$ in the data-misfit function following a variable-projection approach (Golub & Pereyra, 2003). This leads to the monovariate data-fitting problem

$$ \min_{\mathbf{m}} \phi^{FWI}(\mathbf{m}) = \sum_{s=1}^{n_s} \|\mathbf{S}(\mathbf{m})\mathbf{b}_s - \mathbf{d}_s^*\|_2^2, \quad (2.6) $$
where $S(\mathbf{m}) = \mathbf{P}\mathbf{A}^{-1}(\mathbf{m})$ is the forward modeling operator. As already mentioned, this data-fitting misfit function is highly multimodal due to its sensitivity to kinematic errors.

In ES-FWI, the search space of classical FWI is extended by re-introducing the wavefields as optimization variables (van Leeuwen & Herrmann, 2016, 2013). This is implemented by processing the wave equation as a soft constraint with penalty methods, equation 2.7, or augmented Lagrangian methods, equation 2.8:

$$ \min_{\mathbf{u}_s, \mathbf{m}} \phi^P(\mathbf{m}, \mathbf{u}_s) = \sum_{s=1}^{n_s} \|\mathbf{P}\mathbf{u}_s - \mathbf{d}_s^*\|_2^2 + \mu \sum_{s=1}^{n_s} \|\mathbf{A}(\mathbf{m})\mathbf{u}_s - \mathbf{b}_s\|_2^2, \quad (2.7) $$

$$ \min_{\mathbf{u}_s, \mathbf{m}} \max_{\mathbf{v}_s} \phi^{AL}(\mathbf{m}, \mathbf{u}_s, \mathbf{v}_s) = \sum_{s=1}^{n_s} \|\mathbf{P}\mathbf{u}_s - \mathbf{d}_s^*\|_2^2 + \mu \sum_{s=1}^{n_s} \|\mathbf{A}(\mathbf{m})\mathbf{u}_s - \mathbf{b}_s\|_2^2 \\ + \sum_{s=1}^{n_s} \langle \mathbf{v}_s, \mathbf{A}(\mathbf{m})\mathbf{u}_s - \mathbf{b}_s \rangle_{\mathcal{U}}, \quad (2.8) $$

where $\mu \in \mathbb{R}_+$ is a penalty parameter. Compared to the penalty function, the Augmented Lagrangian function incorporates a second leverage (the Lagrangian term) that helps to satisfy the constraint accurately at the convergence point even with a fixed $\mu$ through the defect correction action of the Lagrange multipliers $\mathbf{v}_s$. We will review the theoretical aspects of ES-FWI with the penalty formulation for sake of simplicity and because it is enough to understand the physical principles governing ES-FWI beyond the Born approximation. However, the numerical experiments presented at the end of this study will be performed with the Augmented Lagrangian implementation to illustrate the real potential of ES-FWI. The reader is referred to Aghamiry et al. (2019c) and Gholami et al. (2022) for more details about the frequency-domain and time-domain implementations of augmented-Lagrangian based ES-FWI.

From the algorithmic viewpoint, $\mathbf{u}_s$ and $\mathbf{m}$ are not updated jointly because this would be prohibitively computationally expensive. Instead, they are updated in an alternating mode (van Leeuwen & Herrmann, 2013) or by variable projection of the wavefields in the parameter-estimation subproblem like in classical FWI (van Leeuwen & Herrmann, 2016). In the former case, the $\mathbf{m}$-subproblem is linear by virtue of the bilinearity of the wave equation. The relaxation of the wave equation generated by the penalty function implies that the wavefields in the extended search space, referred to as $\mathbf{u}_s^e$ where the subscript $e$ stands for extended, satisfy the wave equation with some errors $\delta\mathbf{b}_s^e$:

$$ \mathbf{A}(\mathbf{m})\mathbf{u}_s^e = \mathbf{b}_s + \delta\mathbf{b}_s^e. \quad (2.9) $$

These errors have however a clear physical meaning in the framework of diffraction theory as we will explain later. Moreover, the penalty function, equation 2.7, implies that each wavefield $\mathbf{u}_s^e$ is the solution of an overdetermined linear system gathering the observation equation $\mathbf{P}\mathbf{u}_s - \mathbf{d}_s^*$ and the wave equation (van Leeuwen & Herrmann, 2013, their equation 6). Accordingly, $\mathbf{u}_s^e$ are reconstructed with a feedback to the data $\mathbf{d}_s^*$ and, hence $\delta\mathbf{b}_s^*$ should involve $\mathbf{d}_s^*$ in one way or another. Accordingly, we call $\mathbf{u}_s^e$ data-assimilated wavefields.

The estimation of the source extensions $\delta\mathbf{b}_s^e$, their relationship with the recorded data $\mathbf{d}_s^*$ and their key role in the parameter estimation are reviewed in the sequel of this tutorial.

## 2.3 The Lippmann-Schwinger equation and the Born approximation

Understanding how data-assimilated wavefields $\mathbf{u}_s^e$ approximate the true wavefields $\mathbf{u}_s^*$ and the role of the source extensions $\delta\mathbf{b}_s^e$ in these wavefield reconstructions requires to review the Lippmann-Schwinger equation (e.g., Prunty & Snieder, 2020).

Let's denote the ground-truth model and the background model by $\mathbf{m}^*$ and $\mathbf{m}$, respectively. The true wavefield $\mathbf{u}^*$ and the background wavefield $\mathbf{u}$ satisfy the wave equation

$$ \mathbf{A}(\mathbf{m}^*)\mathbf{u}^* = \mathbf{b}, \quad (2.10) $$

$$ \mathbf{A}(\mathbf{m})\mathbf{u} = \mathbf{b}. \quad (2.11) $$

Moreover, $\mathbf{u}^*$ can be written as the sum of $\mathbf{u}$ and the true scattered wavefield $\delta\mathbf{u}^*$ by the true perturbation model or scattering object $\delta\mathbf{m}^* = \mathbf{m}^* - \mathbf{m}$:

$$ \mathbf{u}^* = \mathbf{u} + \delta\mathbf{u}^*. \quad (2.12) $$
Equating equations 2.10 and 2.11 and substituting $\mathbf{u}$ by $\mathbf{u}^* - \delta\mathbf{u}^*$, equation 2.12, gives the Lippmann-Schwinger integral equation satisfied by $\delta\mathbf{u}^*$

$$ 
\begin{aligned}
\mathbf{A}(\mathbf{m})\delta\mathbf{u}^* &= -(\mathbf{A}(\mathbf{m}^*) - \mathbf{A}(\mathbf{m}))\mathbf{u}^* \\
&= -\omega^2\mathrm{diag}(\delta\mathbf{m}^*)\mathbf{u}^* \\
&= -\omega^2\sum_i w_i \delta(\mathbf{x} - \mathbf{x}_i) \\
&= \delta\mathbf{b}^*,
\end{aligned} 
\qquad (2.13) $$

where $w_i = \delta\mathbf{m}_i^*\mathbf{u}_i^*$ and $\delta\mathbf{b}^*$ is the true scattering source. This equation shows that the true scattered wavefield $\delta\mathbf{u}^*$ is the weighted superposition of the impulse response of the background $\mathbf{m}$ over the region containing $\delta\mathbf{m}^*$ where the weights $w_i$ are the point wise product between $\delta\mathbf{m}^*$ and $\mathbf{u}^*$ incident to the scatterers (Prunty & Snieder, 2020). The ingredients involved in this volume integral equation (namely, $\mathbf{u}^*$, $\delta\mathbf{m}^*$ and the Green functions $\mathbf{A}^{-1}(\mathbf{m})$) are sketched in Figure 3 for a two-layer medium.

<figure><img src="image_5.png" /><figcaption>Figure 3: Sketch of the volume integral Lippmann-Schwinger equation. The ground-truth model $\mathbf{m}^*$ contains two homogeneous layers while the background medium $\mathbf{m}$ is homogeneous with the velocity of the top layer. Therefore, the domain of integration is the bottom layer as sketched by the fine grid of scatterers (black circles). The volume source is formed by sampling the true wavefields $\mathbf{u}^*$ with the scatterers $\delta\mathbf{m}_i$. The scattered field is then propagated from this source by the Green functions in $\mathbf{m}$, namely $\mathbf{A}^{-1}(\mathbf{m})$ (green lines). In the Born approximation, the incident field to the scatterers is the background wavefield $\mathbf{u}(\mathbf{m})$ (dot black lines) instead of $\mathbf{u}^*$ (solid black lines) leading to a poor approximation of the scattered field (both from the kinematic and dynamic viewpoints) when strong scattering is generated by the perturbation model. The labels $S$ and $R$ denote the source and the receiver, respectively.</figcaption></figure>

Assuming for now that $\delta\mathbf{b}^*$ is known, we can compute the true wavefield $\mathbf{u}^*$ by solving the wave equation in the background $\mathbf{m}$ with an extended source $\mathbf{b} + \delta\mathbf{b}^*$:

$$ \mathbf{A}(\mathbf{m})\mathbf{u}^* = \mathbf{b} + \delta\mathbf{b}^*. \qquad (2.14) $$

In inverse scattering problems, the scattered wavefields are classically computed with the Born approximation by replacing the true wavefield $\mathbf{u}^*$ by the background wavefield $\mathbf{u}(\mathbf{m})$ in the weights $w_i$ of the scattering source $\delta\mathbf{b}$ (Figure 3). This approximation is only valid for weak scattering. Therefore, ES-FWI aims at extending the linear regime of FWI beyond the Born approximation to avoid cycle skipping by estimating an improved approximation of the scattering source $\delta\mathbf{b}^*$.

## 2.4 The Lippmann-Schwinger equation and classical FWI

In the framework of local-optimization (gradient-based) methods, classical FWI recasts the non linear data fitting problem as a recurrence of surrogate linear problems, whose unknowns are the model perturbations with which the background model $\mathbf{m}$ of the current iteration will be updated and the observables are the data residuals $\delta\mathbf{d}^*(\mathbf{m}) =$
$d^* - d(m)$, where $d(m) = S(m)b$ (Figure 4). In this framework, $\delta d^*(m)$ should be understood as the recorded scattered data by the missing model perturbation $\delta m^*$ in $m$ (Pratt et al., 1998). Accordingly, the data residuals are the restriction at receivers of the true scattered wavefields $\delta u^*$ by $\delta m^*$. Applying the observation operator $P$ to $\delta u^*$ in equation 2.13 gives

$$ \delta d^*(m) = S(m)\delta b^* = -\omega^2 S(m) \sum_i w_i \delta(x - x_i), \qquad (2.15) $$

where it is reminded that $w_i = \delta m_i^* u_i^*$. The ground-truth scattered data are then the recording at receivers of the blended impulse response of the background $m$ over the region containing $\delta m^*$ where the unknown weights depend both on $\delta m^*$ and $u^*$. The nonlinear Equation 2.15 can be rewritten as

$$ \delta d(\delta m) = B(m; u(m + \delta m))\delta m, \qquad (2.16) $$

where $B(m; u(m + \delta m)) = -\omega^2 S(m) \text{diag}(u(m + \delta m))$. Equation 2.15 can be viewed as a reparametrization of the nonlinear forward equation of FWI, $d = S(m)b$. In classical FWI, the simulated scattered data are computed after a linearization of this nonlinear forward equation in the frame of the Born approximation, namely $B(m; u(m)) \approx B(m; u(m + \delta m))$.

<figure><img src="image_6.png" /><figcaption>Figure 4: Sketch of the multimodal nonlinear FWI data fitting misfit function $C^{FWI}(m)$ (blue curve), equation 2.6, and the sequence of surrogate quadratic functions $\tilde{C}_{m_k}(m)$ built at each iteration of Newton-type algorithms (orange curves) by linearization around the background models $m_k$, where $\tilde{C}_{m_k}(m) = C(m_k) + \delta m^T \nabla_m C(m_k) + \frac{1}{2} \delta m^T \nabla_m^2 C(m_k) \delta m$ with $m = m_k + \delta m$ (e.g., Aghamiry et al., 2021b).</figcaption></figure>

The gradient of FWI is formed by the inner product between each column of the sensitivity matrix and the data residuals (e.g., Pratt et al., 1998):

$$ \nabla_m \phi^{FWI}(m) = J(m)^T \delta d(m)^*, \qquad (2.17) $$

where

$$ J = \begin{bmatrix} \frac{\partial d_1}{\partial m_1} & \frac{\partial d_1}{\partial m_2} & \cdots & \frac{\partial d_1}{\partial m_M} \\ \frac{\partial d_2}{\partial m_1} & \frac{\partial d_2}{\partial m_2} & \cdots & \frac{\partial d_2}{\partial m_M} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial d_j}{\partial m_1} & \frac{\partial d_j}{\partial m_2} & \cdots & \frac{\partial d_j}{\partial m_M} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial d_N}{\partial m_1} & \frac{\partial d_N}{\partial m_2} & \cdots & \frac{\partial d_N}{\partial m_M} \end{bmatrix}, \quad \delta d(m)^* = \begin{bmatrix} \delta d_1(m) \\ \delta d_2(m) \\ \vdots \\ \delta d_j(m) \\ \vdots \\ \delta d_N(m) \end{bmatrix}. $$

For one frequency, $M$ and $N = n_s \times n_r$ denote the number of grid points sampling $m$ and the number of data sampling the stationary-recording acquisition, respectively.

The simulated scattered data $\frac{\partial d}{\partial m_i}$ are the restriction at receivers of the so-called partial derivative wavefields, which are the solution of an approximate Lippmann-Schwinger equation for a point source located at the position of the model parameter $i$ with respect to which the partial derivative is computed (Pratt et al., 1998, their equation 15). For column $i$,

$$ \partial_{m_i} d(m) = -\omega^2 S(m) w_i^{FWI} \delta(x - x_i), \qquad (2.18) $$

where $w_i^{FWI} = u_i$. Comparing equations 2.15 and 2.18 shows that the sensitivity matrix can be interpreted as a deblending operator whose columns aim at picking in the weights $w_i$ of the blended scattering source that model
perturbation $\delta \mathbf{m}_i^*$ to be mapped at the position of $\mathbf{m}_i$. The key difference between equations 2.15 and 2.18 lies indeed in the wavefield incident to the scatterers in their scattering sources (the true wavefield and the background wavefield in $w_i$ and $w_i^{FWI}$, respectively). This difference results from the weak-scattering Born approximation with which the simulated scattered data are computed as above mentioned. This approximation is poor when $\delta \mathbf{m}$ generates strong scattering along the incident path connecting the source to the scatterers and prevents accurate picking of the model perturbation by the Born-simulated scattered data as sketched in Figure 3 and illustrated numerically in Figure 5.

<figure><img src="image_7.png" /><figcaption>Figure 5: Scattered data and the Born approximation. Inset: The ground-truth model $\mathbf{m}^*$ (blue curve) is a vertical velocity-gradient model while the background medium $\mathbf{m}$ is homogeneous (red curve). (a) Ground-truth scattered data: $\delta \mathbf{d}^*(\mathbf{m}) = \mathbf{d}^* - \mathbf{d}(\mathbf{m}) = -\omega^2 \mathbf{S}(\mathbf{m})\text{diag}(\mathbf{m}^*)\mathbf{u}^*$. (b) Scattered data in the Born approximation: $\delta \mathbf{d}^{\text{Born}}(\mathbf{m}) = -\omega^2 \mathbf{S}(\mathbf{m})\text{diag}(\mathbf{m}^*)\mathbf{u}$. (c) Scattered data using $\mathbf{u}^e$ in the scattering source: $\delta \mathbf{d}(\mathbf{m}) = -\omega^2 \mathbf{S}(\mathbf{m})\text{diag}(\mathbf{m}^*)\mathbf{u}^e$.</figcaption></figure>

According to the above review, we will show that the key contribution of ES-FWI is simply to implement more accurate wavefields in the virtual source of the partial derivative data beyond the Born approximation.

## 2.5 Estimating more accurate wavefields from the Lippmann-Schwinger equation

We estimate the wavefields that best approximate $\mathbf{u}^*$ with the Lippmann-Schwinger equation because this equation gathers all the available priors about $\mathbf{u}^*$, that is the propagator $\mathbf{m}$ and the recorded data $\mathbf{d}^*$, i.e., the restriction of $\mathbf{u}^*$ at receivers. Let's remind that the true wavefield satisfies the wave equation

$$ \mathbf{A}(\mathbf{m})\mathbf{u}^* = \mathbf{b} + \delta\mathbf{b}^*, \qquad (2.19) $$

where $\delta\mathbf{b}^*$ is the unknown of the following rank-deficient forward equation

$$ \delta \mathbf{d}^* = \mathbf{S}(\mathbf{m})\delta \mathbf{b}^*. \qquad (2.20) $$

We can compute an approximation of $\mathbf{u}^*$ by solving the wave equation, equation 2.19, once we have estimated an approximation of $\delta\mathbf{b}^*$. The best we can do to estimate $\delta\mathbf{b}^*$ is to solve the rank-deficient equation 2.20 in a least-squares sense

$$ \delta \mathbf{b}^e = \arg \min_{\delta \mathbf{b}} \|\delta \mathbf{d}^*(\mathbf{m}) - \mathbf{S}(\mathbf{m})\delta \mathbf{b}\|_2^2 + \mu \|\delta \mathbf{b}\|_2^2, \qquad (2.21) $$

where a damping regularization was introduced for stabilization. This problem is nothing more than the damped least-squared scattered-data fitting problem for scattering source estimation.

The solution of this source-dependent underdetermined problem is given by

$$ \delta \mathbf{b}^e(\mathbf{m}) = \mathbf{S}^T(\mathbf{m}) (\mathbf{S}(\mathbf{m})\mathbf{S}^T(\mathbf{m}) + \mu \mathbf{I})^{-1} \delta \mathbf{d}^* = \mathbf{S}^T(\mathbf{m}) \mathbf{H}_d^{-1}(\mathbf{m}) \delta \mathbf{d}^*, \qquad (2.22) $$

where the Hessian is formulated in the data domain at the expense of the source domain since the dimension of the data space spanned by one source (here, $n_r$ when considering one frequency) is much smaller than the domain where the wavefields are computed. Moreover, let's note that $\delta \mathbf{b}^e$ is related to $\delta \mathbf{b}^*$ by the source resolution matrix $\mathbf{R}(\mathbf{m})$,

$$ \delta \mathbf{b}^e = (\mathbf{S}(\mathbf{m})^T \mathbf{S}(\mathbf{m}) + \mu \mathbf{I})^{-1} \mathbf{S}^T(\mathbf{m}) \delta \mathbf{d}^* = (\mathbf{S}(\mathbf{m})^T \mathbf{S}(\mathbf{m}) + \mu \mathbf{I})^{-1} \mathbf{S}^T(\mathbf{m}) \mathbf{S}(\mathbf{m}) \delta \mathbf{b}^* = \mathbf{R}(\mathbf{m}) \delta \mathbf{b}^*, \qquad (2.23) $$
and, the simulated scattered data $\delta \mathbf{d}^e$ by $\delta \mathbf{b}^e$ are related to the recorded scattered data $\delta \mathbf{d}^*$ by the data-domain resolution matrix $\mathbf{R}_d(\mathbf{m})$:

$$ \delta \mathbf{d}^e = \mathbf{S}(\mathbf{m})\delta \mathbf{b}^e = \mathbf{S}(\mathbf{m})\mathbf{S}^T(\mathbf{m}) (\mathbf{S}(\mathbf{m})\mathbf{S}^T(\mathbf{m}) + \mu \mathbf{I})^{-1} \delta \mathbf{d}^* = \mathbf{R}_d(\mathbf{m})\delta \mathbf{d}^*. \quad (2.24) $$

It is important to stress that the underdetermination of the scattering-source estimation not only results from the incomplete subsurface illumination provided by the receiver layout but more importantly by that fact that the sought $\delta \mathbf{b}^*$ depends on two classes of unknowns, $\delta \mathbf{m}^*$ and $\mathbf{u}^*$, with potential coupling between the two. From the numerical viewpoint, equation 2.22 shows that the scattering sources $\delta \mathbf{b}_s^e$ can be computed in two steps:

For $s = 1, ..., n_s$,

1. Solve approximately or exactly the normal system

$$ (\mathbf{S}(\mathbf{m})\mathbf{S}^T(\mathbf{m}) + \mu \mathbf{I}) \mathbf{y} = \delta \mathbf{d}_s^*. \quad (2.25) $$

2. Solve the adjoint wave equation

$$ \mathbf{A}^T(\mathbf{m})\delta \mathbf{b}_s^e(\mathbf{m}) = \mathbf{P}^T\mathbf{y} \quad (2.26) $$

with the weighted data residuals as source. Note that the left operator in equation 2.25 is independent to the sources for a stationary-recording acquisition.

The scattering sources $\delta \mathbf{b}_s^e$ show obvious analogy with the so-called adjoint wavefield of classical FWI, the difference being related to the weighting operator applied to the data residuals before back-propagation. The role of this weighting operator will be discussed in more detail later. Finally, one can readily see that the $\delta \mathbf{b}$ problem, equation 2.21, is nothing more that a reparametrization of wavefield reconstruction problem, equation 2.7.

$$ [\mathbf{u} \to \delta \mathbf{b}] \Rightarrow \left[ \min_{\mathbf{u}} \|\mathbf{d}^* - \mathbf{P}\mathbf{u}\|_2^2 + \mu \|\mathbf{A}(\mathbf{m})\mathbf{u} - \mathbf{b}\|_2^2 \to \min_{\delta \mathbf{b}} \|\delta \mathbf{d}^* - \mathbf{S}(\mathbf{m})\delta \mathbf{b}\|_2^2 + \mu \|\delta \mathbf{b}\|_2^2 \right]. \quad (2.27) $$

Therefore, the best-fitting wavefield $\mathbf{u}^e$ in the extended space can be inferred from the best-fitting scattering source $\delta \mathbf{b}^e$ by solving the wave equation in the background $\mathbf{m}$ with the following extended source $\mathbf{b}^e$:

$$ \mathbf{A}(\mathbf{m})\mathbf{u}^e = \mathbf{b}^e = \mathbf{b} + \delta \mathbf{b}^e(\mathbf{m}) = \mathbf{b} + \mathbf{R}(\mathbf{m})\delta \mathbf{b}^*(\mathbf{m}). \quad (2.28) $$

Analogy between equations 2.19 and 2.28 shows that $\delta \mathbf{b}^e$ accounts approximately for the scattering generated by the unknown $\delta \mathbf{m}^*$ contained in $\delta \mathbf{b}^*$, equation 2.13, where the approximation is represented by the source-domain resolution matrix $\mathbf{R}(\mathbf{m})$. Figure 5c shows that using $\mathbf{u}^e$ in the scattering source of the Lippmann-Schwinger equation in place of the background wavefield $\mathbf{u}(\mathbf{m})$ allows for the accurate modeling of the ground-truth scattered data in presence of strong scattering unlike the Born approximation (Figure 5b).

## 2.6 Parameter estimation beyond the Born approximation

Now that we saw how to reconstruct wavefields that best match the true wavefields, equation 2.28, we review how they modify the gradient with respect to $\mathbf{m}$ of the classical FWI misfit function. Injecting the expression of $\mathbf{u}^e$, equation 2.28, and $\delta \mathbf{b}^e$, equation 2.22, in the ES-FWI objective function, equation 2.7, and reorganizing the terms recasts the later as a generalized form of the classical FWI misfit function (Gholami et al., 2022; Symes, 2020; van Leeuwen, 2019)

$$ \min_{\mathbf{m}} \sum_{s=1}^{n_s} \|\delta \mathbf{d}_s^*(\mathbf{m})\|_{\mathbf{H}_d^{-1}(\mathbf{m})}^2, \quad (2.29) $$

where $\|\mathbf{y}\|_{\mathbf{Q}}^2 = \mathbf{y}^T \mathbf{Q} \mathbf{y}$.

The data-domain Hessian $\mathbf{H}_d(\mathbf{m})$ acts as a data-domain covariance matrix, which however depends on $\mathbf{m}$. Let's remind that

$$ \mathcal{C}(\mathbf{m}) = \|\delta \mathbf{d}(\mathbf{m})\|_{\mathbf{Q}}^2 \to \nabla \mathcal{C}(\mathbf{m}) = \left[ \frac{\partial \mathbf{d}(\mathbf{m})}{\partial \mathbf{m}} \right]^T \mathbf{Q} \delta \mathbf{d}(\mathbf{m}). \quad (2.30) $$

The above equation shows that introducing a covariance matrix in the FWI data misfit function modifies the adjoint source (i.e., the data residuals) in its gradient. However, as $\mathbf{H}_d^{-1}(\mathbf{m})$ depends on $\mathbf{m}$, it not only modifies the source of the adjoint equation but also the kinematic of the simulated scattered data in the sensitivity matrix. This can be
easily shown by computing the gradient of the ES-FWI objective function with respect to **m** by keeping **u**<sup>e</sup> fixed and independent of **m**. This amounts to minimize the scattering sources δ**b**<sub>s</sub><sup>e</sup>(**m**) = **A**(**m**)<**u**<sub>s</sub><sup>e</sup>** - **b**<sub>s</sub> in a least-squares sense:

$$
\min_{\mathbf{m}} \sum_{s=1}^{n_s} \|\delta \mathbf{b}_s^e(\mathbf{m})\|_2^2. \qquad (2.31)
$$

Remark 1: We drop the penalty parameter µ as a multiplicative factor in the above equation because it is useless since the gradient has not the units of the parameters and hence a metric (step length and/or Hessian) should be defined accordingly to design the descent direction of the m-subproblem (see also Remark 4).

The gradient is given by

$$
\nabla_{\mathbf{m}} \sum_{s=1}^{n_s} \|\delta \mathbf{b}_s^e(\mathbf{m})\|_2^2 = \nabla_{\mathbf{m}} \sum_{s=1}^{n_s} \|\mathbf{A}(\mathbf{m}) \mathbf{u}_s^e - \mathbf{b}_s\|_2^2 = \sum_{s=1}^{n_s} \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right)^T \delta \mathbf{b}_s^e. \quad (2.32)
$$

Substituting **u**<sub>s</sub><sup>e</sup> by **u**<sub>s</sub> + **A**<sup>-1</sup>(**m**)<b>R</b><sub>s</sub>(**m**)<b>δb</b><sub>s</sub><sup>*</sup>(**m**)<sup>*</sup> and <b>δb</b><sub>s</sub><sup>e</sup> by <b>S</b><sup>T</sup>(**m**)<b>H</b><sub>d</sub><sup>-1</sup>(**m**)<b>δd</b><sub>s</sub><sup>*</sup>(**m**)<sup>*</sup> in the above equation highlights
the similarities and differences between FWI and ES-FWI gradients:

$$
\begin{align}
\nabla_{\mathbf{m}} \sum_{s=1}^{n_s} \|\delta \mathbf{b}_s^e(\mathbf{m})\|_2^2 &= \sum_{s=1}^{n_s} \left( \underbrace{\mathbf{S}(\mathbf{m}) \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}}}_{\text{single-scattering}} \right)^T \mathbf{H}_d^{-1}(\mathbf{m}) \delta \mathbf{d}_s^* \nonumber \\
&\quad + \sum_{s=1}^{n_s} \left( \underbrace{\mathbf{S}(\mathbf{m}) \frac{\partial \mathbf{A}(\mathbf{m})}{\partial \mathbf{m}}}_{\text{multi-scattering}} \frac{\mathbf{A}^{-1}(\mathbf{m}) \mathbf{R}_s(\mathbf{m}) \delta \mathbf{b}_s^*}{\delta \mathbf{u}_s^e \approx \delta \mathbf{u}_s^*} \right)^T \mathbf{H}_d^{-1}(\mathbf{m}) \delta \mathbf{d}_s^*. \tag{2.33}
\end{align}
$$

The first term is the gradient of classical FWI except that the data residuals in the adjoint source are weighted by $\mathbf{H}_d^{-1}(\mathbf{m})$, the inverse of the data-domain Hessian of the $\delta\mathbf{b}$ estimation problem, equation 2.21. The term in brackets represents the partial derivative data of FWI computed with the weak-scattering Born approximation. That is, the so-called virtual sources are built by sampling the incident background wavefields $\mathbf{u}_s$ to the scatterers with the radiation pattern matrix $\partial\mathbf{A}(\mathbf{m})/\partial\mathbf{m}$ from which the single-scattered partial derivative wavefields are propagated and recorded at receivers through the forward modeling operator $\mathbf{S}(\mathbf{m})$ (Pratt et al., 1998). The second term supplements the single-scattered component of the partial derivative data with an approximation of the multi-scattered counterpart. This is shown by reading this term from right to left. An approximation $\delta\mathbf{u}_s^e$ of the true scattered wavefields $\delta\mathbf{u}_s^*$,

$$
\delta \mathbf{u}_s^e = \mathbf{A}^{-1}(\mathbf{m}) \delta \mathbf{b}_s^e = \mathbf{A}^{-1}(\mathbf{m}) \mathbf{R}(\mathbf{m}) \delta \mathbf{b}_s^*, \quad (2.34)
$$

are triggered by δ**b**<sub>s</sub><sup>e</sup> = **R**(**m**)<b>δb</b><sub>s</sub><sup>*</sup>, i.e., the approximation of the true scattering sources δ**b**<sub>s</sub><sup>*</sup>. These incident scattered
wavefields are then sampled by the radiation pattern matrix ∂**A**(**m**)/∂**m** to build the multi-scattered component of
the virtual sources from which the multi-scattered component of the partial derivative wavefields are propagated and
recorded at receivers via the forward modeling operator **S**(**m**). The approximation with which scattering generated
by the sought model perturbation is estimated in the incident wavefields is highlighted by the resolution matrix **R**(**m**)
of the scattering-source estimation problem, equation 2.23.

This analogy between FWI and ES-FWI highlights that ES-FWI can be easily implemented in a classical FWI code although the extra ingredients of ES-FWI can generate significant albeit manageable computational overheads. These ingredients are the scattered wavefields added to the background wavefields in the virtual sources and the inverse of the data-domain Hessian in the adjoint source. These two ingredients can be easily turned off in the code to switch from ES-FWI to FWI once the updated model lies in the basin of attraction of the global minimizer.

Remark 2: The attentive reader will have noted that the estimation of the true scattering sources, namely δb<sub>s</sub><sup>e</sup>, are involved two times in the ES-FWI algorithm. First, they are used to supplement the background wavefield with the scattered wavefield by the missing model perturbation in the background model to build an estimation of the true wavefields, equation 2.28. These wavefields are then used to build more accurate virtual sources of the partial derivative data, term in brackets in equation 2.32. Second, δb<sub>s</sub><sup>e</sup> are used as the so-called adjoint wavefields, i.e., the estimation of the true scattering sources, which are correlated with the simulated virtual sources to build the gradient of ES-FWI.
## 2.7 Parameter estimation as a linear problem

In the previous section, we formulate the parameter-estimation subproblem in a way that would highlight the similarities and differences between the steepest-descent direction of FWI and ES-FWI. Alternatively, when the wavefields and the model parameters are updated in alternating mode, the subsurface model can be obtained directly taking advantage of the bilinearity of the wave equation reviewed at the beginning of this tutorial. Assuming that the wavefields are known and independent of $\mathbf{m}$, the objective function of ES-FWI for model estimation can be recast as a quadratic function

$$ \min_{\mathbf{m}} \sum_{s=1}^{n_s} \|\mathbf{L}(\mathbf{u}_s^e)\mathbf{m} - \mathbf{y}_s(\mathbf{u}_s^e)\|_2^2, \qquad (2.35) $$

where $\mathbf{L}(\mathbf{u}_s^e) = \omega^2 \mathrm{diag}(\mathbf{u}_s^e)$ and $\mathbf{y}_s = \mathbf{b}_s - \nabla^2 \mathbf{u}_s^e$, equation 2.2 (Aghamiry et al., 2019c; van Leeuwen & Herrmann, 2013).

The solution of this overdetermined problem reads

$$ \mathbf{m} = \frac{\sum_{s=1}^{n_s} \mathbf{L}^T(\mathbf{u}_s^e)\mathbf{y}_s(\mathbf{u}_s^e)}{\sum_{s=1}^{n_s} \mathbf{L}^T(\mathbf{u}_s^e)\mathbf{L}(\mathbf{u}_s^e)}, \qquad (2.36) $$

which is indeed the least-squares analogue of the expression of $\mathbf{m}$ inferred from the bilinearity of the wave equation when a monochromatic wavefield triggered by a single source is known everywhere in the subsurface, equation 2.3. The denominator of the right-hand side of equation 2.36 shows that the Gauss-Newton Hessian is diagonal and is formed by the zero-lag autocorrelation of the virtual sources in the extended search space. It corrects the descent direction of ES-FWI from geometrical spreading along the incident paths connecting the sources to the scatterers (Shin et al., 2001). The missing component of the Gauss-Newton FWI Hessian accounting for the paths connecting the scatterers to the receivers has been moved in the wavefield estimation subproblem through the data-domain Hessian of the scattering-source estimation problem. Disregarding the damping term, this Hessian is formed by the normal matrix $\mathbf{S}\mathbf{S}^T$ where $\mathbf{S}$ describes the propagation from the scatterers to the receivers. Therefore, $\mathbf{S}\mathbf{S}^T$ describes the summation over scatterers of the correlation in time and space between Green functions recorded at different receivers. Applying the inverse of this operator to the recorded scattered data aims at deconvolving the later from the limited bandwidth effects generated by these correlations (Gholami et al., 2022, Their appendix A). We will show however later that this deconvolution has not the same effect on the recorded data and the simulated data in the data residuals when the background $\mathbf{m}$ is highly inaccurate.

## 2.8 Summary of the overall ES-FWI workflow

As a summary, we recap the two step-procedure of ES-FWI. The starting point is to recognize that the true wavefields satisfy a wave equation that can be written in two different forms:

$$ \mathbf{A}(\mathbf{m}^*)\mathbf{u}^* = \mathbf{b}, \qquad \text{Classical wave equation} \qquad (2.37) $$

$$ \mathbf{A}(\mathbf{m})\mathbf{u}^* = \mathbf{b} + \delta\mathbf{b}^*, \qquad \text{Lippmann-Schwinger volume integral equation} \qquad (2.38) $$

where $\delta\mathbf{b}^*$ is the unknown of the rank-deficient forward-problem equation

$$ \delta\mathbf{d}^* = \mathbf{S}(\mathbf{m})\delta\mathbf{b}^*. \qquad (2.39) $$

*Wavefield reconstruction: Pushing $u$ towards $u^*$ with Lippmann-Schwinger equation*

According to the Lippmann-Schwinger equation satisfied by $\mathbf{u}^*$ and $\mathbf{u}^e$,

$$ \mathbf{A}(\mathbf{m})\mathbf{u}^* = \mathbf{b} + \delta\mathbf{b}^*, \qquad (2.40) $$

$$ \mathbf{A}(\mathbf{m})\mathbf{u}^e = \mathbf{b} + \delta\mathbf{b}^e, \qquad (2.41) $$

we seek $\mathbf{u}^e$ that best fits $\mathbf{u}^*$ by looking for $\delta\mathbf{b}^e$ that best fits $\delta\mathbf{b}^*$. From equation 2.39, $\delta\mathbf{b}^e$ is defined as the damped least-squares solution of the scattered-data fitting problem:

For $s = 1, \dots, n_s$,

$$ \delta\mathbf{b}_s^e = \arg\min_{\delta\mathbf{b}_s} \|\delta\mathbf{d}_s^*(\mathbf{m}) - \mathbf{S}(\mathbf{m})\delta\mathbf{b}_s\|_2^2 + \mu\|\delta\mathbf{b}_s\|_2^2, \qquad (2.42) $$

which gives

$$ \delta\mathbf{b}_s^e = \mathbf{S}^T(\mathbf{m})\mathbf{H}_d^{-1}(\mathbf{m})\delta\mathbf{d}_s^*(\mathbf{m}). \qquad (2.43) $$
Parameter estimation: Pushing m towards m* by minimizing the scattering sources

According to

$$A(\mathbf{m}^*)\mathbf{u}^* = \mathbf{b} + \mathbf{0}, \qquad (2.44)$$

$$A(\mathbf{m})\mathbf{u}^e = \mathbf{b} + \delta\mathbf{b}^e, \qquad (2.45)$$

and assuming that $\mathbf{u}^e$ is a good approximation of $\mathbf{u}^*$, pushing $\mathbf{m}$ towards $\mathbf{m}^*$ amounts to jointly minimize $\delta\mathbf{b}^e$ for each source $\mathbf{b}$ of the experiment while minimizing the residual data misfit:

$$\min_{\mathbf{m}} \sum_{s=1}^{n_s} \|\mathbf{d}_s^*(\mathbf{m}) - \mathbf{P}\mathbf{u}_s^e(\mathbf{m})\|_2^2 + \mu \sum_{s=1}^{n_s} \|\mathbf{A}(\mathbf{m})\mathbf{u}_s^e(\mathbf{m}) - \mathbf{b}_s\|_2^2, \qquad (2.46)$$

which is equivalent to

$$\min_{\mathbf{m}} \sum_{s=1}^{n_s} \|\delta \mathbf{d}_s^*(\mathbf{m})\|_{\mathbf{H}_d^{-1}(\mathbf{m})}^2. \qquad (2.47)$$

If $\mathbf{u}_s^e$ and $\mathbf{m}$ are updated in alternating mode rather than through a variable projection method, the $\mathbf{m}$-subproblem reduces to the minimization of the scattering sources (adjoint wavefields),

$$\min_{\mathbf{m}} \sum_{s=1}^{n_s} \|\delta \mathbf{b}_s^e(\mathbf{m})\|_2^2, \qquad (2.48)$$

which is quadratic in $\mathbf{m}$ by virtue of the bilinearity of the wave equation, equation 2.46, and amounts to minimize the data residuals $\delta\mathbf{d}^*(\mathbf{m})$ according to equation 2.43.

The gradient of the two objective functions, equations 2.47 and 2.48, are the same, while the Hessian is not (Gholami et al., 2022; van Leeuwen & Herrmann, 2016).

## 2.9 On the accuracy of the reconstructed wavefields in extended search space

Understanding the limits of ES-FWI and designing relevant heuristics to drive it toward accurate minimizer require additional insights on the accuracy with which the reconstructed wavefields match the true counterparts. Let's first remind the expression of $\mathbf{u}^e$ from equations 2.28 and 2.22:

For $s = 1, ..., n_s$,

$$A(\mathbf{m})\mathbf{u}_s^e = \mathbf{b}_s^e(\mathbf{m}) = \mathbf{b}_s + \mathbf{S}^T(\mathbf{m}) (\mathbf{S}(\mathbf{m})\mathbf{S}^T(\mathbf{m}) + \mu\mathbf{I})^{-1} \delta \mathbf{d}_s^*(\mathbf{m}) \qquad (2.49)$$

and let's denote the simulated data in the extended space by $\mathbf{d}_s^e = \mathbf{P}\mathbf{u}_s^e$.

Reordering the terms of equation 2.49 provides the following identities:

For $s = 1, ..., n_s$,

$$\mathbf{b}_s^e(\mathbf{m}) = (\mathbf{I} - \mathbf{R}(\mathbf{m}))\mathbf{b}_s + \mathbf{S}(\mathbf{m})^T \mathbf{H}_d^{-1}(\mathbf{m}) \mathbf{d}_s^*, \qquad (2.50)$$

$$\mathbf{u}_s^e(\mathbf{m}) = \mathbf{A}^{-1}(\mathbf{m})(\mathbf{I} - \mathbf{R}(\mathbf{m}))\mathbf{b}_s + \mathbf{A}^{-1}(\mathbf{m})\mathbf{A}(\mathbf{m})^{-T}\mathbf{P}^T\mathbf{H}_d^{-1}(\mathbf{m})\mathbf{d}_s^*, \qquad (2.51)$$

$$\mathbf{d}_s^e(\mathbf{m}) = (\mathbf{I} - \mathbf{R}_d(\mathbf{m}))\mathbf{d}_s(\mathbf{m}) + \mathbf{R}_d(\mathbf{m})\mathbf{d}_s^*, \qquad (2.52)$$

$$\Delta \mathbf{d}_s^e(\mathbf{m}) = \mu (\mathbf{S}(\mathbf{m})\mathbf{S}^T(\mathbf{m}) + \mu \mathbf{I})^{-1} \delta \mathbf{d}_s^*(\mathbf{m}) = (\mathbf{I} - \mathbf{R}_d(\mathbf{m})) \delta \mathbf{d}^*, \qquad (2.53)$$

where $\Delta \mathbf{d}_s^e(\mathbf{m}) = \mathbf{d}_s^e - \mathbf{d}_s^e(\mathbf{m})$ are the data residuals in the extended search space, which should not be confused with the simulated scattered data in the extended search space, $\delta \mathbf{d}_s^e = \mathbf{d}_s^e(\mathbf{m}) - \mathbf{d}(\mathbf{m}) = \mathbf{S}(\mathbf{m})\delta \mathbf{b}_s^e = \mathbf{R}_d(\mathbf{m})\delta \mathbf{d}^*$, equation 2.24.

The extended sources $\mathbf{b}_s^e(\mathbf{m})$, equation 2.50, are the weighted sum of $\mathbf{b}_s$ and the back-propagated weighted recorded data $\mathbf{H}_d^{-1}(\mathbf{m})\mathbf{d}_s^*$. When $\mu$ is set to a small value, the source resolution matrix $\mathbf{R} \approx \mathbf{I}$ and $\mathbf{b}_s^e(\mathbf{m})$ mostly correspond to the back propagated weighted recorded data. The reconstruction of the extended source $\mathbf{b}_s^e$ is illustrated in Figure 6 for the velocity models shown in Figure 5. Figure 6b shows the back propagated simulated data. Due to the kinematic consistency of the simulated data $\mathbf{d}(\mathbf{m})$ with the background $\mathbf{m}$, the back-propagated simulated data optimally match the point source $\mathbf{b}$ in space and time. Conversely, the back-propagated recorded data are smeared in time and space as the witness of the inaccuracy of $\mathbf{m}$ (Figure 6a). When the weighted data residuals are back-propagated and summed with $\mathbf{b}$, the sign-reversed simulated component of the back-propagated data residuals cancels out with $\mathbf{b}$ by subtraction mostly leaving the back propagated recorded data in $\mathbf{b}^e$ (Figure 6c).

**Remark 3:** The back-propagation of the recorded data to estimate the source in the extended search space is more directly highlighted when ES-FWI is parametrized with the source rather than with the wavefield or the scattering
<figure><img src="image_8.png" /><figcaption>Figure 6: Scattering source. Snapshots of (a) the weighted back-propagated recorded data, (b) the weighted back-propagated simulated data, and (c) the weighted back-propagated data residuals.</figcaption></figure>

source (Huang et al., 2018a, their equations 10 and 11). However, the connection with inverse scattering theory and hence classical FWI is less obvious with the source parametrization.

From the expression of the extended sources, the extended wavefields $\mathbf{u}_s^e$, equation 2.51, can be formulated as the weighted sum of the background wavefields and the migrated-demigrated weighted recorded data (the migration-demigration operator is given by $\mathbf{A}^{-1}(\mathbf{m})\mathbf{A}(\mathbf{m})^{-T}\mathbf{P}^T$). When $\mu$ is set to a small value, $\mathbf{R}_d(\mathbf{m}) \approx \mathbf{I}$ and the extended wavefields mostly reduce to the migrated-demigrated weighted recorded data. Accordingly, they match the recorded data at receivers by virtue of wave invariance for time reversal. However, they accumulate kinematic errors as they propagate away from the receivers in the background $\mathbf{m}$. This is illustrated in Figure 7, which shows snapshots of the true wavefield $\mathbf{u}^*$, the background wavefield $\mathbf{u}$ and the data-assimilated wavefield $\mathbf{u}^e$ for the velocity models of Figure 5. One can see that $\mathbf{u}^e$ matches $\mathbf{u}^*$ only near the receivers when $\mathbf{m}$ differs significantly from $\mathbf{m}^*$ while it fits $\mathbf{u}$ in the deep part with a complex transition between the shallow to the deep regions.

<figure><img src="image_9.png" /><figcaption>Figure 7: Data-assimilated wavefield. Snapshots of (a) the true wavefield $\mathbf{u}^*$, (b) the background wavefield $\mathbf{u}$, and (c) the data-assimilated wavefield $\mathbf{u}^e$.</figcaption></figure>

From the decomposition of $\mathbf{u}_s^e$, equation 2.51, the simulated data $\mathbf{d}_s^e$ in the extended space are the convex combination of the simulated data in $\mathbf{m}$, $\mathbf{d}_s(\mathbf{m}) = \mathbf{S}(\mathbf{m})\mathbf{b}_s$, and the recorded data $\mathbf{d}_s^*$ where the weighting operator is the data-resolution matrix. Again, setting $\mu$ to a small value pushes $\mathbf{d}^e$ towards $\mathbf{d}_s^*$ at the expense of $\mathbf{d}_s(\mathbf{m})$ to prevent cycle skipping. The ability to match the recorded data from arbitrary background model $\mathbf{m}$ is further checked numerically
in Figure 8 with the velocity models shown in Figure 5. The idea to push the simulated data towards the recorded data (e.g., Baek et al., 2014) or the reverse (e.g., Yao et al., 2019) is shared by almost all the approaches aiming at avoiding cycle skipping with the distinctiveness of ES-FWI to be physics driven.

Finally, equation 2.53 shows that the data residuals in the extended space are linearly related to the recorded scattered data via the data-domain Hessian of the $\delta b$ subproblem and further scaled by the penalty parameter $\mu$. This scaling highlights how the amplitudes of $\Delta d_s^e$ decrease (i.e., how $d_s^e$ can arbitrarily fit the recorded data) as $\mu$ decreases (Figure 8). The effect of the data-domain Hessian on the data residuals is discussed in the next section.

**Remark 4:** The fact that $d_s^e$ closely match $d_s^*$ via the migration-demigration of the FWI data residuals, hence $\Delta d_s^e$ are small, doesn't indeed mean that $\delta b_s^e$, namely the scattering sources of the $u^e$-subproblem and the adjoint wavefields of the $m$-subproblem, are small. The sources of the adjoint equation satisfied by $\delta b_s^e$ are the weighted data residuals of classical FWI by the inverse of the data domain Hessian, equation 2.22, that is $(1/\mu) \Delta d_s^e$ according to equation 2.53. The division of $\Delta d_s^e$ by the small $\mu$ restores the true strength of the scattering source. Similarly, the adjoint source of ES-FWI in the gradient of the ES-FWI misfit function is $H_d^{-1}(m)^{-1}\delta d^*$ and not $\Delta d^e$ as highlighted by the generalized form of the ES-FWI misfit function, equation 2.46. This can be further checked by computing explicitly the Gauss-Newton descent direction of the $m$-subproblem for $u^e$ fixed and independent of $m$. The gradient of the penalty objective function, equation 2.7, is given by

$$ \nabla_m C^P(\mathbf{u}_s^e, \mathbf{m})|_{\mathbf{u}_s^e} = \mu \sum_s \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right)^T \delta \mathbf{b}_s^e \quad (2.54) $$

$$ = \mu \sum_s \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right)^T \mathbf{S}^T(\mathbf{m}) \mathbf{H}_d^{-1}(\mathbf{m}) \delta \mathbf{d}_s^*(\mathbf{m}), \quad (2.55) $$

and the Gauss-Newton Hessian reads

$$ \nabla_m^2 C^P(\mathbf{u}_s^e, \mathbf{m})|_{\mathbf{u}_s^e} = \mu \sum_s \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right)^T \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right). \quad (2.56) $$

Hence, the Gauss-Newton descent direction $p_k$ is given by

$$ p_k(\mathbf{m})|_{\mathbf{u}_s^e} = (\nabla_m^2 C^P(\mathbf{u}_s^e, \mathbf{m})|_{\mathbf{u}_s^e})^{-1} \nabla_m C^P(\mathbf{u}_s^e, \mathbf{m})|_{\mathbf{u}_s^e} = \frac{\sum_s \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right)^T \mathbf{S}^T(\mathbf{m}) \mathbf{H}_d^{-1}(\mathbf{m}) \delta \mathbf{d}^*(\mathbf{m})}{\sum_s \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right)^T \left( \frac{\partial \mathbf{A}(\mathbf{m}) \mathbf{u}_s^e}{\partial \mathbf{m}} \right)}. \quad (2.57) $$

Note how $\mu$ is removed by division of the gradient by the Hessian. Remembering the expressions of $u_s^e$, equation 2.28, and $\delta b_s^e$, equation 2.22, one can readily see that the penalty parameter $\mu$ is involved in the descent direction of the $m$-subproblem only as a damping regularization term of the data-domain Hessian of the scattering source estimation problem consistently with the expression of the generalized FWI misfit function, equation 2.46.

<figure><img src="image_10.png" /><figcaption>Figure 8: (a) Recorded data. (b) Simulated data in the reduced parameter space, $d(\mathbf{m}) = S(\mathbf{m})b$. (c) Simulated data in the extended space, equation 2.52, as the convex combination of the recorded data and simulated data. (d) Data residuals in the reduced parameter space (the recorded scattered data). (e) Data residuals in the extended space.</figcaption></figure>
## 2.10 On the role of the data-domain Hessian in the adjoint source

It remains to us to discuss the role of the data-domain Hessian in the adjoint source of ES-FWI in the light of the former analysis of the data-assimilated wavefields $u_s^e$. An heuristic interpretation is given hereafter. Let's first remind that this data-domain Hessian is the Hessian of a source location problem. Accordingly, let's imagine that we want to estimate the source $b$ either from the pair $(d^*, m^*)$ or $(d(m), m)$ by solving the damped least-squares data-fitting problems

$$ \min_{b} \|S(m)b - d(m)\|_2^2 + \mu \|b\|_2^2, \qquad (2.58) $$

$$ \min_{b} \|S(m^*)b - d^*\|_2^2 + \mu \|b\|_2^2. \qquad (2.59) $$

The minimizers are

$$ b = S^T(m) (S(m)S^T(m) + \mu I)^{-1} d(m) \qquad (2.60) $$

$$ b = S^T(m^*) (S(m^*)S^T(m^*) + \mu I)^{-1} d^*. \qquad (2.61) $$

We can expect the two reconstructed sources to be similar as well as the effect of the data-domain Hessian on $d(m)$ and $d^*$ because the data are kinematically consistent with the model in both cases although a more focused reconstruction of $b$ is expected from that model providing the best illumination of its wavenumber spectrum (for example, if $m^*$ contains reflectors as opposed to a smooth $m$, $b$ may be better reconstructed from $m^*$ due to the improved wavenumber illumination provided by back-scattered waves). Now, let's imagine that we wish to locate $b$ from the recorded data and the background $m$

$$ \min_{b} \|S(m)b - d^*\|_2^2 + \mu \|b\|_2^2. \qquad (2.62) $$

The reconstructed source will be smeared in space and time due to the kinematic inconsistency between $d^*$ and $m$ and the focusing effect of the inverse of the data-domain Hessian $S(m)S^T(m)$ on $d^*$ may be less effective accordingly.

<figure><img src="image_11.png" /><figcaption>Figure 9: Effect of the data-domain Hessian on the data residuals. (a) Data residuals $\delta d^*(m)$ for the velocity models of Figure 5. (b) Data residuals after deconvolution by $H_d^{-1}(m)$. (c) Recorded data $d^*$.</figcaption></figure>

The contrasted effect of the data-domain Hessian on $d(m)$ and $d^*$ is shown in Figure 9 for the velocity models shown in Figure 5. The direct arrival in the background model has been focused at short offsets and near onset times such that the adjoint operator $S^T(m)$ in equation 2.22 optimally focuses the back-propagated simulated data at the position of $b$ and at onset times. Conversely, the data-domain Hessian has less visible effects on the recorded data. Therefore, we may conclude that the descent direction of the ES-FWI iterations is mostly driven by $d^*$ at the expense of $d(m)$ when $m$ is strongly inaccurate due to the kinematic consistency of the incident wavefields $u_s^e$ in the virtual sources of the simulated scattered data and the back-propagated recorded data, while the opposite behavior occurs in classical FWI due to the kinematic consistency between the incident background wavefield and the back-propagated simulated data. The action of the data-domain Hessian in the adjoint source as a damping of the simulated component of the data residuals accounts for the fact that the true scattering sources are modeled approximately in the sensitivity matrix
of ES-FWI as highlighted by the resolution matrix $R(m)$, equation 2.22-2.23. Accordingly, the simulated scattered data $\delta d^e$ (and the partial derivative data) cannot predict well enough the recorded scattered data to capture the full differential information they carry out as highlighted by the data resolution matrix $R_d(m)$ in equation 2.24. This is accounted for by weighting the data residuals in the adjoint source by the inverse of the Hessian of the scattering source location problem. This contrasts with the adjoint source of classical FWI where the data residuals are not weighted. This results because the weak-scattering Born approximation implicitly assumes that the simulated scattered data in the sensitivity matrix can capture the full differential information contained in the data residuals generated by weak scattering. Therefore, as ES-FWI converges toward an accurate $m$, the scattering generated by $\delta m^*$ becomes weaker and weaker, the effect of the data-domain Hessian on the simulated and recorded data becomes similar, and ES-FWI tends asymptotically toward FWI.

# 3 Few remarks on algorithmic aspects

Discussing in depth algorithmic aspects of ES-FWI is beyond the scope of this tutorial, which is focused on the physical principles governing ES-FWI. Here, we only briefly review for sake of completeness the algorithmic strategies that were addressed in previous studies or that will deserve further investigations to design workable ES-FWI algorithms with the best trade-off between computational efficiency and solution accuracy.

## 3.1 Optimization algorithm and regularization

### 3.1.1 Why Augmented Lagrangian method?

As mentioned in the section *FWI as a constrained optimization problem: full space versus reduced space methods*, ES-FWI can be implemented with penalty method or augmented Lagrangian method. It is well acknowledged that one issue of penalty methods is related to the adaptive tuning of the penalty parameter such that a sufficient constraint relaxation is generated during the early iterations while the constraint is satisfied at the convergence point (Fu & Symes, 2017). Augmented Lagrangian method allows for a constant penalty parameter to be used because the Lagrange multipliers will record the history of the constraint errors in iterations to progressively remove their footprint in iterations following an iterative refinement procedure (i.e., the iterative solution of an ill-posed linear problem) (Aghamiry et al., 2019c, Their Appendix A). As an illustration of this recording, one can readily check from equation 2.8 that the Lagrange multipliers $v_s$, also referred to as the dual variables, reduce to the running sum of the constraint violations in iterations, when [i] the primal variables, namely $u_s$ and $m$, and the dual variables $v_s$ are updated in alternating mode, and [ii] $v_s$ are updated with basic gradient ascent steps. The reader is referred to Nocedal & Wright (2006, Chapter 17) for a comparative analysis of penalty methods and augmented Lagrangian methods from a mathematical perspective, and Aghamiry et al. (2019c) for their assessment in the context of ES-FWI.

### 3.1.2 Implementing regularization with the alternating-direction method of multipliers (ADMM)

ES-FWI is an underdetermined optimization problem as highlighted by the ill-posedness of the scattering source estimation problem, which requires some regularization. Another advantage of the augmented Lagrangian method (or method of multipliers) is to provide a versatile framework to implement various kinds of regularization including nonsmooth (non differential) ones with the alternating-direction method of multipliers (ADMM) (Esser et al., 2018; Goldstein & Osher, 2009). In this approach, some auxiliary variables are introduced in the regularized optimization problem via an additional constraint to decouple the least-squares problem from the $l_1$-norm problem and recast the later as a denoising problem, which can be solved efficiently with proximal algorithms (Parikh & Boyd, 2013). The resulting constrained problem is solved with ADMM by updating the various classes of primal variables in alternating mode. Moreover, the primal and the dual variables are also updated in alternating mode in the frame of the method of multipliers. The reader is referred to Aghamiry et al. (2019b, 2020b, 2021b) for the implementation of various kinds of regularization in ES-FWI and FWI with ADMM.

## 3.2 ES-FWI in the time domain

ES-FWI was originally formulated in the frequency domain because it was unclear how to compute the data-assimilated wavefields with explicit time stepping schemes since these wavefields are the solution of an overdetermined linear system, (e.g. van Leeuwen & Herrmann, 2013, Their equation 6) and hence of a normal equation (e.g. Aghamiry et al., 2020a, Their equation 4). This tutorial has reviewed how to tackle this issue by solving the wave equation in the background model $m$ with an extended source in the right hand side formed by the physical source plus a scattering source, equation 2.28 and equations 2.25-2.26. Another issue is related to the management of the Lagrange multipliers in the time domain since they have the size of the wavefields. Gholami et al. (2022) showed how
to project these Lagrange multipliers in the data domain through the relationship between $\delta\mathbf{b}^e$ and $\mathbf{H}_d^{-1}\delta\mathbf{d}^*$ via the adjoint of the wave-equation operator. Finally, compared to penalty methods, implementing ES-FWI with the augmented Lagrangian method just requires adding to the weighted data residuals $\mathbf{H}_d^{-1}\delta\mathbf{d}^*$ their running sum in iterations in the adjoint source of ES-FWI.

### 3.3 Inversion preconditioning

We have shown how the accuracy of the reconstructed wavefields decreases away from the receivers. It remains unclear under which conditions this varying accuracy will drive ES-FWI toward accurate or spurious subsurface models. To illustrate this point, inaccurate reconstructed wavefields sounds inconsistent with the linearization of the *m*-subproblem around these wavefields. Therefore, one may investigate depth continuation or layer stripping approaches to further extend the linear regime of ES-FWI by updating first the shallow part of the subsurface before considering the deeper part. This might be implemented with appropriate covariance matrices in the data domain and/or in the source domain (Lee & Pyun, 2020).

### 3.4 Waveform inversion workflow

We have shown how ES-FWI evolves toward FWI as the inversion approaches the convergence point. Accordingly, one may view to degrade the accuracy with which some ingredients of ES-FWI are estimated during the late iterations to mitigate its computational burden or switch to FWI during the latest stages of the inversion. On this subject, the data-domain Hessian of the scattering-source reconstruction subproblem generates the main computational overhead compared to classical FWI when its effects are taken into account accurately. Therefore, selecting the most suitable low-rank approximation of this data domain Hessian in iterations (diagonal approximation (Gholami et al., 2022), matching filters (Guo et al., 2022), truncated Gauss-Newton approximation (Guo et al., 2022), point-spread function (Lin et al., 2022), receiver encoding (Aghamiry et al., 2022)) to find the best compromise between computational cost and quality of the solution is one of the main issue of ES-FWI.

## 4 Numerical examples

We conclude this tutorial with two numerical illustrations corresponding to the Camembert model and the 2004 BP salt model. We would like to stress that the results obtained hereafter have been obtained with the augmented-Lagrangian implementation of ES-FWI. It would be very challenging if not impossible to obtain such results with the penalty method implementation even when adaptive tuning of the penalty parameter is used. Before showing the results of this test, we briefly discuss with a numerical example the convexity of the ES-FWI objective function.

### 4.1 On the convexity of the ES-FWI objective function

To the best of our knowledge, there is no mathematical proof of the convexity of the ES-FWI objective function. Anyway, it is easy to design a numerical experiment that would trap ES-FWI in a spurious minimum suggesting that this convexity cannot be proved in situations met in exploration geophysics. The best we can propose is to illustrate with numerical examples that the objective function of ES-FWI is more convex than that of classical FWI. To do so, we consider a 2D laterally homogeneous model of horizontal and vertical dimensions 8 km × 3 km where the velocity linearly increases with depth, $v(z) = v_0 + \alpha z$ (Mulder & Plessix, 2008). The acquisition consists of a single source located at x = 200 m, and a line of receivers spaced 50 m apart along the surface. The data are generated with a velocity $v0 = 2000$ m/s and a gradient $\alpha = 0.65 s^{-1}$. The cost function of FWI and ES-FWI when $\mu = 1e-2$ are computed in a 41 × 41 grid covering the 2D search space spanned by the two parameters $v0$ and $\alpha$: $v0 \in [1750, 2250]$ and $\alpha \in [0.4, 0.9]$ with steps of 12.5 m/s and 0.0125 1/s for $v0$ and $\alpha$, respectively. The FWI and ES-FWI cost functions are shown in Figure 10. As expected the FWI cost function shows several local minima unlike the ES-FWI cost function. A similar test is shown in van Leeuwen & Herrmann (2013) where $v0$ and $\alpha$ are processed separately.

### 4.2 The Camembert test

The Camembert model (Figure 11) contains a large circular inclusion in an homogeneous background (Gauthier et al., 1986). The wavespeeds in the homogeneous background and in the inclusion are $V_0 = 3.2$ km/s and $V_{inc} = 3.5$ km/s, respectively, and the radius of the inclusion is 3.5 km. A circular acquisition formed by 60 sources and receivers surrounds the inclusion to provide a full scattering-angle illumination of the target. Compared to the experimental setup designed by Gauthier et al. (1986), we use an ideal acquisition to focus on the resilience of
<figure><img src="image_12.png" /><figcaption>Figure 10: (a) FWI and (b) ES-FWI objective functions for the velocity gradient test.</figcaption></figure>

ES-FWI to cycle skipping without acquisition footprint generated by incomplete illumination. We perform a mono-frequency inversion for the 5 Hz frequency to precisely quantify the amount of cycle skipping. Moreover, we apply total-variation regularization to manage the ill-posedness of the scattering-source estimation. The starting model is the homogeneous background model ($V_0 = 3.2$ km/s). Ten wavelengths are propagated in the inclusion ($N_\lambda = 10$) and the strength of the contrast $\chi = (V_{inc} - V_0)/V_0$ is 20%. The traveltime mismatch $\delta t$ accumulated during the propagation across the inclusion is 0.1875 s, which is almost two times greater than half the period (0.1 s). The path followed by ES-FWI is illustrated in Figure 11 at iterations 10, 30, 100. The top row illustrates, for the left most source, the true wavefield (left panel) and the simulated wavefields in the estimated models $\mathbf{m}_k^e$ by ES-FWI, namely $\mathbf{u}(\mathbf{m}_k^e) = \mathbf{A}^{-1}(\mathbf{m}_k^e)\mathbf{b}$. We show these wavefields in the equivalent form of their exact scattering source $\delta\mathbf{b}(\mathbf{m}; \mathbf{m}_b) = -\omega^2\mathrm{diag}(\mathbf{m} - \mathbf{m}_b)\mathbf{u}(\mathbf{m})$ when the background model $\mathbf{m}_b$ is set to $\mathbf{m}_0$, namely

$$ \mathbf{u}^* = \mathbf{A}^{-1}(\mathbf{m}^*)\mathbf{b} = \mathbf{A}^{-1}(\mathbf{m}_0)(\mathbf{b} + \delta\mathbf{b}(\mathbf{m}^*; \mathbf{m}_0)), \quad (4.1) $$

$$ \mathbf{u}(\mathbf{m}_k^e) = \mathbf{A}^{-1}(\mathbf{m}_k^e)\mathbf{b} = \mathbf{A}^{-1}(\mathbf{m}_0)(\mathbf{b} + \delta\mathbf{b}(\mathbf{m}_k^e; \mathbf{m}_0)). \quad (4.2) $$

The exact scattering source depends on the full medium and the background medium and should not be confused with the estimated scattering source from the background medium and the recorded data $\mathbf{d}^*$, equation 2.22. Accordingly, $\mathbf{u}(\mathbf{m}_k^e)$ should not be confused with the estimated $\mathbf{u}_k^e$ by ES-FWI from $\mathbf{m}_{k-1}^e$ and $\mathbf{d}^*$, equation 2.28. The representation of the reconstructed wavefields by their scattering source provides a more direct assessment of their accuracy as manifested by the smearing of $\delta\mathbf{b}(\mathbf{m}_k^e; \mathbf{m}_0)$ compared to $\delta\mathbf{b}(\mathbf{m}^*; \mathbf{m}_0)$, the spatial support of the later matching that of the true inclusion (left panel). We see that $\delta\mathbf{b}(\mathbf{m}_{100}^e; \mathbf{m}_0)$ matches fairly well $\delta\mathbf{b}(\mathbf{m}^*; \mathbf{m}_0)$ at the convergence point, which correlates with the good match between $\mathbf{m}_{100}^e$ and $\mathbf{m}^*$ (second row) according to the bilinearity of the wave equation. We also see that $\delta\mathbf{b}(\mathbf{m}_{10}^e; \mathbf{m}_0)$ still lacks the interior part of the true scattering source, which corresponds to the area away from the receivers. This correlates with the fact that $\mathbf{m}_{10}^e$ mostly contains the boundary of the inclusion (second row). The scenario where the short wavelengths are reconstructed before the longer ones (as opposed to the multiscale approach from long to short wavelengths implemented in classical FWI by frequency, time and offset continuations (e.g. Górszczyk et al., 2017)) generally traps FWI in spurious minima (Gauthier et al., 1986, Their Figure 9). We note however that $\delta\mathbf{b}(\mathbf{m}_{10}^e; \mathbf{m}_0)$ matches $\delta\mathbf{b}(\mathbf{m}^*; \mathbf{m}_0)$ also at receivers located far away from the source, hence highlighting the resilience of ES-FWI to cycle skipping. The reconstructed models in iterations (second row) show how the inclusion is progressively reconstructed from its boundary (near the receivers) to its center (far away from the receivers). Transposing this statement to surface acquisition would lead to the conclusion that ES-FWI reconstructs the subsurface from the shallow parts to the deeper parts as the accuracy of the wavefields improved away from receivers in iterations following some kinds of layer stripping. The two bottom rows show the sought scattering sources $\delta\mathbf{b}(\mathbf{m}^*; \mathbf{m})$ and the estimated scattering sources $\delta\mathbf{b}^e(\mathbf{m})$ by ES-FWI for $\mathbf{m} = \mathbf{m}^*$, $\mathbf{m}_{k=10,30,100}^e$. Indeed,
the left panels ($\mathbf{m} = \mathbf{m}^*$) are nil since the contrast are zero and hence the scattering source and the data residuals are zero too. The third row shows how $\delta\mathbf{b}(\mathbf{m}^*; \mathbf{m})$ vanishes as $\mathbf{m}$ tends to $\mathbf{m}^*$. Moreover, comparing the third and fourth rows highlights how the relative error between $\delta\mathbf{b}(\mathbf{m}^*; \mathbf{m}_k^e)$ and $\delta\mathbf{b}^e(\mathbf{m}_k^e)$ decreases as $\mathbf{m}_k^e$ tends to $\mathbf{m}^*$, i.e., as $\mathbf{u}_k^e$ tends to $\mathbf{u}^*$. This trend highlights how the ill-posedness of the scattering source reconstruction decreases as the estimated wavefield becomes closer to the true one, remembering that these scattering sources depend on two classes of unknowns, namely the sought contrast and the sought wavefield. To summarize, the two-step workflow of ES-FWI (section *Summary of the overall ES-FWI workflow*) is illustrated by noting how the reconstructed model and the wavefield reconstruction improve at the same time while the reconstructed scattering source vanishes. The central issue of ES-FWI is to mitigate the ill-posedness of the scattering source estimation problem due to the potential trade-off between the wavefields and the model perturbations.

<figure><img src="image_13.png" /><figcaption>Figure 11: Camembert test. Top row: from left to right, $\delta\mathbf{b}(\mathbf{m}; \mathbf{m}_0)$ for $\mathbf{m} = \mathbf{m}^*$, $\mathbf{m}_{k=10,30,100}^e$ as an alternative representation of the ground-truth wavefield $\mathbf{u}^*$ and the simulated wavefields in $\mathbf{m}_{k=10,30,100}^e$. Second row: From left to right, true model $\mathbf{m}^*$ and ES-FWI models $\mathbf{m}_{k=10,30,100}^e$. Third and fourth rows: Sought scattering sources $\delta\mathbf{b}(\mathbf{m}^*; \mathbf{m})$ and estimated scattering source $\delta\mathbf{b}^e(\mathbf{m})$ by ES-FWI, equation 2.22, for $\mathbf{m} = \mathbf{m}^*$, $\mathbf{m}_{10}$, $\mathbf{m}_{30}$, $\mathbf{m}_{100}$.</figcaption></figure>
4.3 The 2004 BP salt benchmark

We conclude with the challenging 2004 BP salt model to illustrate the potential of ES-FWI for imaging complex subsurface models from inaccurate (smooth) starting model and ultra long-offset data when salt bodies generate a quite uneven illumination of the subsalt structure (Figure 12a).

<figure><img src="image_14.png" /><figcaption>Figure 12: 2004 BP salt model benchmark. (a) Ground-truth model **m**\*. (b) Starting model **m**₀. (c) Final ES-FWI model **m** ES-FWI.</figcaption></figure>

The dimensions of the model are 67 km × 12 km. The ultra long-offset stationary-recording acquisition consists of 67 hydrophones spaced 1 km apart on the seabed, and 450 pressure sources spaced 150 m apart at 25-m depth. For the sake of computational efficiency, we use the spatial reciprocity of Green’s functions to process sources as receivers and vice versa. A free-surface boundary condition is used on top of the grid, hence the surface multiples are involved in ES-FWI, which improves illumination while increasing nonlinearity due to more complex data anatomy (Figure 13). The source signature is a 4 Hz Ricker wavelet. The starting model **m**₀ is a crude velocity-gradient model (Figure 12b). Direct comparison between the simulated seismograms in the ground-truth model **m**\* and **m**₀, namely d\* = S(**m**\*)b and d(**m**₀) = S(**m**₀)b, respectively, shows obvious cycle skipping (Figure 13a), while direct comparison between d\* and the simulated data-assimilated seismograms in **m**₀, namely dᵉ(**m**₀) = S(**m**₀) (b + δbᵉ(**m**₀)), highlights the ability to match arbitrarily well the recorded data provided that the normal system for δdᵉ, equation 2.25, is solved accurately (Figure 13b). Note that we compute dᵉ in the time domain (Figure 13b) by solving the normal system "exactly" in the frequency domain with a direct solver for each frequency sampling the source bandwidth before taking the inverse Fourier transform of the monochromatic data.

We perform ES-FWI in the frequency domain in the 1.5 Hz - 5 Hz frequency band with a frequency interval of 0.25 Hz. Following a multiscale frequency continuation strategy, we proceed over small frequency batches from low frequencies to higher ones when each batch contains two frequencies with one frequency overlap. We perform three paths through the frequencies, using the final model of one path as the initial model of the next one. The starting and finishing frequencies of the three paths are [1.5, 2], [1.5, 4], [3, 5] Hz, respectively. The stopping criterion of iterations is a maximum of ten iterations per batch in all cases. We stabilize inversion with bound constraints and compound regularization formulated as the infimal convolution of TV and Tikhonov regularizations (Aghamiry et al., 2020b).

The final ES-FWI model denoted by **m**<sub>*f*</sub><sup>*e*</sup> matches fairly well the ground-truth model with some mild ringing artifacts (Figure 12a,c). Accordingly, the seismograms computed in this final model d(**m**<sub>*f*</sub><sup>*e*</sup>) = S(**m**<sub>*f*</sub><sup>*e*</sup>)b match fairly well the recorded counterparts d<sup>*</sup> (Figure 13c). Since ES-FWI aims at fitting both the data and the source, it is also worth
checking how the wave equation constraint is relaxed at the beginning of the inversion and how well it is satisfied at the convergence point. Example of scattering sources at the first and final iterations of the first frequency-batch inversion, namely $\sum_s \delta b_s^e(\mathbf{m}_{0,1.5}^e)$ and $\delta b_s^e(\mathbf{m}_{10,1.5}^e)$ respectively, are shown in Figure 14 where the second subscript denotes the frequency. Matching the data with the highly-inaccurate model shown in Figure 12b generates strong scattering sources with large spatial support covering the full subsurface domain (Figure 14a). At iteration 10, we show how these scattering sources vanish as the multivariate inversion for wavefield and parameters approaches the convergence point, equations 2.44-2.45. The fact the scattering sources vanish means that both the physical source and the data were matched since the scattering sources are computed by back-propagating the data residuals. Another example is shown in Figure 15 for the 5 Hz frequency and the first frequency path. We show also the true scattering source $\delta b(\mathbf{m}^*, \mathbf{m}_{0,5})$ in Figure 15b for sake of comparison with $\delta b_s^e(\mathbf{m}_{0,5})$ even if it cannot be matched since $\mathbf{m}^*$ is beyond the resolution power of FWI at 5 Hz. The fact that $\delta b(\mathbf{m}^*, \mathbf{m}_{0,5})$ mainly contains the edges of the salt bodies suggests that their long wavelengths have been already reconstructed at previous frequencies. Also, the amplitudes of $\delta b_s^e(\mathbf{m}_{0,5})$ (Figure 15a) look smaller than those of the 1.5 -Hz scattering source (Figure 14a) suggesting that the weak scattering approximation may be valid at this stage. Therefore, we may switch to FWI at this stage or use a low-rank approximation of the data-domain Hessian in ES-FWI. These numerical strategies will be discussed in more details in future studies. Finally, we see that $\delta b_s^e(\mathbf{m}_{10,5})$ vanishes at the convergence point (Figure 15c).

## 5 Conclusion

ES-FWI modifies the gradient of FWI with more accurate wavefields in the virtual sources of the partial derivative data. Beyond the Born approximation, this allows the partial derivative data to account approximately for scattering generated by the sought model perturbation along the incident wavepath to the scatterers. These more accurate wavefields are formulated as the sum of the background wavefields plus an approximation of the scattered wavefields by the sought model perturbation. These scattered wavefields are computed by solving the wave equation in the background model with a scattering source in the right-hand side (namely, by solving the Lippmann-Schwinger equation). These scattering sources are estimated by solving an underdetermined scattered data fitting problem, where the scattered data are the data residuals at the current iteration. Accordingly, the scattering sources are computed by propagating backward in time the weighted data residuals by the inverse of the data-domain Hessian. When the wavefields and the model are updated in alternating mode, the latter is updated by minimizing the scattering sources in a least-squares sense, which amounts to minimize the data residuals while pushing the updated model toward the ground truth model. In the end, the imaging principle of FWI and ES-FWI are similar: they rely on the correlation between the simulated virtual sources and the estimated scattering sources from the observables (generally, referred to as the adjoint wavefields in the FWI framework). However, these two ingredients have different expressions in FWI and ES-FWI. The virtual sources are built with data-assimilated wavefields and background wavefields in ES-FWI and FWI, respectively. Moreover, the estimated scattering sources are the least-squares solutions of the scattered-data fitting problem in ES-FWI, while they are the adjoint approximation of this problem in FWI, a consequence of the weak-scattering approximation underlying FWI. These modifications allow ES-FWI to match the recorded data arbitrarily well hence avoiding cycle skipping. However, the recorded scattered data (i.e., the data residuals) cannot be matched accurately by the simulated counterpart (i.e., the partial derivative data) for crude initial models due to the limited resolution with which the scattering sources are estimated. In this case, the data-domain Hessian accounts for the limited resolution with which the simulated scattered data can match the recorded counterpart by reshaping the data residuals in the so-called adjoint source to give more weight to the recorded data at the expense of the simulated data in the weighted data residuals. That is, the descent direction of ES-FWI is mostly driven by the recorded data at the expense of the simulated ones for crude background models. Then, the full differential information contained in the data residuals is captured more accurately as the inversion converges toward the minimizer, the accuracy of the scattering source estimation improves, and the data-domain Hessian has a similar effect on the recorded and simulated data, accordingly. With these remarks, the conditions that are required to converge toward accurate models remain unclear. The accuracy of the reconstructed wavefields degrades away from the receivers as they are built by propagating backward in time the recorded data in the background model. The spatial support and the amplitudes of the scattering sources increase accordingly and it becomes more challenging to minimize their footprint. Therefore, it is likely that heuristic approaches based on layer stripping combined with regularizations are useful to guide the inversion toward accurate models when the inversion starts from crude initial models. The computational burden of ES-FWI generated by the data-domain Hessian of the scattering-source estimation problem is another issue. However, taking into account for the Hessian accurately is probably not necessary during the late stages of the inversion when the weak scattering approximation becomes valid. Therefore, designing numerical strategies providing low-rank approximation of the data-domain Hessian and tuning the accuracy with which this Hessian should be approximated in iterations is the aim of ongoing work. The resilience of the method to noise in the frame of real data applications still needs to be assessed.
# 6 Acknowledgments

This study has been performed in the frame of the WIND project supported by Petrobras, Shell, Total. The authors are grateful to the OPAL infrastructure from Observatoire de la Côte d'Azur (CRIMSON) for providing resources and support. This work was granted access to the HPC resources of IDRIS under the allocation 0596 made by GENCI.

# References

Abubakar, A., Hu, W., van den Berg, P. M., & Habashy, T. M., 2008. A finite-difference contrast source inversion method, *Inverse problems*, **24**(6), 1–17.

Abubakar, A., Hu, W., Habashy, T. M., & van den Berg, P. M., 2009. Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data, *Geophysics*, **74**(6), WCC47–WCC58.

Abubakar, A., Pan, G., Li, M., Zhang, L., Habashy, T. M., & van den Berg, P., 2011. Three-dimensional seismic full-waveform inversion using the finite-difference contrast source inversion method, *Geophysical prospecting*, **59**, 874–888.

Aghamiry, H., Gholami, A., & Operto, S., 2019a. ADMM-based multi-parameter wavefield reconstruction inversion in VTI acoustic media with TV regularization, *Geophysical Journal International*, **219**(2), 1316–1333.

Aghamiry, H., Gholami, A., & Operto, S., 2019b. Implementing bound constraints and total-variation regularization in extended full waveform inversion with the alternating direction method of multiplier: application to large contrast media, *Geophysical Journal International*, **218**(2), 855–872.

Aghamiry, H., Gholami, A., & Operto, S., 2019c. Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers, *Geophysics*, **84**(1), R139–R162.

Aghamiry, H., Gholami, A., & Operto, S., 2020a. Accurate and efficient wavefield reconstruction in the time domain, *Geophysics*, **85**(2), A7–A12.

Aghamiry, H., Gholami, A., & Operto, S., 2020b. Compound regularization of full-waveform inversion for imaging piecewise media, *IEEE Transactions on Geoscience and Remote Sensing*, **58**(2), 1192–1204.

Aghamiry, H., Gholami, A., & Operto, S., 2020c. Multi-parameter wavefield reconstruction inversion for wavespeed and attenuation with bound constraints and total variation regularization, *Geophysics*, **85**(4)(4), 1–16.

Aghamiry, H., Gholami, A., & Operto, S., 2020d. Robust wavefield inversion with phase retrieval, *Geophysical Journal International*, **221**(2), 1327–1340.

Aghamiry, H., Gholami, A., & Operto, S., 2021a. Complex-valued imaging with total variation regularization: An application to full-waveform inversion in visco-acoustic, *SIAM Journal on Imaging Sciences (SIIMS)*, **14**(1), 58–91.

Aghamiry, H., Gholami, A., & Operto, S., 2021b. Full Waveform Inversion by Proximal Newton Methods using Adaptive Regularization, *Geophysical Journal International*, **224**(1), 169–180.

Aghamiry, H., Gholami, A., & Operto, S., 2021c. On efficient frequency-domain full-waveform inversion with extended search space, *Geophysics*, **86**(2)(2), R237.

Aghamiry, H., Gholami, A., & Operto, S., 2022. Highly-accurate wavefield reconstruction inversion using convergent born series, in *Expanded Abstracts*, EAGE.

Aghamiry, H. S., Mamfoumbi-Ozoumet, F. W., Gholami, A., & Operto, S., 2021d. Efficient extended-search space full-waveform inversion with unknown source signatures, *Geophysical Journal International*, **227**(1), 257–274.

Aghazade, K., Gholami, A., Aghamiry, H., & Operto, S., 2022a. Randomized source sketching for full waveform inversion, *IEEE Transactions on Geoscience and Remote Sensing*, **60**, 1–12.

Aghazade, K., Gholami, A., Aghamiry, H., & Operto, S., 2022b. Anderson accelerated augmented lagrangian for extended waveform inversion, *Geophysics*, **87**, R79–R91.

Baek, H., Calandra, H., & Demanet, L., 2014. Velocity estimation via registration-guided least-squares inversion, *Geophysics*, **79**(2), R79–R89.

Baeten, G., de Maag, J. W., Plessix, R.-E., Klaassen, R., Qureshi, T., Kleemeyer, M., ten Kroode, F., & Rujie, Z., 2013. The use of low frequencies in a full-waveform inversion and impedance inversion land seismic case study, *Geophysical Prospecting*, **61**(4), 701–711.

Banerjee, B., Walsh, T. F., Aquino, W., & Bonnet, M., 2013. Large scale parameter estimation problems in frequency-domain elastodynamics using an error in constitutive equation functional, *Computational Methods in Applied Mechanics Engineering*, **253**, 60–72.
Barnier, G., Biondi, E., Clapp, R. G., & Biondi, B., 2012. Full waveform inversion by model extension: theory, design and optimization, arXiv:2205.14341v1.

Biondi, B. & Almomin, A., 2014. Simultaneous inversion of full data bandwidth by tomographic full-waveform inversion, *Geophysics*, **79**(3), WA129–WA140.

Brenders, A., Dellinger, J., Ahmed, I., Díaz, E., Gherasim, M., Jin, H., Vyas, M., & Naranjo, J., 2022. The wolfspar experience with low-frequency seismic source field data: Motivation, processing, and implications, *The Leading Edge*, **41**(1), 9–18.

eds Chavent, G. & Sabatier, P. C., 1996. *Inverse problems of wave propagation and diffraction*, Lecture Notes in Physics, Proceedings, Aix-les-Bains, France, Springer.

Engquist, B., Froese, B. D., & Yang, Y., 2016. Optimal transport for seismic full waveform inversion, *Communications in Mathematical Sciences*, **14**(8), 2309–2330.

Epanomeritakis, I., Akcelik, V., Ghattas, O., & Bielak, J., 2008. A Newton-CG method for large-scale three-dimensional elastic full waveform seismic inversion, *Inverse Problems*, **24**, 1–26.

Esser, E., Guasch, L., van Leeuwen, T., Aravkin, A. Y., & Herrmann, F. J., 2018. Total variation regularization strategies in Full-Waveform Inversion, *SIAM Journal Imaging Sciences*, **11**(1), 376–406.

Fang, Z., Wang, R., & Herrmann, F. J., 2018. Source estimation for wavefield-reconstruction inversion, *Geophysics*, **83**(4), R345–R359.

Fiddy, M. A. & Pommet, D. A., 1996. *Recovery of strongly scattering permitivity distributions from limited backscattered data using a nonlinear filtering technique*, vol. 15B, pp. 58–70, Lecture Notes in Physics, Springer.

Fu, L. & Symes, W. W., 2017. A discrepancy-based penalty method for extended waveform inversion, *Geophysics*, **R282-R298**, 78–82.

Gauthier, O., Virieux, J., & Tarantola, A., 1986. Two-dimensional nonlinear inversion of seismic waveforms: numerical results, *Geophysics*, **51**(7), 1387–1403.

Gholami, A., Aghamiry, H. S., & Operto, S., 2022. Extended full waveform inversion in the time domain by the augmented Lagrangian method, *Geophysics*, **87**(1), R63–R77.

Goldstein, T. & Osher, S., 2009. The split Bregman method for L1-regularized problems, *SIAM Journal on Imaging Sciences*, **2**(2), 323–343.

Golub, G. & Pereyra, V., 2003. Separable nonlinear least squares: the variable projection method and its applications, *Inverse problems*, **19**(2), R1.

Górszczyk, A., Operto, S., & Malinowski, M., 2017. Toward a robust workflow for deep crustal imaging by FWI of OBS data: The eastern nankai trough revisited, *Journal of Geophysical Research: Solid Earth*, **122**(6), 4601–4630.

Guo, G., Operto, S., Gholami, A., & Aghamiry, H. S., 2022. A practical implementation of data-space hessian in the time-domain formulation of source extended full-waveform inversion, in *Second International Meeting for Applied Geoscience & Energy*, pp. 757–761, Society of Exploration Geophysicists and American Association of Petroleum ....

Hajjaj, R. F., de Ridder, S. A., Livermore, P. W., & Ravasi, M., 2022. Wavefield reconstruction inversion modelling of marchenko focusing functions, arXiv preprint arXiv:2210.14570.

Huang, G., Nammour, R., & Symes, W. W., 2018a. Source-independent extended waveform inversion based on space-time source extension: Frequency-domain implementation, *Geophysics*, **83**(5), R449–R461.

Huang, G., Nammour, R., & Symes, W. W., 2018b. Volume source-based extended waveform inversion, *Geophysics*, **83**(5), R369–387.

Lee, D. & Pyun, S., 2020. Seismic full-waveform inversion using minimization of virtual scattering sources, *Geophysics*, **85**, R299–R311.

Lin, Y., Liu, H., Xing, L., & Lin, H., 2022. Time-domain wavefield reconstruction inversion solutions in the weighted full waveform inversion form, *IEEE Transactions on Geoscience and Remote Sensing*, **60**, 1–14.

Ma, Y. & Hale, D., 2013. Wave-equation reflection traveltime inversion with dynamic warping and full waveform inversion, *Geophysics*, **78**(6), R223–R233.

Métivier, L., Allain, A., Brossier, R., Mérigot, Q., Oudet, E., & Virieux, J., 2018. Optimal transport for mitigating cycle skipping in full waveform inversion: a graph space transform approach, *Geophysics*, **83**(5), R515–R540.

Mora, P. R., 1987. Nonlinear two-dimensional elastic inversion of multi-offset seismic data, *Geophysics*, **52**, 1211–1228.
Mulder, W. & Plessix, R. E., 2008. Exploring some issues in acoustic full waveform inversion, *Geophysical Prospect-ing*, **56**(6), 827–841.

Nocedal, J. & Wright, S. J., 2006. *Numerical Optimization*, Springer, 2nd edn.

Parikh, N. & Boyd, S., 2013. Proximal algorithms, *Foundations and Trends in Optimization*, **1**(3), 123–231.

Peters, B. & Herrmann, F. J., 2019. A numerical solver for least-squares sub-problems in 3d wavefield reconstruction inversion and related problem formulations, in *SEG International Exposition and Annual Meeting, OnePetro*.

Pratt, R. G., Shin, C., & Hicks, G. J., 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion, *Geophysical Journal International*, **133**, 341–362.

Prunty, A. C. & Snieder, R. K., 2020. An acoustic lippmann-schwinger inversion method: applications and comparison with the linear sampling method, *Journal of Physics Communications*, **4**(015007), 1–14.

Rizzuti, G., Louboutin, M., Wang, R., & Herrmann, F. J., 2021. A dual formulation of wavefield reconstruction inversion for large-scale seismic inversion, *Geophysics*, **86**(6), R879–R893.

Shin, C., Jang, S., & Min, D. J., 2001. Improved amplitude preservation for prestack depth migration by inverse scattering theory, *Geophysical Prospecting*, **49**, 592–606.

Sirgue, L., Barkved, O. I., Dellinger, J., Etgen, J., Albertin, U., & Kommedal, J. H., 2010. Full waveform inversion: the next leap forward in imaging at Valhall, *First Break*, **28**, 65–70.

Symes, W. W., 2008. Migration velocity analysis and waveform inversion, *Geophysical Prospecting*, **56**, 765–790.

Symes, W. W., 2020. Wavefield reconstruction inversion: an example, *Inverse Problems*, **36**(10), 105010.

Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation, *Geophysics*, **49**(8), 1259–1266.

Tarantola, A. & Valette, B., 1982. Generalized nonlinear inverse problems solved using the least square criterion, *Reviews of Geophysical and Space Physics*, **20**, 219–232.

van den Berg, P. M. & Kleinman, R. E., 1997. A contrast source inversion method, *Inverse Problems*, **13**(6), 1607.

van Leeuwen, T., 2019. A note on extended full waveform inversion, arXiv preprint arXiv:1904.00363.

van Leeuwen, T. & Herrmann, F., 2016. A penalty method for PDE-constrained optimization in inverse problems, *Inverse Problems*, **32**(1), 1–26.

van Leeuwen, T. & Herrmann, F. J., 2013. Mitigating local minima in full-waveform inversion by expanding the search space, *Geophysical Journal International*, **195**(1), 661–667.

van Leeuwen, T., Herrmann, F., & Peters, B., 2014. A new take on FWI - wavefield reconstruction inversion, in *76<sup>th</sup> Annual EAGE Meeting (Amsterdam)*.

Virieux, J. & Operto, S., 2009. An overview of full waveform inversion in exploration geophysics, *Geophysics*, **74**(6), WCC1–WCC26.

Wang, C., Yingst, D., Farmer, P., & Leveille, J., 2016. Full-waveform inversion with the reconstructed wavefield method, in *SEG Technical Program Expanded Abstracts*, pp. 1237–1241.

Wang, C., Yingst, D., Farmer, P., Jones, I., Martin, G., & Leveille, J., 2017. Reconstructed full-waveform inversion with the extended source, in *SEG Technical Program Expanded Abstracts*, pp. 1449–1453.

Warner, M. & Guasch, L., 2016. Adaptive waveform inversion: Theory, *Geophysics*, **81**(6), R429–R445.

Yang, Y., Engquist, B., Sun, J., & Hamfeldt, B. F., 2018. Application of optimal transport and the quadratic Wasserstein metric to full-waveform inversion, *Geophysics*, **83**(1), R43–R62.

Yao, G., da Silva, N. V., Warner, M., Wu, D., & Yang, C., 2019. Tackling cycle skipping in full-waveform inversion with intermediate data, *Geophysics*, **84**(3), R411–R427.
<figure><img src="image_15.png" /><figcaption>Figure 13: 2004 BP salt model benchmark. (a-c) Direct comparison between (a) $d^* = S(m^*)b$ and $d(m_0) = S(m_0)b$; (b) $d^*$ and $d^e(m_0)$; (c) $d^*$ and $d(m_f^e) = S(m_f^e)b$, where $m_f^e$ denotes the final model found by ES-FWI (Figure 12). In (a-c), $d(m_0)$, $d^e(m_0)$ and $d(m_f^e)$ are superimposed in transparency on $d^*$. The seismograms $d(m_0)$, $d^e(m_0)$ and $d(m_f^e)$ are plotted with black-white-yellow color scale (right scale) while $d^*$ are plotted with a red-white-blue color scale (middle scale). The two sets of seismograms match when their superimposition appears as brown-white-olive colored seismograms (left scale), because black+red=brown and yellow+blue=olive.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 14: 2004 BP salt model benchmark. (a-b) Estimated scattering source during the first frequency batch inversion of the first path (a) iteration 1. (b) iteration 10. The physical source is located at (x,z)=(7500,1350) m.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: 2004 BP salt model benchmark. (a,c) Estimated scattering source during the last frequency batch inversion (5 Hz) of the first path. (a) Iteration 1. (c) Iteration 10. (b) Exact scattering source, $\delta\mathbf{b}(\mathbf{m}^*, \mathbf{m}_{0,5})$ where $\mathbf{m}_{0,5}$ denotes the starting model of the frequency batch inversion involving the 5 Hz frequency. The physical source is located at (x,z)=(7500,1350) m.</figcaption></figure>