Robust Wavefield Inversion via Phase Retrieval

TL;DR

This paper introduces Wavefield Inversion with Phase Retrieval (WIPR), a robust method for seismic imaging that enhances traditional wavefield inversion by incorporating phase correction, improving accuracy in challenging initial models.

Contribution

The paper proposes a novel phase retrieval-based cost function for wavefield inversion, combined with an alternating-direction algorithm, to improve robustness against phase errors and local minima.

Findings

Successfully reconstructs the 2004 BP salt model from sparse data.

Demonstrates robustness with a simple initial velocity model.

Achieves accurate imaging at low frequencies (3 Hz).

Abstract

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

Figures

Tables

Table 1. Table 1: Model error (ME) for the velocity models and least-squares data and source misfit at the end of the first frequency batch inversion and at the convergence point of the four inversion tests (with or without phase retrieval, with or without TT regularization). Here, ME is defined as 100 ​ ‖ 𝕞 k − 𝕞 ∗ ‖ 1 ‖ 𝕞 ∗ ‖ 1 100 subscript norm subscript 𝕞 𝑘 subscript 𝕞 1 subscript norm subscript 𝕞 1 100\frac{\|\mathbb{m}_{k}-\mathbb{m}_{*}\|_{1}}{\|\mathbb{m}_{*}\|_{1}} , where 𝕞 ∗ subscript 𝕞 \mathbb{m}_{*} denotes the true model.

			ME	${\Vert 𝐀\left(𝐦\right)𝐮-𝐛\Vert }_{\mathrm{𝟐}}$	${\Vert \mathrm{𝐏𝐮}-𝐝\Vert }_{\mathrm{𝟐}}$
Iteration	No.45	IR-WRI	23.23	0.5954	7.8550
		WIPR	13.97	0.0672	1.8450
		IR-WRI+TT	19.81	0.5640	4.4630
		WIPR+TT	9.49	0.04926	0.6706
Final	iteration	IR-WRI	7.57	0.2046	2.8911
		${\text{IR-WRI}}_{pr}$	6.50	0.0647	1.6009
		IR-WRI+TT	4.89	0.1272	1.6470
		${\text{IR-WRI}}_{pr}$+TT	3.16	0.0359	0.6350

Equations

$$\underset{x}{\min }f(x)=\underset{x}{\min }\frac{1}{2}\Vert |Lx|-|y|{\Vert }_2^2.$$

$$x^{k+1}=\underset{x}{\mathrm{argmin}}\frac{1}{2}\Vert \mathbb{L}x-|y|e^{j\angle {\mathbb{L}x}^k}{\Vert }_2^2.$$

$$\underset{m}{\min }\Vert \mathbb{P}u(m)-d{\Vert }_2^2,$$

$$u(m)={\mathbb{A}}^{-1}(m)b,$$

$$\mathbb{A}(m)=\Delta +{\omega }^2\mathbb{C}\text{diag}(m)\mathbb{B},$$

$$\underset{m}{\min }\Vert \mathbb{A}(m)u(m)-b{\Vert }_2^2$$

$$u(m)={\left(\lambda \mathbb{A}(m{)}^T\mathbb{A}(m)+{\mathbb{P}}^T\mathbb{P}\right)}^{-1}\left(\lambda \mathbb{A}(m{)}^{T}b+{\mathbb{P}}^{T}d\right).$$

$$m^{k+1}=\underset{m}{\mathrm{argmin}}\Vert \mathbb{A}(m)u^k-b{\Vert }_2^2,$$

$$m^{k+1}=\underset{m}{\min }\Vert \mathbb{L}(u^k)m-y(u^k){\Vert }_2^2,$$

$$\{\begin{array}{l}\mathbb{L}(u^k)={\omega }^2\mathbb{C}\text{diag}(\mathbb{B}{u}^k),\\ y(u^k)=b-\Delta u^k.\end{array}$$

$$b^{k+1}$$

$$d^{k+1}$$

$$m^{k+1}=\underset{m}{\mathrm{argmin}}\Vert |\mathbb{L}(u^k)m|-|y(u^k)|{\Vert }_2^2.$$

$$m^{k+1}=\underset{m}{\mathrm{argmin}}\Vert \mathbb{L}(u^k)m-|y(u^k)|e^{j\angle \mathbb{L}(u^k)m^k}{\Vert }_2^2.$$

$${\mathbb{A}(m}^k)u^k\approx b$$

$$\mathbb{L}(u^k)m^k\approx y(u^k),$$

$$\underset{\varphi }{\min }\Vert \mathbb{L}(u^k)m^k-|y(u^k)|e^{j\varphi }{\Vert }_2^2.$$

$$m^{k+1}=\underset{m\in \mathcal{C}}{\mathrm{argmin}}\Vert m{\Vert }_{\text{TT}}+\lambda \Vert \mathbb{L}(u^k)m-|y(u^k)|e^{j\angle \mathbb{L}(u^k)m^k}{\Vert }_2^2,$$

$$\Vert m{\Vert }_{\text{TT}}=\underset{m=m_1+m_2}{\min }\Vert m_1{\Vert }_{\text{TV}}+\alpha \Vert m_2{\Vert }_{\text{Tikh}},$$

$$\{\begin{array}{l}\Vert m{\Vert }_{\text{TV}}=\sum \sqrt{|{\nabla }_{x}m{|}^2+|{\nabla }_{z}m{|}^2},\\ \Vert m{\Vert }_{\text{Tikh}}=\sum \left(|{\nabla }_{xx}m{|}^2+2|{\nabla }_{xz}m{|}^2+|{\nabla }_{zz}m{|}^2\right),\end{array}$$

$$(\Vert {\mathbb{A}(m}^k)u^k-b{\Vert }_F\le 10^{-3}\text{and}\Vert {\mathbb{P}u}^k-d{\Vert }_F\le 10^{-5}),$$

$$\Delta u+{\omega }^2\text{diag}(m)u=b$$

$$\underset{\mathbb{L}(u)m}{\underbrace{{\omega }^2\text{diag}(m)u}}=\underset{y(u)}{\underbrace{b-\Delta u}}.$$

$$\Vert \mathbb{L}(u)m-y(u){\Vert }_2^2=\Vert (|\mathbb{L}(u)m|-|y(u)|)e^{j\angle \mathbb{L}(u)m}{\Vert }_2^2$$

$$=\Vert \mathbb{L}(u)m-|y(u)|e^{j\angle \mathbb{L}(u)m}{\Vert }_2^2.$$

$$m=\frac{1}{{\omega }^2}\left(\frac{1}{u}\right)\circ y(u),$$

$$m=\frac{1}{{\omega }^2}\left(\frac{1}{|u|}\right)\circ |y(u)|,$$

$$\underset{x}{\min }f(x)=\underset{x}{\min }\frac{1}{2}\Vert |Lx|-|y|{\Vert }_2^2,$$

$$\tilde{f}(z)=\frac{1}{2}\Vert \sqrt{z{z}^{\ast }}-|y|{\Vert }_2^2.$$

$$g(x,x^k)=f(x^k)+\nabla f(x^k{)}^T\Delta x+\Delta x^T\mathbb{H}(x^k)\Delta x,$$

$$\nabla f(x)$$

Full text

1 INTRODUCTION

Full waveform inversion (FWI) is a waveform matching procedure which provides a subsurface model with a resolution reaching to the half of the smallest propagated wavelength (Virieux & Operto, 2009). Ultra-long offset wide-azimuth stationary-recording acquisitions provide a varied angular illumination of the subsurface, which is amenable to the development of the broadband velocity model (Sirgue & Pratt, 2004). However, the large number of propagated wavelengths generated by these acquisitions makes FWI prone to cycle skipping or phase wrapping (Shah, 2014; Choi & Alkhalifah, 2015). Various continuation strategies in frequencies, travel-times and offsets can be used to mitigate this pathology (e.g., Shipp & Singh, 2002; Górszczyk et al., 2017). However, these multi-scale approaches can be tedious to implement and still require quite accurate initial velocity models and low frequencies to prevent cycle skipping at long offsets.

To increase the resilience of FWI to cycle skipping, frequency-domain wavefield reconstruction inversion (WRI) has been proposed by van Leeuwen & Herrmann (2013, 2016). The governing idea of the method relies on the following statement: if we would be able to record a monochromatic wavefield generated by a single source everywhere in the subsurface, then it would be straightforward to reconstruct the subsurface model exactly in virtue of the wave equation bilinearity (Fig. 1). Since we cannot record wavefields everywhere, the best we can do to approach the true wavefields is to find those that best match the sparse observations and satisfy the wave equation in a least-squares sense. This least-squares problem can be implemented by solving an overdetermined linear system gathering the weighted wave equation and the observation equation (the equation relating the data to the wavefield through a sparse sampling operator). The observation equation generates a wave equation relaxation whose strength depends on the accuracy of the available subsurface model and the weight assigned to the wave equation. Then, the subsurface parameters are updated from the reconstructed wavefields by the least-squares minimization of the wave equation errors (source residuals), and this two-step cycle is iterated until both the observation equation and the wave equation are satisfied with a prescribed accuracy. Typically, a small weight is assigned to the wave equation during early iterations to foster data fitting and prevent cycle skipping accordingly. This weight is progressively increased to guarantee that the wave equation is fulfilled near the convergence point. The parameter-estimation subproblem can be solved through a variable projection approach (van Leeuwen & Herrmann, 2016) or through alternating optimization (van Leeuwen & Herrmann, 2013). In this latter case, the parameter-estimation subproblem is linearized around the reconstructed wavefield according to the wave-equation bilinearity.

WRI has been originally implemented with a penalty method, which requires a tedious dynamic control of the penalty parameter as above mentioned (Fu & Symes, 2017). To overcome this issue, Aghamiry et al. (2019d) proposed to replace the penalty method with an augmented Lagrangian method equipped with operator splitting, namely the alternating-direction method of multiplier (ADMM) (Boyd et al., 2010). Contrary to penalty methods, augmented Lagrangian methods converge to an accurate solution with a fixed penalty parameter by iteratively updating the Lagrange multipliers with the running sum of the constraint violations (Nocedal & Wright, 2006, Chapter 17). This defines what is sometimes referred to as an iterative refinement or defect correction algorithm. Accordingly, Aghamiry et al. (2019d) called their method iteratively-refined (IR-) WRI.

ADMM also provides a suitable framework to implement bound constraints and $\ell 1$-based nonsmooth regularizations (such as total variation (TV) regularization) in IR-WRI (Aghamiry et al., 2019c). IR-WRI with TV regularization was further improved by using more versatile Tikhonov-TV (TT) regularization combining second-order Tikhonov and TV regularizers through infimal convolution (Aghamiry et al., 2019b). Infimal convolution means that the TT regularization explicitly decomposes the model into two components of different statistical properties (a smooth one and a blocky one) such that a suitable regularization can be tailored to each component (Tikhonov and TV) (Gholami & Hosseini, 2013).

Although IR-WRI, when equipped with basic frequency continuation strategies and advanced regularization techniques, mitigates cycle skipping and helps to overcome nonlinearity issues broadly speaking, still it can be trapped in a local minimum when the initial model is very far from the true model. In this case, the phase and amplitude of the reconstructed wavefields match those of the true wavefields only in the vicinity of the receivers, while they can be quite inaccurate elsewhere. In IR-WRI, the inaccuracies of the reconstructed wavefields are directly mapped into the subsurface model through a linear deconvolution-like processing. This implies that only the shallow part of the subsurface model is expected to be reliably updated during the early iterations of IR-WRI for surface acquisitions. It remains unclear from a theoretical viewpoint how these inaccuracies are progressively canceled out by IR-WRI in iterations through a depth-continuation process.

It is well acknowledged that the phase often has a more dominant role than amplitude in image processing and waveform inversion methods (Oppenheim & Lim, 1981; Shechtman et al., 2015). Therefore, the inaccuracies in the phase of the reconstructed wavefields may be those which can more likely drive the inversion toward spurious local minima. Accordingly, the objective of this study is to assess whether a reformulation of the parameter-estimation subproblem that mitigates the role of the phase during the early stages of the IR-WRI at the benefit of the more robust amplitude counterpart contributes to stabilizing the inversion.

To achieve this goal, we recast the parameter-estimation subproblem of IR-WRI as a phase retrieval problem (Fienup, 1982). Phase retrieval considers the fundamental problem of how to reconstruct a signal from the magnitude of linear measurements (e.g., the magnitude of its Fourier transform). The name phase retrieval arose because if the phase of the linear measurements can be retrieved, then it is easy to reconstruct the signal by solving a linear problem (Waldspurger et al., 2015). It has a long history of applications in science and engineering when the phase of the linear measurements cannot be recorded or is inaccurate., e.g., optics (Walther, 1963), X-ray crystallography (Millane, 1990; Harrison, 1993), astronomy (Fienup, 1982). Also, it has recently attracted renewed interest in many imaging problems (Fogel et al., 2016), such as X-ray medical imaging (Pfeiffer et al., 2006; Burvall et al., 2011), optical imaging (Shechtman et al., 2015) and seismic processing (Gholami, 2014).

In IR-WRI, the linear operator of the parameter estimation subproblem is formed by the so-called virtual sources (Pratt et al., 1998), while the right-hand sides (the linear measurements) depend on the source and the Laplacian of the wavefields (considering the Helmholtz equation as wave equation). According to the above definition of phase retrieval, we update the subsurface parameters from the magnitude of the right-hand sides rather than from their phase and amplitude. A potential downside of the phase-retrieval algorithm is ill-posedness generated by the rapidly-decreasing sensitivity of the inversion with depth, in particular in large contrast or attenuating media when a small amount of the seismic energy is transmitted in the deep subsurface. This makes the use of efficient regularization necessary, as we will show. As in IR-WRI, we implement the iterative refinement procedure in the phase retrieval algorithm to design an adaptive control on the regularization parameters and achieve a faster convergence to the desired model. Once the velocity model estimation has been stabilized by phase retrieval during a first low-frequency band inversion, we switch to classical IR-WRI to process the higher frequencies with the phase information.

We assess our phase-retrieval based waveform inversion, referred to as WIPR, against two complex synthetic case studies, the Marmousi II and the large contrast 2004 BP salt models with a sparse long-offset fixed-spread acquisition. We show how WIPR helps to reconstruct more accurately the phase of the reconstructed wavefields, which is translated into more accurate velocity estimation. Then, we show the impact of this more accurate early-stage velocity reconstruction on the convergence of IR-WRI at higher frequencies.

This paper is organized in a method, numerical, and discussion sections. In the method section, we first introduce phase retrieval in a general context. Second, we review the principles of WRI with iterative refinement and show how we implement phase retrieval in the underlying parameter-estimation subproblem. Then, we briefly review how compound regularizations can be implemented in WIPR. In the numerical experiment section, we show the application on the Marmousi and 2004 BP salt models. In the discussion section, we further discuss the resolution power of amplitudes in WIPR and show its impact on the phase reconstruction.

2 NOTATION

The mathematical symbols adopted in this paper are as follows. We use italics for scalar quantities, boldface lower-case letters for vectors, and boldface capital letters for matrices and tensors. We use the superscript $T$ to denote the adjoint of an operator and $*$ to show the complex conjugate. We use diag$(\bullet)$ to show a square diagonal matrix which includes vector $\bullet$ on its main diagonal. The $i$th component of the column vector $\mathbb{x}$ is shown by $\it{x}_{i}$. For a complex number $\it{x}=\it{a}e^{j\it{b}}$ with $j=\sqrt{-1}$, $|\it{x}|=\it{a}$ denotes the magnitude of $\it{x}$ and $\angle\it{x}=\it{b}$ denotes its phase. For the $n$-length column vectors $\mathbb{x}$ and $\mathbb{y}$, the dot product is defined by $\langle\mathbb{x},\mathbb{y}\rangle=\mathbb{x}^{T}\mathbb{y}=\sum_{i=1}^{n}\it{x}_{i}^{*}\it{y}_{i}$ and their Hadamard product, denoted by $\mathbb{x}\circ\mathbb{y}$, is another vector made up of their component-wise products, i.e. $(\mathbb{x}\circ\mathbb{y})_{i}=\it{x}_{i}\it{y}_{i}$. The $\ell_{1}$- and $\ell_{2}$-norms of $\mathbb{x}$ are, respectively, defined by $\|\mathbb{x}\|_{1}=\sum_{i=1}^{n}|\it{x}_{i}|$ and $\|\mathbb{x}\|_{2}=\sqrt{\langle\mathbb{x},\mathbb{x}\rangle}=\sqrt{\sum_{i=1}^{n}|\it{x}_{i}|^{2}}$.

3 METHOD

In this section, we briefly review the phase retrieval problem (Gerchberg, 1972) and a simple majorization-minimization (MM) (Lange, 2016) algorithm to solve this non-convex problem. Then we review the frequency-domain FWI (Pratt et al., 1998), its WRI alternative (van Leeuwen & Herrmann, 2016), and its improved version called iteratively refined WRI (IR-WRI) (Aghamiry et al., 2019d). Finally, we introduce phase retrieval in the parameter-estimation subproblem of IR-WRI. In the following of this study, we refer to IR-WRI with phase retrieval as WIPR. WIPR is solved efficiently with the alternating direction method of multipliers (ADMM) (Boyd et al., 2010) and MM algorithms, which provide a suitable framework to implement bound constraints and TT regularization in the parameter estimation subproblem and tackle large data sets.

3.1 Phase retrieval

Phase retrieval refers to the reconstruction of a complex signal from the amplitude of linear measurements. Consider the complex-valued linear problem $\mathbb{Lx}=\mathbb{y}$ with $\mathbb{L}\in\mathbb{C}^{m\times n}$, $\mathbb{x}\in\mathbb{C}^{n\times 1}$, and $\mathbb{y}\in\mathbb{C}^{m\times 1}$ and assume that only the amplitude of the right hand side, $|\mathbb{y}|$, can be measured or is reliable. The corresponding phase-retrieval problem can be written as (Gholami, 2014)

[TABLE]

Problem (3.1) is non-convex where this non-convexity finds its root in removing the phase of the right hand side ($|\mathbb{Lx}|=|\mathbb{y}|$ is non-unique and there are many $\mathbb{x}$ that can fit the magnitude of $\mathbb{y}$).

Several algorithms have been proposed to solve the problem ranging from the Gerchberg-Saxton algorithm (Gerchberg, 1972) as a first algorithm to more modern optimization techniques (Candes et al., 2013; Netrapalli et al., 2013; Eldar et al., 2016). In this paper, we present a simple algorithm based on the MM technique (Appendix References) to find a local minimum of $f(\mathbb{x})$ iteratively. Based on the proposed algorithm, the minimizer of the convex surrogate function

[TABLE]

converges to the minimizer of $f(\mathbb{x})$ as iteration number, $k$, tends to infinity. The quadratic problem 3.2 is the same as the least-squares norm of the original system $\mathbb{Lx=y}$, when the phase of $\mathbb{y}$ is replaced with the phase of $\mathbb{Lx}^{k}$, which is extracted from the solution of $k$’th iteration.

3.2 Full-waveform inversion versus wavefield inversion

The reduced formulation of frequency-domain FWI can be written as (Pratt et al., 1998; Plessix, 2006)

[TABLE]

where $\mathbb{d}$ denotes the recorded seismic data, $\mathbb{P}$ is the linear observation operator that samples the wavefield $\mathbb{u(m)}$ at the receiver positions, and $\mathbb{m}$ denotes the subsurface parameters. The wavefield $\mathbb{u}(\mathbb{m})$ is the solution of the wave equation

[TABLE]

where $\mathbb{b}$ is the source term and $\mathbb{A(m)}\in\mathbb{C}^{n\times n}$ is the discretized wave-equation operator. In this study, the wave equation is the Helmholtz equation, whose operator $\mathbb{A(m)}$ is given by

[TABLE]

where $\omega$ is the angular frequency, $\mathbb{m}$ contains the squared slownesses, $\mathbb{\Delta}$ is the discretized Laplace operator, $\mathbb{C}$ introduces boundary conditions (e.g., sponge-like absorbing boundary conditions such as perfectly-matched layers (Bérenger, 1994)) and $\mathbb{B}$ spreads the “mass” term $\omega^{2}\mathbb{C}\text{diag}(\mathbb{m})$ over all the coefficients of the stencil to improve its accuracy following an anti-lumped mass strategy (Chen et al., 2013).

A main drawback of the reduced formulation in equation (3.3) is that the objective function depends on the model parameters $\mathbb{m}$ through the oscillating inverse operator $\mathbb{A}^{-1}(\mathbb{m})$. This makes the inverse problem highly nonlinear, and hence prone to convergence to inaccurate minimizer when the initial $\mathbb{m}$ is not accurate enough.

An alternative objective function as proposed by van Leeuwen & Herrmann (2016) is defined as

[TABLE]

in which the wavefield is given by

[TABLE]

The wavefield $\mathbb{u}(\mathbb{m})$ is the solution of an augmented wave equation system, which gathers the wave equation $\mathbb{A(m})\mathbb{u}(\mathbb{m})=\mathbb{b}$ weighted by the scalar penalty parameter $\lambda>0$ and the observation equation $\mathbb{Pu}(\mathbb{m})=\mathbb{d}$ (van Leeuwen & Herrmann, 2013, their equation 6).

A main advantage of equation (3.6) over equation (3.3) is that its objective function depends on the model parameters $\mathbb{m}$ through a regularized inverse of the wave-equation operator $\left(\lambda\mathbb{A}(\mathbb{m})^{T}\mathbb{A}(\mathbb{m})+\mathbb{P}^{T}\mathbb{P}\right)^{-1}$, where the regularizer injects the prior information on the true wavefields, namely their sparse measurements $\mathbb{d}$, through the sampling operator $\mathbb{P}$. This wavefield reconstruction driven by the observations introduces a relaxation of the wave equation, namely source residuals, which are minimized to update the model parameters, equation 3.6, such that the wavefields, equation 3.7, are pushed back toward the wave equation.

The minimization problem, equation (3.6), can be solved with a variable projection approach by enforcing the closed-form expression of $\mathbb{u}$, equation 3.7, in the objective function, equation 3.6 (van Leeuwen & Herrmann, 2016), or with an alternating-direction strategy for $\mathbb{u}$ and $\mathbb{m}$ through a Gauss-Seidel like iteration (van Leeuwen & Herrmann, 2013). We will follow the latter option, where the model estimation problem, equation (3.6), reduces to a quadratic optimization problem

[TABLE]

which, in virtue of the wave-equation bilinearity (van Leeuwen & Herrmann, 2013; Aghamiry et al., 2019d), can be recast as

[TABLE]

with

[TABLE]

and $\mathbb{u}^{k}\equiv\mathbb{u}(\mathbb{m}^{k})$.

During wavefield reconstruction, equation 3.7, a small value of $\lambda$ is used during the early iterations to foster the data fit, hence mitigating cycle skipping, at the expense of the fidelity with which the wave equation is satisfied. Then, $\lambda$ is progressively increased in iterations such that the wave equation is satisfied at the convergence point (Fu & Symes, 2017). To avoid the tedious dynamic control of penalty parameters and converge in an automatic way toward accurate solution with fixed $\lambda$, Aghamiry et al. (2019d) introduced iteratively-refined WRI (IR-WRI). IR-WRI relies on an iterative refinement procedure, which finds its root in augmented Lagrangian method (Nocedal & Wright, 2006, Chapter 17). In the end, it simply consists of updating at the end of each iteration the right-hand sides of the wave equation and the observation equation, namely the data $\mathbb{d}$ and the source $\mathbb{b}$, with the running sum of the data and source residuals in iteration

[TABLE]

In the framework of linear inverse problems, $\mathbb{d}^{k+1}$ and $\mathbb{b}^{k+1}$ record the wavefield and model refinements performed at previous iterations to refine the wavefields and the model parameters at the current iteration from the residual errors only. This iterative refinement or defect correction (Böhmer & Stetter, 1984) procedure, which is lacking in WRI, is the key ingredient which controls the accuracy of the wavefield and model solutions at the convergence point when fixed $\lambda$ is used. The reader is also referred to Gholami et al. (2018); Gholami & Naeini (2019) for other recent applications of iterative refinement in the field of geophysics.

IR-WRI has shown promising results even for complicated velocity models when equipped with frequency continuation strategies and appropriate regularization (Aghamiry et al., 2019c, b). However, when the initial model is far from the true model, the phase of the reconstructed wavefields becomes inaccurate in areas that are located far away from the receivers. These phase inaccuracies can lead to an inaccurate model reconstruction of these regions, hence trapping the inversion into a local minimum. In order to mitigate this issue, we show in the next section how to solve the quadratic optimization problem, equation 3.9, with phase retrieval.

3.3 Wavefield inversion with phase retrieval (WIPR)

To estimate the model parameters without involving inaccurate phase information, we propose to replace the IR-WRI parameter-estimation problem given by equation 3.9 by the following phase retrieval problem

[TABLE]

In order to solve (3.12), we replace it by a surrogate majorizing function (as described in section 3.1 and Appendix References)

[TABLE]

This objective function is now quadratic and admits a closed form solution.

In order to better understand why phase retrieval helps to improve the solution, we give a tentative interpretation hereafter. From equation (3.7) it is seen that, at each iteration, the wavefield approximately satisfies the wave equation

[TABLE]

or

[TABLE]

where the approximation level is controlled by the parameter $\lambda$ and the approximation becomes an equality as $\lambda\rightarrow\infty$. The next iterate $\mathbb{m}^{k+1}$ is then found such that it minimizes the wave-equation errors generated by the wave-equation relaxation, equation (3.8). In this way, simple minimization via equation (3.8) forces the model to match both the amplitude and phase information of $\mathbb{y}(\mathbb{u}^{k})$. In order to reduce the imprint of the phase errors, a phase correction or phase alignment step can be applied to equation (3.15) before updating the model parameters:

[TABLE]

The optimal minimizer of equation (3.16) is given by $\mathbb{\phi}=\angle\mathbb{L}(\mathbb{u}^{k})\mathbb{m}^{k}$. Using the optimal $\mathbb{\phi}$, the next iterate is found by minimizing the corrected quadratic problem, as presented in equation (3.13). In this way, the solution of the inverse problem is less affected by the phase errors in $\mathbb{u}$ and thus hopefully less prone to convergence to a local minimum. The reader is referred to Jiang et al. (2017); Qian et al. (2017) for other applications of phase retrieval with alternating optimization.

Furthermore, we apply the iterative refinement procedure to WIPR in the same way as for IR-WRI by updating $\mathbb{b}$ and $\mathbb{d}$ according to equations 3.11a and 3.11b as outlined in the Algorithm 1. This implies that, although we recast the parameter-estimation subproblem as a phase retrieval problem (Algorithm 1, Line 5), we use both the phase and amplitude of the source and data residuals to update the right-hand sides of the wavefield-reconstruction and parameter-estimation subproblems (Algorithm 1, Lines 8 and 9).

3.4 Tikhonov-total variation (TT) regularized WIPR with bounding constraints

In this section, we implement TT regularization and bounding constraints in WIPR to stabilize the updates of the model parameters. By adding a convex regularization term $\|\mathbb{m}\|_{\text{TT}}$ to the cost function of WIPR, we get the following optimization for the $\mathbb{m}$-subproblem

[TABLE]

where $\mathcal{C}=\{\mathbb{x}\in\mathbb{R}^{n\times 1}~{}|~{}\mathbb{m}_{{l}}\leq\mathbb{x}\leq\mathbb{m}_{{u}}\}$ is the set of all feasible models bounded by the lower bound $\mathbb{m}_{{l}}$ and the upper bound $\mathbb{m}_{{u}}$ and $\|\mathbb{m}\|_{\text{TT}}$ is the TT regularization functional as defined by Aghamiry et al. (2019c, b)

[TABLE]

with

[TABLE]

in which the sum runs over all element, $\nabla_{\!\!i}$ are first order difference operators in direction $i$ and $\nabla_{\!\!ij}$ are second-order differential operators in directions $i$ and $j$.

The TT regularization is an infimal convolution-based combination of the second-order Tikhonov and TV regularizations and is suitable for recovering piecewise-smooth models (Gholami & Hosseini, 2013; Aghamiry et al., 2019b). It decomposes the model $\mathbb{m}$ explicitly into two components, $\mathbb{m}_{1}$ and $\mathbb{m}_{2}$, which have different statistical properties. Here, the blocky component $\mathbb{m}_{1}$ and the smooth component $\mathbb{m}_{2}$ are captured by the TV and the Tikhonov regularizations, respectively.

Subproblem (3.17) is solved with ADMM and the split-Bregman scheme, which de-couples the $\ell_{1}$ and $\ell_{2}$ components and bound constraints of the function through the introduction of auxiliary variables and solves each related subproblem in sequence (Goldstein & Osher, 2009). By doing this, we come up with a least-squares problem with a closed-form expression to update $\mathbb{m}_{1}$ and $\mathbb{m}_{2}$, a proximity subproblem to update the auxiliary variables and a gradient ascent step to update dual variables. The reader is referred to Aghamiry et al. (2019b) and Aghamiry et al. (2019c, Section 2.2.2) for the detailed derivation of the algorithm.

4 NUMERICAL EXAMPLES

We assess our workflow against the Marmousi II model (Martin et al., 2006) and the large-contrast 2004 BP salt model (Billette & Brandsberg-Dahl, 2004) when the inversion is started from a homogeneous velocity model and a 3 Hz frequency such that a considerable phase error is generated. In both cases, we consider a long-offset fixed-spread acquisition and perform the frequency-domain inversion in the 3 Hz - 13 Hz frequency band with a frequency interval of 0.5 Hz. Small batches of two frequencies with one frequency overlapping between two consecutive batches are successively inverted following a classical frequency continuation strategy. We apply the phase retrieval algorithm only during the first frequency batch processing to stabilize the inversion. Then, we proceed with regular IR-WRI at higher frequencies. The IR-WRI results obtained from an early WIPR stage will be referred to as IR-WRIpr to prevent confusions with results obtained when IR-WRI is applied from the first frequency-batch inversion. During IR-WRI/IR-WRIpr, we perform three paths through the frequency batches to improve the inversion results, using the final model of one path as the initial model of the next one (these cycles can be viewed as outer iterations of the algorithms). The starting and finishing frequencies of the paths are [3.5, 6], [4, 8.5], [6, 13] Hz, respectively, where the first element of each pair shows the starting frequency and the second one is the finishing frequency. We will compare the results of IR-WRIpr and IR-WRI to assess the improvements resulting from the initial WIPR step. We will also compare the results of IR-WRI and IR-WRIpr when regularization is applied or not such that the relative role of phase retrieval and regularization can be discriminated. Namely, we will illustrate that the role of phase retrieval is to extend the linear regime of the waveform inversion further when an inaccurate starting model is used, while regularization injects suitable prior to manage non-uniqueness of the solution in the poorly illuminated area. When regularization is applied, the same tuning is used during WIPR and IR-WRI, and during IR-WRI and IR-WRIpr. For all the numerical tests, we use a 9-point finite-difference staggered-grid stencil with PML boundary condition along the four edges of the model (no free-surface boundary condition is used) and anti-lumped mass to solve the Helmholtz equation, where the stencil coefficients are optimized to the frequency (Chen et al., 2013).

4.1 Marmousi II model

The Marmousi II model covers a $11500~{}m\times 4200~{}m$ spatial domain (Fig. 1a). The fixed-spread surface acquisition consists of 56 sources spaced 200 m apart, and 230 receivers spaced 50 m apart at the surface. The source signature is a 10 Hz Ricker wavelet, and the wavespeed in the homogeneous starting model is 3 km/s.

We first compare the results of WIPR and IR-WRI after 45 iterations of the [3,3.5] Hz frequency batch inversion when bound constraints are applied without any additional regularization (Fig. 2a-b). Bound constraints are used for all of the tests of this section from the first iteration. WIPR manages to reconstruct an accurate long-wavelength approximation of the Marmousi model (Fig. 2b), while IR-WRI clear fails to capture the kinematic trend of the model (Fig. 2a). This statement is clearly illustrated by the direct comparison between the true model, the initial model, and the reconstructed IR-WRI and WIPR models along with three vertical logs at horizontal distances of 4, 8.2, and 10 km (Fig. 3a). We also apply TT regularization during WIPR and IR-WRI (Figures 2c-d and 3b). The TT regularization improves the WIPR and IR-WRI results in the poorly-illuminated area, although it does not remove high-velocity artefacts near the left-bottom end of the IR-WRI velocity model (Fig. 3b, x=4 km).

Since the model parameters are updated from the reconstructed wavefields, the accuracy of these wavefields is a visible indicator of the reliability of the estimated model parameters. We compare the real part of the 3 Hz wavefields reconstructed by IR-WRI and WIPR without and with TT regularization (Fig. 4) with the true wavefield (Fig. 1b) for a source located at the horizontal distance of 6 km. Both IR-WRI and WIPR wavefields are in phase with the true wavefield along the receiver line since we tune the penalty parameter such that the data are matched from the first iteration. However, the WIPR wavefield matches much better the true wavefield than the IR-WRI counterpart in depth. When regularization is used, the accuracy of the wavefield reconstruction is mostly improved in the poorly-illuminated zones, mainly near the bottom ends of the model.

We continue with IR-WRI at higher frequencies using the final IR-WRI and WIPR models of the [3, 3.5] Hz inversion (Fig. 2) as initial models. The stopping criterion of iterations is a maximum of 30 iterations per batch or

[TABLE]

where $F$ refers to the Frobenius norm. The final IR-WRI and IR-WRIpr velocity models are shown in Fig. 5, while a direct comparison between them and the true model are plotted in Fig. 6 at distances 4, 8.2 and 10 km. IR-WRI performed a total of 382 and 369 iterations without and with TT regularization, respectively (Fig. 5a,c), while IR-WRIpr performed a total of 201 and 234 iterations (Fig. 5b,d). The smaller number of iteration required to satisfy the stopping criterion of iterations performed by IR-WRIpr compared to IR-WRI highlights the role of the initial model built by WIPR during the [3, 3.5] Hz inversion (Fig. 2b,d). Without regularization, IR-WRIpr outperforms IR-WRI in the deep part, namely near the reservoir at around 8 km distance and near the salt layers near the ends of the models. Regularization improves significantly both IR-WRI and IR-WRIpr results in these areas. However, IR-WRIpr still provides superior results in the deep part.

4.2 The large contrast 2004 BP salt model

We consider now a target of the large contrast 2004 BP salt model located on the left side of the model (Fig. 7a). The 2004 BP salt model mainly consists of a simple sediment background with a complex rugose multi-valued salt body, sub-salt slow velocity anomalies related to over-pressure zones, and a fast velocity anomaly to the right of the salt body. The selected subsurface model is 16250 m wide and 5825 m deep and is discretized with a 25 m grid interval. The surface fixed-spread acquisition consists of 66 sources spaced 250 m apart, and 131 receivers spaced 125 m apart. We perform forward modelling using a 10 Hz Ricker wavelet as a source signature. The starting velocity model for inversion is homogeneous with a velocity of 3 km/s.

4.2.1 IR-WRI versus WIPR at low frequencies

We first perform 45 iterations of bound-constrained WIPR and IR-WRI for a starting [3,3.5] Hz frequency batch. Note that we turn on the bound constraints after 20 iterations to discriminate their role from those of the TT regularization and optimization scheme (see later discussion in the text). The velocity models reconstructed by bound-constrained IR-WRI and WIPR with and without TT regularization are shown in Fig. 8, while the direct comparison between the true model, the initial model, and the models reconstructed by IR-WRI and WIPR are shown in Fig. 9 along three vertical logs at horizontal distances of 3.5 km, 7.7 km, and 12 km. The most striking feature is that IR-WRI overestimates velocities above the salt and mispositions the top of the salt layers accordingly, while WIPR captures the top salt much more accurately (Fig. 9a, green versus red curves, x=3.5 km and 7.7 km). This overestimation of the shallow velocities by IR-WRI is also shown in the right part of the model away from the salt (Fig. 9a, red curve, x=12 km). On the other hand, WIPR fails to reconstruct the bottom part of the thick salt body and overprints unrealistic variations on the subsalt area when TT regularization is not used (Fig. 8b and Fig. 9a, green curve, x=3.5 km). These artefacts are likely generated by the increased ill-posedness of WIPR with depth, which itself results from the lack of phase information. This ill-posedness is however, efficiently mitigated by the TT regularization (Fig. 8d). The TV component of the TT regularization helps to sharpen the top and bottom of the salt layers (Fig. 9b, green curve, x=7.7km), recovers better the geometry of the salt layer (Fig. 9b, green curve, x=3.5km), while the Tikhonov component stabilizes the inversion in-depth and drives it toward smooth reconstruction of the poorly-illuminated subsalt structures (Fig. 9b, green curve, x=12km). During IR-WRI, TT regularization fails to achieve these goals (Fig. 8c) because the mispositioning in depth of the velocity structures generated by the increased nonlinearity of IR-WRI are too heavy to be corrected by the regularization (Fig. 9b, red curve).

We compare the true 3-Hz wavefield shown in Fig. 7b with the wavefields reconstructed by IR-WRI and WIPR without and with TT regularization (Fig. 10). Like the Marmousi test, both IR-WRI and WIPR wavefields are in phase with the true wavefield along the receiver line since we tune the penalty parameter such that the data are matched from the first iteration. However, the WIPR wavefields match much better the true wavefield than the IR-WRI counterpart in the first two kilometers of the model. This is pointed by the arrows in Fig. 10, which show that the IR-WRI wavefield transmitted into the salt body is out-of-phase with respect to the true wavefield, unlike the WIPR wavefield.

We also computed time-domain seismograms in the TT-regularized IR-WRI and WIPR models (Fig. 8c,d) with a 10-Hz Ricker wavelet for a source located at 15.8 km (Fig. 11). We show quite significant travel-time mismatches at long offsets between the true seismograms and those computed in the IR-WRI model (Fig. 11a), which result from the overestimated velocities above the salt shown in Fig. 8c, while the match between the true seismograms and those computed in the WIPR model is far better. This poor data fit achieved by IR-WRI highlights both the inaccuracies of the reconstructed velocity model and the limited accuracy with which the wave equation has been satisfied at the convergence point.

This is further supported by Fig. 12, which shows the relative model error (Fig. 12a) and the joint evolution of the data misfit ($\|\mathbb{Pu}^{k}-\mathbb{d}\|_{2}$) and wave-equation error ($\|\mathbb{A(m}^{k})\mathbb{u}^{k}-\mathbb{b}\|_{2}$) in iterations (Fig. 12b). Also, the evolution of the data misfit and wave-equation error in iterations are shown separately in Figs. 12c and 12d, respectively. First, these error curves show that WIPR converges toward a more accurate velocity model than IR-WRI (Fig. 12a). Accordingly, it shows that WIPR satisfies more accurately the wave-equation constraint than IR-WRI (Fig. 12d), although the parameter estimation has been recast as a phase retrieval problem rather than as a wave-equation error minimization. Fig. 12a also shows that bound constraints have the least impact on the TT regularized WIPR (blue curve) since the descent direction is rather continuous around iteration 20 (when bound constraints are activated) . This highlights the more stable descent direction followed by TT regularized WIPR relative to IR-WRI.

4.2.2 IR-WRI versus $\text{IR-WRI}_{\textit{pr}}$

We continue with IR-WRI at higher frequencies using the final IR-WRI and WIPR models of the [3, 3.5] Hz inversion (Fig. 8) as initial models. The stopping criterion of iterations is 10 iterations per batch. The final velocity models built by IR-WRI and IR-WRIpr are shown in Fig. 13, while a direct comparison between them and the true model are plotted in Fig. 14 at distances 3.5 km, 7.7 km and 12 km.

Let’s first focus on the results obtained without TT regularization (Fig. 13a,b). IR-WRI apparently reconstructs the salt body with good accuracy (Fig. 13a), although the shallow kinematic inaccuracies of the starting IR-WRI model shown in Figs. 8a and 9, red curves. This highlights the potential of the search-space expansion implemented in IR-WRI to manage cycle skipping. The final IR-WRIpr velocity model shows low-velocity artefacts in the salt body at around 3 km distance (Fig. 13b), which are probably inherited from the underestimated velocities below the top of the salt shown in Figs. 8b and 9a, green curve, x=3.5 km. For both inversions, the reconstruction of the subsalt structures are quite noisy. To gain more quantitative insights on the relative accuracy of the IR-WRI and IR-WRIpr models, we show the difference between these velocity models and the true model in Fig. 15(a-b) and outline the model error (ME) for the velocity models as well as the least-squares data and source misfit at the convergence point in Table 1. Figure 15(a-b) clearly shows that the IR-WRIpr model is as a whole more accurate than the IR-WRI counterpart, although the above-mentioned local artefacts.

The TT regularization removes the local artefacts shown in the salt body of the IR-WRIpr model and improves the subsalt imaging (Fig. 13d) more significantly than for the IR-WRI model (Fig. 13c). The difference between the true model and the IR-WRI/IR-WRIpr models (Fig. 15(c-d)) and the model error (ME) provided in Table 1 outline the superior accuracy of the IR-WRIpr model relative to the IR-WRI counterpart. These relative model accuracies are consistent with the data and source misfit achieved by IR-WRIpr and IR-WRI at the convergence point (Table 1).

5 DISCUSSION

With the large contrast BP salt model, we have shown that the phase inaccuracies of the reconstructed wavefields prevent accurate estimation of the shallow sediments and accurate positioning in depth of the top of the salt after the first frequency batch inversion. In contrast, WIPR reconstructs more accurately the velocity gradient in the sedimentary cover as well as the sharp velocity contrast on top of the salt. This accurate positioning of the top salt is indeed important to speed up the convergence and improve the solution at higher frequencies.

To further illustrate the sensitivity of IR-WRI to wavefield amplitudes and phases in large contrast media, we show in Figs. 16 and 17 the amplitude and phase of $\Delta\mathbb{u}$ for the 3 Hz frequency and a source located at 7 km of distance. In these figures, $\Delta\mathbb{u}$ have been computed in the true velocity model, and in the TT regularized IR-WRI and WIPR velocity models obtained after the first frequency batch inversion (Figs. 8c,d). We remind that $\Delta\mathbb{u}$ is one of the terms forming the right-hand side $\tilde{\mathbb{y}}$ of the quadratic parameter-estimation subproblem together with the updated source $\mathbb{b}+{\mathbb{b}^{k}}$, Algorithm 1, Line 5. As such, $\Delta\mathbb{u}$ can be viewed as the observable of the parameter-estimation subproblem.

The amplitudes of $\Delta\mathbb{u}$ computed in the true velocity model show the dominant imprint of the reflection and the refraction from the top of the salt (Fig. 16a, shallow black zone). Also, the geometry of the seismically-transparent salt body is almost perfectly delineated. This shows that the amplitudes of $\Delta\mathbb{u}$ intrinsically embeds high-resolution information on the subsurface.

This high-resolution potential can be simply illustrated with the wave-equation bilinearity as discussed below. Let us remind the wave equation for a single source and a single frequency

[TABLE]

which leads to the following linear equation of the parameter-estimation subproblem:

[TABLE]

Compared to equation 3.5, we drop the matrix $\mathbb{B}$ and $\mathbb{C}$ which are not important for the discussion. The equality implies that $|\mathbb{L(u})\mathbb{m}|=|\mathbb{y}(\mathbb{u})|$ and $\angle\mathbb{L(u})\mathbb{m}=\angle\mathbb{y}(\mathbb{u})$. If $\mathbb{u}$ is the true wavefield (if it would be recorded everywhere) then we would be able to recover the true model exactly by minimizing

[TABLE]

It is seen that, basically, the true model is the solution of a weighted phase retrieval problem (as indicated by the middle term) where the complex exponential of the phase serves as the weight.

Indeed, for an approximate wavefield which is obtained from eq. (3.7) with a rough velocity model or a large penalty parameter the equality in eq. (5) is not necessarily satisfied because, in this case, $\angle\mathbb{L(u})\mathbb{m}\neq\angle\mathbb{y}(\mathbb{u})$. The proposed WIPR algorithm simply forces this equality to be satisfied at each iteration by aligning the phase of $\mathbb{y}$ with that of $\mathbb{L(u})\mathbb{m}$ (right term in equation 5).

The ability to reconstruct the true model from amplitude information only can be illustrated more explicitly if we assume that $\mathbb{m}$ is real in equation 5.1 (non attenuating medium). The matrix $\mathbb{L(u})$ is diagonal and contains the so-called virtual source (equation 5.2).

[TABLE]

where $\circ$ denotes element wise multiplication. Indeed, we can equally reconstruct the true $\mathbb{m}$ from the magnitude of $\mathbb{y}$ as

[TABLE]

where $1/\mathbb{u}$ and $1/|\mathbb{u}|$ are element-wise reciprocal of $\mathbb{u}$ and $|\mathbb{u}|$, respectively. The Marmousi II model computed with equations 5.4 and 5.5 by using the true wavefield (Fig. 1b) is shown in Fig. 1c and matches exactly the true one (Fig. 1a).

The Fig. 16c, as well as the vertical profile extracted at 6.5 km of distance, confirm that WIPR manages to match the amplitudes of the true $\Delta\mathbb{u}$ down to the top of the salt at 1.25 km depth, while IR-WRI (Fig. 16b) allows for an amplitude fit down to a maximum depth of 0.5 km only, consistently with the accuracy of the velocity fields reconstructed by IR-WRI and WIPR (Fig. 8c,d).

It is also instructive to look at the phase of $\Delta\mathbb{u}$ (Fig. 17). To facilitate the comparison of the phases computed in the true model (Fig. 17a) and in the reconstructed models (Fig. 17b,c), we superimpose in Fig. 17b,c the contours of the phase computed in the true model at the spatial positions where the phase wraparounds, namely when it jumps from $\pm\pi$ to $\mp\pi$ (Fig. 17b,c). We clearly show that WIPR reconstructs the phase of $\Delta\mathbb{u}$ much more accurately than IR-WRI around the salt body. This is further illustrated by the difference between the phases computed in the true model and the reconstructed models along the vertical profile at 6.5 km distance (Fig. 17d). The phase difference associated with IR-WRI (solid red curve) reaches a value of $-\pi$ at the top of the salt (1.25 km depth) showing that at these depths the two wavefields are cycle skipped, while the phase differences associated with WIPR (solid green curve) remains far below this limit down to around 2.7 km depth.

Conversely, it is worth noting that IR-WRI matches better the phase of the true wavefield than WIPR in the smooth bottom-right part of the model (Compare Figs. 17b and 17c). This highlights that the sensitivity of the inversion to amplitudes decreases rapidly with depth due to lack of illumination, geometrical spreading effects and energy partitioning at interfaces. We also want to stress that, although we update the wavespeeds from the amplitudes of $\mathbb{y}$, we update $\mathbb{y}$ at each iteration with the phase and amplitude of the wave equation error (the source residuals) in the framework of the ADMM optimization (Algorithm 1). This right-hand side updating re-injects at each iteration the phase error as a defect correction term in the optimization.

For sake of completeness, we show also the amplitude and phase of $\Delta\mathbb{u}$ for the Marmousi case study in Figs. 18 and 19. Fig. 18 shows how WIPR reproduces more accurately the wide-angle reflection from the dipping layer on the left side of the source than IR-WRI. Fig. 19 shows that WIPR better reconstructs the true phase of $\Delta\mathbb{u}$ than IR-WRI.

Although the application of the phase retrieval algorithm has been limited to the scalar Helmholtz equation in this study, the method has been formulated such that more complex physics and multi-parameter reconstruction can be viewed (the reader is referred to Aghamiry et al. (2019e) and Aghamiry et al. (2019a) for application of IR-WRI in visco-acoustic VTI media). The key feature allowing for these extensions is the linearization of the non-convex phase retrieval problem performed by the MM approach, as reviewed in Appendix References.

6 CONCLUSIONS

We present a preliminary application of bound-constrained and TT regularized WIPR. Using the large contrast BP salt model, we show that, when the inversion is started from scratch (homogeneous initial velocity model), WIPR improves significantly the reconstruction of the sedimentary background and the top of the salt during the early iterations of the inversion (that is, when the wavefields are not yet accurately reconstructed far away from the receivers). Sparsity-promoting regularization is a necessary ingredient to stabilize the phase retrieval inversion at greater depths where the amplitude information becomes more challenging to extract. Sparsity-promoting regularization combined with phase retrieval maybe also two complementary tools to relax the need for finely-sampled stationary-recording acquisitions. Phase retrieval can also be implemented in classical FWI. Also, WIPR can be theoretically extended to multi-parameter reconstruction in attenuating media. Further theoretical and numerical works are, however, necessary to flesh out the potential and limits of phase retrieval in the FWI technology.

7 ACKNOWLEDGEMENTS

We would like to thank editors Huajian Yao, Louise Alexander, and reviewers Hejun Zhu and an anonymous for their comments which help improving the manuscript. This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.fr), sponsored by AKERBP, CGG, CHEVRON, EQUINOR, EXXON-MOBIL, JGI, PETROBRAS, SCHLUMBERGER, SHELL, SINOPEC, and TOTAL. This study was granted access to the HPC resources of SIGAMM infrastructure (https://www.oca.eu/fr/ mesocentre-sigamm) and CINES/IDRIS/TGCC under the allocation A0050410596 made by GENCI.