A priori estimates of attraction basins for velocity model reconstruction by time-harmonic full-waveform inversion and data-space reflectivity formulation

The full-waveform inversion (FWI) approach to seismic inversion consists of minimizing the misfit between the recorded and synthetic data with respect to the earth parameters. It uses the full waveform of the recorded waves as data, which makes it free of horizon picking, and hence a good candidate for automatic inversion. It was first introduced by Bamberger et al. (1977, 1979) for the 1D wave equation, and it was followed, in two dimensions, by the work of Lailly (1983) and Tarantola (1984). In particular, Lailly (1983) links the gradient of the misfit functional, computed by the adjoint-state method, with the migration imaging operator used in geophysics. Initially carried out in the time domain, FWI can also be used in the frequency domain, as promoted by Pratt and Worthington (1990), Pratt and Goulty (1991), and Pratt et al. (1996). The use of complex frequencies was later introduced by Shin and Cha (2008, 2009), giving rise to the so-called Laplace and Laplace-Fourier domain FWI. Early numerical experimentations by Tarantola (1987) and Gauthier et al. (1986) already show the strength and the weakness of FWI: It can produce very clear and detailed images of the subsubsurface (strength), under the imperative condition that a good background velocity model is available, i.e., a priori information (weakness). In fact, the misfit functional with respect to background parameters shows local minima due to phase shift/cycle skipping, which hamper the reconstruction by local algorithms, unless the minimization starts inside the attraction basin.

In recent years, with the increase in computational power, the solution of the FWI approach to seismic imaging by local optimization algorithms has been extensively developed. The iterative minimization is usually conducted with Newton-type algorithms (Pratt et al., 1998). To avoid the formation of the dense Hessian (and its inversion), multiscale Newton methods are studied by Akcelik et al. (2002). For similar reasons, the pseudo-Hessian is studied by Shin et al. (2001) and the iterative conjugate gradient approach is implemented by Métivier et al. (2013). Estimating the step length is another challenge for the iterative minimization problem, and it is often conducted using the line search method (Nocedal and Wright, 2006). Alternatives for a more precise step exist, such as the trust region method (Eisenstat and Walker, 1994) or the maximum projected curvature method of Chavent (2010), but with increased computational complexity. The performance of FWI is further compared in the context of inverse scattering by Barucq et al. (2018), depending on the choice of search direction and the line search method.

However, none of these local optimization methods can avoid local minima. The determination of the full velocity model, including the low spatial frequency background, requires very low-frequency data (as unrealistic as 1 Hz) in marine seismic applications, coupled with a resolution from low- to high-frequency content to avoid “cycle skipping” (Bunks et al., 1995; Sirgue and Pratt, 2004). To overcome the presence of local minima when solving the FWI problem by local optimization algorithms, alternative parameterizations have been proposed. They are based on a decomposition of the earth parameter into background and reflectivity parameters, resulting in an increase of the computational complexity. The differential semblance optimization by Symes and Carazzone (1991) extends the depth reflectivity model to account for the various illuminations in the data, defining a semblance criterion whose minimization produces the desired background model; the migration-based traveltime (MBTT) approach of Chavent and Clément (1992), Clément et al. (2001), and Chavent et al. (2015) parameterizes the depth reflectivity by a data-space reflectivity, which is used to store the reflectivity during the update of the background. Note that the so-called reflection FWI (Xu et al., 2013) retains from MBTT the migration-demigration process, but it differs in that it does not use the data-space reflectivity concept and it uses the gradient of the classical FWI data misfit as the search direction.

In this paper, we investigate the optimizability of FWI in the frequency domain, using the geometric tools set up by Barucq et al. (2019b). We analyze the size of the attraction basin, in which the misfit functional is guaranteed to have no local minimum, and the associated maximal distance allowed between the data and the attainable set (i.e., the allowed noise level).

We first present the FWI minimization problem for a time-harmonic acoustic model, and its MBTT reformulation. The second section analyzes the influence of (complex) frequency and search geometry on the size of attraction basins and tolerable error levels for the classical FWI formulation. We also demonstrate the accuracy with experiments of FWI reconstruction. Then, we compare the optimizability of classic FWI and the MBTT reformulation. We quantify the improvement obtained with the MBTT approach which, thanks to the migration-demigration that eliminates phase shift/cycle-skipping problems, leads to larger attraction basins. We further indicate the guidelines for an efficient implementation of MBTT.

The main difficulty with FWI is the determination of the background velocity part of $m$ because the misfit functional suffers from many local minima caused by the phase shifts between the observed and simulated signals. These can only be avoided if the data contain extremely low frequencies, as illustrated in our experiments below. As an alternative, the following MBTT reformulation of the problem is designed to alleviate these difficulties.

We shall refer to $s$ as the data-space reflectivity, in opposition to the depth reflectivity $r$⁠. The data-space reflectivity $s$ is one of the unknowns that we search for, together with $p$⁠, in the MBTT formulation. At the optimum, $s$ can be understood as a demultiplied and deconvolved version of the data themselves, which hence provides a first approximation of $s$⁠.

As indicated above, the motivation for parameterizing the reflectivity by an element $s$ of the data space is to eliminate the phase shifts induced in the simulations by changes in the background $p$⁠: The events in the synthetics are obtained from the data-space reflectivity $s$ by migration, followed by simulation with the same kinematics. Hence, they are expected to have the same phase as those of $s$⁠, as illustrated in our experiments. In addition to controlling the phase, the stack involved in the migration turns, for a fixed data-space reflectivity, the data misfit into a coherency measure for the current background (Chavent, 2017). On the other hand, the MBTT model decomposition induces greater computational cost because of the dependency on $DF(p)$ in equation 9. Evaluations are derived from adjoint-state formula, and we refer to Barucq et al. (2019b) in their Appendix C for more details.

The iterative reconstruction has been shown successful by Chavent et al. (2015), starting with 5 Hz data, and it follows alternative updates of the two model components:

1. 1)

   Start the reconstruction with the recovery of the reflectivity part $s$⁠, using a fixed $p$⁠.

2. 2)

   Fix the recovered reflectivity $s$⁠, and iterate to reconstruct the background part of the model $p$⁠.

3. 3)

   Continue with alternate updates, in a succession of increasing frequency.

The optimizability study that we provide below quantifies the enlargement of the attraction basins provided by the MBTT decomposition, in particular for the second step.

We review the properties that allow us to represent the optimizability of least-squares problems such as in equations 5 and 12. We refer to optimizability as the absence of local minimum in the misfit functional, which ensures the convergence of the minimization algorithm toward the appropriate (hence global) minimum, provided that the data (observations) are “not too far” from the attainable set $F(M)$⁠. The precise analysis can be found in Barucq et al. (2019b) (section 2). Here, we want to give a comprehensive overview of the tools that we use to analyze least-squares problems. It is important to note that the estimates to quantify optimizability are defined a priori; i.e., they only depend on the forward problem (not on the data) and can be computed without — and preferably before — performing any inversion.

In seismic applications, the model space $M$ does not usually contain much information because subsurface characterization is mostly unknown. In particular, the data phase shifts resulting from an error in the long spatial wavelength of the initial guess produce local minima for the problem stated in equation 5, hampering the search for the global minimum. Therefore, our idea is to analyze, from an initial model $m0$⁠, what is, for a given direction $u$⁠, the largest interval in which the misfit functional is guaranteed to contain no parasitic stationary points — and hence no parasitic local minima — beside the global minimum, provided that the error in the data is smaller than some tolerable level. We refer to such an interval $[m0−Δm0uu,m0+Δm0uu]$as a directional attraction basin of size $Δm0u$ centered around $m0$ in the normalized direction $u$⁠, i.e., with $‖u‖=1$⁠.

The numerical estimation of $RG(P)$ will be the basis for the a priori determination of attraction basins in the next sections.

All the above formulas have been illustrated with a plane curve $P$⁠. They remain true when $P$ is a curve of a Hilbert data space (at the exception of equation 25).

We apply the methodology to the a priori determination of directional attraction basins $[m0−Δm0uu,m0+Δm0uu]$⁠, which will allow us to analyze the behavior of seismic inverse problems and to compare formulations: The size $Δ$ of a directional attraction basin in a descent direction tells us how far we can move away in this direction without being stopped by parasitic local minima. Our objective is hence to determine (see the illustration in Figure ²)

1. 1)

   The size$Δm0u$ of the directional attraction basin centered at $m0$⁠. The larger $Δm0u$ is, the better the least-squares problem is amenable to minimization by local algorithm because we allow a larger area for investigation. In our numerical experiments, we shall scale the estimate with the norm of the nominal model $‖m0‖E$ to provide the relative (to the model) quantity.

2. 2)
   The associated tolerable error level $RG,m0u$⁠. It is the largest tolerable error on the data $d$⁠, which ensures the absence of parasitic local minima for the least-squares objective function

   $t∈[0,1]↝12‖F(m0+(2t−1)Δm0uu)−d‖D2,$

   (27)

   over $[0,1]$⁠. The larger $RG,m0u$ is, the better is the robustness of the minimization procedure to noise in the data. In our numerical experiments, we divide the estimates with the norm of the synthetic data $d0=F(m0)$ to provide relative (to the data) quantity.

We shall use two types of estimates:

1. 1)

   $Θ$ estimates of $Δm0u$⁠, based on equation 23, where optimizability is obtained by ensuring the sufficient condition $Θ(P)≤π/2$⁠. In this case, $RG(P)=R(P)$⁠, so the tolerable error level is given by the minimum over the $[0,1]$ interval of $R(t)$ given in equation 18.

2. 2)

   $RG$ estimates of $Δm0u$⁠, based on equation 17, where optimizability is obtained by satisfying the necessary and sufficient condition $RG>0$⁠. In this case, the associated tolerable error level $RG$ has to be computed by numerically evaluating the infimum of $RG(t,t′)$ in equation 15.

To compute the exact estimates, deflection maps and global radius maps are computed, which display the values of $Θ(t,t′)$ and $RG(t,t′)$ between the points $t,t′$ of $[−Δ,+Δ]×[−Δ,+Δ]$⁠. On the diagonal of the maps, where $t=t′$⁠, $RG$ is not defined in equations 30 and 31, and we indicate instead the values of $R(t)$ given in equations 18 and 19, which represent the limits of $RG(t,t′)$ when $t′→t$⁠. We review the different steps for the computation of the maps in Algorithm 1.

One can then read on these maps (see the numerical experiments below):

1. 1)

   The exact $Θ$ estimate of the attraction basin size $Δm0u$⁠, given by the largest square centered at $(0,0)$ where $Θ(t,t′)≤π/2$ for all $t$⁠, $t′$⁠.

2. 2)

   The exact $RG$ estimate of the attraction basin size $Δm0u$⁠, given by the largest square centered at $(0,0)$ where $RG(t,t′)>0$ for all $t$⁠, $t′$⁠.

During this process, when the size $Δ$ of the investigated square increases from zero to the exact $Θ$ estimate, the associated exact tolerable error $Rm0u=inf−Δ≤t≤ΔR(t)$ decreases from its value $R0$ at $m0$ to the tolerable error $Rm0u$ of the $Θ$-attraction basin. When the size $Δ$ of the square increases further to the exact $RG$ estimate, the tolerable error $RG=inf−Δ≤t,t′≤ΔRG(t,t′)$ continues to decrease, until it reaches the value zero, which gives the limit of the $RG$-attraction basins. The exact estimates are relatively demanding because they require evaluation of $RG$ and $Θ$ for many couples $(t,t′)$⁠, and we now introduce the cheaper local estimates.

Our objective in this section is to give a comprehensive understanding of the behavior of the seismic FWI problem equation 5 with respect to parameters such as the frequency, the geometry of the target, and the parameterization. We select a nominal model $m0$ and compute for different (normalized) directions $u$⁠:

1. 1)

   The size $Δm0u$ of an attraction basin, which represents the interval, in the model space, where the misfit functional is guaranteed to have only one minimum.

2. 2)

   The associated tolerable error level $Rm0u$⁠, which is the maximal distance allowed between the data and the attainable set (i.e., the maximal tolerable error level).

When $m0$ is smooth, the direction $u$ can be seen as the search direction or the geometry of the unknown parameter to be reconstructed in the framework of inversion.

We consider the 2D geophysical setup for the acoustic Helmholtz equation 1, and we take the squared slowness $m≔c−2$ as unknown parameter, expressed in $s2m−2$⁠. In our experiments, we consider a $9.2×3km$ domain, and we take $nsrc=19$ sources and $nrcv=183$ receivers associated with each source. Both are located near the surface: The sources are at 10 m depth, and the receivers, which remain in the same position for all sources, are at 40 m in depth.

We first investigate individual frequency, where the forward map $Fω$ at frequency $ω$ for parameter $m$ is a vector of $C(nrcv×nsrc)$⁠. The nominal model $m0$ is the smooth background velocity pictured in Figure ³, in which the velocity varies linearly with depth, according to the fact that no information is known on the subsurface structures to reconstruct.

We investigate the behavior of attraction basins for different directions $u$⁠, which are chosen to be based on straight reflectors, or extracted from the Marmousi model, a synthetic wavespeed model designed by the Institut Français du Pétrole in the late 1980s; we refer to Versteeg (1994) for additional details. They are shown in Figure ⁴: $u1$ has a single, slim reflector (Figure ⁴), $u2$ has two slim reflectors (Figure ⁴), $u3$ has a single large reflector (Figure ⁴), and $um$ is obtained by high-pass filtering from the Marmousi model of Figure ⁴. The amplitude for the directions is selected such that $uk$ is normalized: $‖uk‖=1$⁠, for all $k$⁠.

Given a perturbation size $Δ$⁠, one can think of $m1=m0+Δu$ as being the target velocity, which one tries to retrieve starting from the smooth model $m0$⁠, by minimization of the FWI misfit function in direction $u$ by a local algorithm (it is analogous to the update formula given in equation 6 in Newton-type algorithm). The question is: How large can $Δ$ be while still guaranteeing that the algorithm will not stop before the global minimum, provided that the error in the data is smaller than some $R$ or $RG$?

The determination of $Δm0u$ and $Rm0u$ with the local $Θ$ estimate, equations 33 and 34, and $Rm0u$ with equation 34, only requires the computation of the first- and second-order derivatives in the direction $u$⁠, obtained by differentiation of the Helmholtz equation (Barucq et al., 2019b, Appendix A). Figure ⁵ shows the evolution of $Δm0u$ (Figure ⁵) and $Rm0u$ (Figure ⁵) for frequencies between 0.5 and 15.0 Hz, for the four directions $u$ of Figure ⁴.

We observe that, in all cases, the size $Δm0u$ of the directional attraction basin decreases with increasing frequency, meaning that the largest basin of attraction is obtained for the lowest frequencies. This is the expected behavior that indicates that lower frequencies give a better chance for the iterative minimization to converge, especially when no prior information is known for the initial model.

Concerning the effect of the perturbation direction on attraction basins, Figure ⁵ shows that the size of the basins decreases in the order $u1→u2→um→u3$⁠. This can be understood as follows: Directions $u1$⁠, $u2$⁠, and $um$ are essentially reflective and contain no low-frequency components. Hence, increasing $Δ$ will not change the phase of backscattered data, but because of multiple reflections, some nonlinearity will occur, which limits the size of the attraction basins. The directions $u1$ and $u2$ with fewer multiple reflections have the largest basins of attraction. Then, comes $um$⁠, which has many reflectors and hence many multiples. On the contrary, direction $u3$ has a strong low-frequency component besides its reflectivity. Hence, increasing $Δ$ will generate a phase shift in the backscattered data, which is likely to create local minima more quickly, and one expects the size $Δm0u$ of the attraction basin in the direction $u3$ to be much smaller than in the three other directions.

Concerning the largest tolerable error, Figure ⁵ shows that $Rm0u$ increases with the frequency, meaning that the low frequencies are more likely to be affected in the presence of noise. We notice the exception of $u3$ for which the allowed error remains relatively similar (and worse than the others) with frequency.

These results clearly show that it is harder to converge toward the salt dome than it is to recover Marmousi-like structures. Because of its low spatial frequency component, the former reduces the size of the attraction basin and is less tolerant to noise in the data. This situation is further illustrated in the context of model reconstruction in the subsection below. Note that the order of magnitude in Figure ⁵ and ⁵ is approximately $10−1$ of the norms of the model and data, respectively.

As an alternative to confirm these preliminary observations, we compute $Δm0u$ using the exact $Θ$ and $RG$ estimate. It requires the computation of deflection and global radius maps on the investigated interval, from which the estimate is extracted. In Figures ⁶ and ⁷, we show the resulting maps in the direction of the Marmousi reflectors $um$ of Figure ⁴ and the single large reflector $u3$ of Figure ⁴, using frequencies of 4 and 7 Hz. In Table 1, we show the numerical values of the interval size centered around $m0$⁠, depending on the estimation.

The maps confirm the results of the previous experiment: The exact $Θ$ estimates of $Δm0u$ (Figure ⁶ and ⁶) are of the same order of magnitude, but they are slightly larger than the local $Θ$ estimates of Figure ⁵ (see Table 1). Similarly, the exact $RG$ estimates of $Δm0u$ (Figure ⁶ and ⁶) are larger (Table 1), but at the price of a smaller tolerable error on the data. We also note a much more consistent pattern in the deflection for $um$ compared to $u3$⁠. Comparing the formulas, we observe that the inexpensive local $Θ$ estimates already provide accurate values for classic FWI problems, close to the more refined (and computationally more intensive) exact $Θ$ and $RG$ estimates. Therefore, we shall only use the local $Θ$ estimate in the remainder of this section.

We now investigate the influence of the initial model $m0$ on the size of attraction basins. We consider two initial models: the acoustic Marmousi medium $mm$ and a model $mc$ encompassing objects of high contrast. These are pictured, together with their respective local $Θ$ estimates $Δmju$ given in equation 33, in Figure ⁸, for the two directions $u1$ (one reflector) and $u3$ (one large reflector). We see that $Δmju$ is barely affected by the choice of initial model for the computation. We still observe the decrease in the size of the interval with increasing frequency, and the deterioration in the direction of the large reflector $u3$⁠. Clearly, when comparing amplitudes of our estimates, the direction $u$ has much more influence than the initial model $m0$⁠.

From this first set of experiments, we clearly see that the frequency should evolve from low to high values to facilitate the convergence. When no prior information is known for the initial model, one should take advantage of the larger interval size $Δm0u$ given by low frequencies. Here, we motivate the well-known frequency progression of, e.g., Bunks et al. (1995) and Sirgue and Pratt (2004), from quantitative estimates of the size of the basin of attraction, which allows the quantification of the benefits. However, the low frequencies require accurate data, i.e., close to the $F(M)$⁠, as illustrated with the estimates of $Rm0u$⁠. In addition, we mention the fact that seismic data are particularly affected by noise at low frequencies, which are sometimes unusable. This highlights where resides the difficulty: The low frequencies that are required for convergence may be unusable due to the noise, or the inaccuracy may violate the curvature condition.

We have seen that sequential progression in frequency must be conducted from the low to the high regime. We now investigate the case of frequency bandwidth: when the forward problem encompasses not only one, but several frequencies. This approach is particularly natural for the time domain inverse problem in which several frequencies are simultaneously contained in the data. The forward problem is now a vector of $C(nrcv×nsrc×nω)$⁠, where $nω$ is the number of frequencies taken in the subgroup.

We remain with the smooth background of Figure ³ for $m0$⁠. We select groups of 10 consecutive frequencies. The initial group encompasses the frequencies from 0.1 to 1.0 Hz using 0.1 Hz increment. Similarly, we design 15 groups so that the last one is from 14.1 to 15.0 Hz, using an identical increment. The local $Θ$ estimate of the attraction basin size is shown in Figure ⁹, together with the largest tolerable error $Rm0u$⁠. We investigate two directions: the Marmousi one of Figure ⁴ and the large reflector of Figure ⁴. We picture curves and histograms that correspond to sequential or group of frequencies, respectively. Numerical values for the attraction basis are prescribed in Table 2.

The comparisons provided in Figure ⁹ and ⁹ and Table 2 highlight that the size of the attraction basin is dictated by the largest frequency in the group when using the group of frequencies. Therefore, it is better to use the lowest frequency sequentially to have the largest size possible. For instance, using the single frequency 0.1 Hz provides a size of interval one order of magnitude larger than when including frequencies up to 1.0 Hz. Then, at a fixed frequency, taking into account lower frequencies does not affect the size of the interval.

This preference on sequential over multiple frequency is actually advocated by Brenders et al. (2012), motivated, as the frequency progression (Bunks et al., 1995), by the intuition that it reduces the cycle-skipping effect. Here, we justify this choice by a quantitative estimation of the size of the attraction basin. Moreover, we show that once the lowest frequencies are taken individually, incorporating them in higher frequency content does not modify the size of the attraction basin.

Regarding the largest tolerable error on data, Figure ⁹ shows that adding frequency increases the robustness in the direction $um$⁠, except for low frequency. Then, for the direction of the large contrasting object $u3$⁠, we do not observe any improvement, and the allowed distance remains relatively stable, as it was observed in Figure ⁵.

We compute the estimation of $Δm0u$ (using the local $Θ$ estimate) and $Rm0u$ with complex frequencies to see how it affects the basin of attraction and robustness to noise. We consider $m0$ the smooth background velocity. The estimates are shown in Figures ¹⁰ and ¹¹, in which we compare the behaviors on a logarithmic scale for the 2D complex frequency plane. Here, $f$ and $σ$ vary from 0 to 15, with the frequency given in Hz and the Laplace component in $s−1$⁠, excluding the case in which $f$ and $σ$ are zero.

It is clear that the incorporation of a damping coefficient in the frequency is beneficial for the convergence because it increases the size of the attraction basin of several orders of magnitude. We observe

1. 1)

   For a fixed $f$⁠, the interval size increases with increasing damping, with the exception of the case $f=0Hz$ (see Figure ¹⁰ and ¹⁰).

2. 2)

   At fixed damping $σ$⁠, the interval size decreases with increasing $f$⁠.

3. 3)

   Following a selection from the large to small attraction basin gives

   • •

     Initial iterations should use $f=0Hz$ (the largest interval), and a damping coefficient progression from low to high (see Figure ¹⁰ and ¹⁰).

   • •

     Fix $f≠0Hz$ as small as possible, the damping coefficient should evolve from high to low (see Figure ¹⁰ and ¹⁰), and then $f$ increases progressively.

   • •

     Finally, one can take $σ=0s−1$ and $f$ evolves from low to high values.

Regarding the largest tolerable data error $Rm0u$⁠, we see that at fixed $σ$⁠, increasing $f$ allows an increase in the distance between the data and the attainable set. Then, the pattern strongly depends on the type of the search direction geometry because it has already been observed in the previous experiment (and contrary to $Δm0u$⁠). The distance is allowed to increase with $σ$ for the direction $um$⁠, but it decreases in the direction of the strong reflection, $u3$⁠. Also in the direction of the Marmousi reflectors, Figure ¹¹, we observe a band in which the allowed distance decreases, around $σ=1s−1$⁠. On the other hand, for both directions, the case in which $σ=0s−1$ appears the most robust (i.e., where the allowed distance between the data and the attainable set is maximal). In terms of the magnitude, we clearly see that the high contrast object has globally much less allowance, except when the damping is close to zero.

It is important to notice that the use of complex frequencies, especially when $f=0Hz$⁠, improves the convergence by increasing notably the size of the attraction basin. Therefore, it provides a fundamental benefit for the reconstruction, in particular when no information is initially known. Note that one can also use very low frequencies (actually, the magnitude of the estimates for complex frequencies eventually matches the very low-frequency ones; we can see a continuation in the magnitudes). On the other hand, such frequencies are more sensitive to noise and require accuracy in the data. We refer to Shin and Cha (2008), Shin et al. (2010), Petrov and Newman (2014), and Faucher (2017) for illustrations of the reconstruction.

We illustrate the influence of the subsurface geometry with an acoustic FWI experiment, in which we demonstrate the behavior observed in the previous analysis (Figures ⁵ and ¹⁰). In particular, we have shown that a high contrast object complicates the process of reconstruction compared to a Marmousi-like subsurface, by reducing the size of the attraction basin. In addition, we have seen that the very low or complex frequencies, by increasing the size of the basin of attraction, provide benefits on the reconstruction.

Therefore, we perform FWI for 2D geophysical experiments in the acoustic settings. We consider a domain of size $9.2×3km$⁠, which is the size of the models presented in Figure ⁸. We use the algorithm to recover the wavespeed $c$ in the Helmholtz equation 1, from measurements of pressure fields (see equation 3). For the seismic acquisition, we consider 91 sources and 183 receivers per source, which are both located at 20 m depth. The measurements consist of synthetic frequency-domain data with added noise.

We apply the subsurface reconstruction FWI procedure, for two targeted subsurface media of different geometries. We consider the reconstruction of the Marmousi model of Figure ⁸, and the reconstruction of the medium encompassing salt domes of Figure ⁸. For both experiments, we take the same initial (starting) model, which does not assume any information on the subsurface (see Figure ¹²).

The iterative minimization is conducted using frequencies from 2 to 10 Hz, with 1 Hz increment. There are 20 iterations per frequency, which are taken sequentially (from low to high). In Figure ¹², we compare the reconstructions of the Marmousi model and the salt dome, using a similar setup: the same model size, the same frequency progression, and the same initial model with no a priori information (actually, the starting model shown in Figure ¹² is relatively close to the background of the salt medium shown in Figure ⁸ in terms of amplitude).

We observe that only the experiment with the Marmousi medium is able to accurately capture the subsurface structures, in which only the deepest parts are missing in Figure ¹². On the other hand, none of the contrasting salt domes are retrieved, see Figure ¹², with only artifacts in the reconstruction. It demonstrates the results of the previous analysis that, by reducing the size of the basin of attraction, high-contrast objects are more complicated to reconstruct than Marmousi-like media when using a similar frequency range in the data. Note that to obtain a better reconstruction, one could use an initial model less challenging than the one of Figure ¹², i.e., by incorporating a priori information.

It has been shown in Figure ⁵ that low frequencies enlarge the size of the basin. As an alternative, complex frequencies produce the same behavior (see Figure ¹⁰). To demonstrate the validity of these observations, we repeat the experiment of salt dome reconstruction, incorporating low or complex frequencies:

1. 1)

   In Figure ¹³, we show the reconstruction that starts with lower frequencies, using the sequence ${0.5,0.75,1,2,3,4,5,6,7,8,9,10}Hz$⁠, hence incorporating three (unrealistically) low frequencies compared to the reconstruction of Figure ¹².

2. 2)

   In Figure ¹³, we first use the following set of complex frequencies (labeled as $(σ,f)$ in equation 35): ${(7,2),(5,2),(2,2),(1,2),(7,3),(5,3),(2,3),(1,3)}$⁠, with $f$ given in Hz and $σ$ in $s−1$⁠, before using the same sequence as for Figure ¹²: from 2 to 10 Hz (where $σ=0s−1$⁠). Note that our selection of complex frequencies follows our deduction (low to high for $f$⁠, high to low for $σ$⁠).

The FWI algorithm benefits from the incorporation of low or complex frequencies, and it is now able to accurately recover the contrasting objects, as expected. We note that the complex frequency experiment produces less accurate results; in particular, it misses the deepest boundary of the salt dome. However, it does not use a frequency of less than than 2 Hz. Then, both experiments show drastic improvement compared with the results of Figure ¹².

We have used synthetic experiments to validate the basin of attraction estimates and its dependency on the geometry of the target and on the frequency. The use of very low frequencies (as in the reconstruction of Figure ¹³) is efficient but unrealistic in seismic application. The complex frequency alternative is probably a better candidate, but also relies on efficient data preprocessing (e.g., noise removal).

To overcome the issue of low-frequency requirement, the MBTT model decomposition has been introduced. The decomposition is expected to ameliorate the robustness with respect to missing low frequencies, in particular for the recovery of the background velocity. In the next section, we investigate its performance, by the means of attraction basin size, and we quantify the benefits of this approach.

In this section, we apply the optimizability analysis tools to compare the MBTT formulation of the FWI presented in equation 12 (where the unknown model is parameterized by a smooth background $p$ and a data-space reflectivity $s$⁠), with the original FWI problem stated in equation 5 (where the unknown model is the squared slowness $m$⁠). The MBTT approach was originally designed to overcome the phase shift problem that hampers classic FWI, so the optimizability study of this section is expected to quantify the gain — if any — of the MBTT formulation over classic FWI in this respect.

For the perturbation of the background part of the velocity model, we consider the direction of a 1D ramp $u$⁠, given in Figure ¹⁵. In Figure ¹⁶, we show the forward problem (i.e., the solution of the Helmholtz equation at the receivers’ location, deprived from the direct arrivals) corresponding to unperturbed and perturbed models, for a single, centrally located source at 4 Hz. The figure shows the effect of applying the perturbation onto the global model (as used in classic FWI) or onto the background unknown $p$ of the MBTT formulation: When the perturbation is applied to the global model $m$⁠, the phases of the signal are shifted, and the amplitudes are incorrect, but when the perturbation is applied to the background MBTT variable $p$⁠, the phases of the signal are preserved, and only the amplitudes are modified.

For the perturbation direction $us$ associated with the data-space reflectivity $s$⁠, we select a random vector. Note that we expect large basins of attraction with respect to $s$ in the MBTT formulation because the FWI misfit functional is known to be nearly quadratic with respect to the reflectivity, i.e., to the high spatial frequency part, and consequently for the MBTT misfit functional with respect to the data-space reflectivity $s$⁠.

From Figure ¹⁷, we observe that the size of the attraction basins decreases with frequency, and it is only slightly larger when the perturbation is applied onto $p$ (MBTT) instead of $m$ (classic FWI). The gain appears minimal, and it is not up to the expectations raised as the MBTT parameterization allows us to overcome the phase-shift problem, as shown in Figure ¹⁶. Similarly, the attraction basin for $s$ is relatively small compared to what is expected. However, we shall recall that these are only local $Θ$ estimates that can represent a coarse approximation. Therefore, we compute exact $Θ$ estimates in the following, to investigate if those preliminary observations still hold.

To compute the exact $Θ$ and $RG$ estimates, the perturbation of the model is either applied onto the global model $m$ (FWI) or onto the parameters $p$⁠, $s$ of the MBTT decomposition. Therefore, the deflection and global radius of the curvature maps are computed between the following points, for all $−Δ≤t,t′≤Δ$⁠:

1. 1)

   FWI (attraction basin for $m$⁠): $F(m0+tu)$ and $F(m0+t′u)$⁠;

2. 2)

   MBTT (attraction basin for $p$⁠): $F(p0+tu,s0)$ and $F(p0+t′u,s0)$⁠;

3. 3)

   MBTT (attraction basin for $s$⁠): $F(p0,s0+tus)$ and $F(p0,s0+t′us)$⁠.

In Figures ¹⁸ and ¹⁹, we show the deflection and global radius of curvature maps at 7 Hz, for $m$⁠, $p$⁠, and $s$⁠. The corresponding values for the estimates are provided in Table 3, in which we also give the results for the frequency 4 Hz.

We now observe a clear difference between the MBTT parameterization and the global model one, and it shows the gain (in terms of attraction basin size) of the former approach (see Table 3). In the context of global model FWI (Figures ¹⁷ and ¹⁸) the local and exact estimates are of the same size, as observed in the above section. In sight of the upper bound equation 32, on which the local $Θ$ estimate is based, one can think that the FWI formulation corresponds to the “worst” situation, in which the image of a segment in the background space is a curve close to an arc of circle in the data space. However, in the MBTT formulation (Figures ¹⁷ and ¹⁸), the exact $Θ$ estimate is approximately 10 times larger than its local $Θ$ estimate. Therefore, the exact estimates are needed for a proper comparison of the FWI and MBTT approaches. Then, we observe that the values of deflection are much lower when the background perturbation $u$ is applied to $p$ (MBTT) rather than to $m$ (FWI): Few portions reach $π/2$ with MBTT (Figure ¹⁸), whereas it is rapidly attained with classic FWI (Figure ¹⁸). Consequently, the MBTT formulation produces larger $Θ$-attraction basins than the standard FWI formulation, roughly by a factor of 10 (Table 3). Regarding the global radius of curvature $RG$⁠, it is strictly positive for the complete investigated interval in the case of MBTT, whereas, in classic FWI, it rapidly decreases to zero outside the diagonal. Therefore, the $RG$-attraction basins are much larger for MBTT, by a factor of 4–8 (Table 3). In addition, the amplitude of $RG$⁠, whose minimum over the attraction basin gives a tolerable error level, is, near the diagonal, larger for FWI than for MBTT.

Regarding the data-space reflectivity $s$⁠, Figure ¹⁹, the exact $Θ$-attraction basin is particularly large (23–54 times the norm of $s0$ depending on frequency), and it is $105$ times larger than its local estimate. We observe the expected behavior because the forward map $F$ is nearly linear with respect to $s$⁠, hence resulting in large attraction basins. It also confirms the necessity of exact estimates for accurate characterization of attraction basins in the context of MBTT parameterization.

To summarize, MBTT significantly extends the size of attraction basins with respect to background perturbations, explaining its success with relatively high-frequency content in the data (Chavent et al., 2015). However, it costs a reduction in the admissible error level in the data.

The previous comparisons have been performed on a nominal model that has a reflectivity level $β$ of 1% (see equation 38). We repeat the same computations, at 7 Hz, changing only the reflectivity of the model, which is fixed to a low level, 0.2% in Figure ²⁰ (close to the Born approximation), and to a stronger value, 5% (close to Marmousi reflectivity) in Figure ²¹. Table 4 summarizes the sizes and corresponding tolerable error of the attraction basins.

Increasing the reflectivity increases the relative importance of multiple reflections — whose phases are not controlled by the migration/demigration included in the MBTT approach — with respect to primary reflections, whose phases are controlled. It is then expected that the performance of MBTT, when applied to a full-wave model as is the case here, will deteriorate when the reflectivity level increases. This is confirmed in Table 4: The MBTT attraction basin for a reflectivity level of 5% is three times smaller than at 1%, which itself is two times smaller than at 0.2%. On the other hand, the FWI attraction basin remains essentially unaffected by the reflectivity level, and the MBTT basins remain always larger, by a factor of 4–15 with respect to the $Θ$ estimates and a factor of 8 for the $RG$ estimates.

Here, we illustrate an important condition that has to be satisfied for the MBTT decomposition to work. When a background perturbation is applied to $p0$⁠, whereas $s0$ is maintained fixed, the phase shifts of the reflections are controlled, provided these reflections are generated only by $s0$⁠. Any backscattered energy generated by the updated background $p0+tu$ — or worse, by $p0$ itself — will lead to reflections whose phase is not controlled by the migration-demigration process included in MBTT, thus annihilating its benefits. This appears already in Figure ¹⁸, in which one sees that the deflection reaches $π/2$ (respectively, the global radius of curvature equates zero) much closer to the diagonal at the top-right corner than at the bottom-left one. This zone (upper right corner of Figure ¹⁸) corresponds to large positive values of $t$⁠, $t′$⁠, which lead to velocity models (⁠$p0+tu$ and $p0+t′u$⁠) made of strongly increasing background, consequently generating backscattered energy. At the opposite corner of the diagrams (bottom left), which corresponds to mildly increasing velocities, the amount of scattered energy is smaller and the limiting values are further apart from the diagonal. We have considered for simplicity only attraction basins that are centered at the nominal model $m0$⁠, but on most of the deflection and curvature maps, the attraction basin could, in fact, be extended on the side of negative $t$ (where backgrounds are less steep).

In the case in which the unperturbed background scatters back some energy, as the one of Figure ²², the situation is even worse, as illustrated in the corresponding deflection and curvature maps of Figure ²², computed at 4 Hz. These maps show no enlargement of the attraction basin in favor of MBTT; in fact, it is even smaller than for FWI. That is why it is fundamental to ensure the smoothness of the background in the MBTT decomposition by using an adapted parameterization, e.g., spline functions are used by Chavent et al. (2015) to represent $p$⁠.

We illustrate in this section the behavior of the FWI and MBTT misfit functions in the model space direction $u$ at the nominal model $m0=m(p0,s0)$ for some seismic data $d$⁠, where $p0$ is the smooth background of Figure ³. The cost functions are, respectively, $J(m0+tu)$ and $J(p0+tu,s0)$ for FWI and MBTT, according to equations 5 and 12. This will allow us to verify the theoretical analysis: If the minimum of the misfit function is smaller than the tolerable error $RGu$ associated to an attraction basin, the misfit functional should not have any local minimum over the basin. However, when the data exhibit an error larger than the tolerable error, local minima may but do not necessarily appear.

When $η$ increases from zero to one, the distance of $d$ to the attainable set increases from zero to $‖dm−d0‖$⁠, which is larger than the tolerable error associated to the attraction basin of interest (FWI or MBTT).

The computations are done for the nominal model with a 1% reflectivity level, and the background perturbation direction $u$ of Figure ¹⁵. Figure ²³ shows the cost functions for FWI and MBTT at 7 Hz for different values of $η$⁠. We also show by vertical lines the exact $RG$-attraction basin associated to zero tolerable error, and we show by horizontal lines the levels $1/2RG2$ corresponding to the tolerable error level associated to the exact $Θ$-attraction basin. In Figure ²³, we observe that the FWI and MBTT attraction basins are essentially independent of the error level $η$ in the data. The size of the FWI $RG$-attraction basins is slightly larger than the one anticipated by $Δm0u$⁠. When $η=0$⁠, $d=d0$ is attainable, with a zero error level, and one sees that the first local minimum is just after the right limit of the exact $RG$-attraction basin; hence, it is consistent with our analysis. Note that the tolerable error level $RG,p0u$ on $d$ for the $Θ$-attraction basin is attained at $η=0.1$⁠. Clearly, when one moves away from the nominal value $t=0$⁠, the local minima appear much earlier with FWI than with MBTT and the size of intervals in which there is no local minimum is well anticipated by the estimates $Δu$⁠. In the numerical setting considered here, the MBTT approach appears to be robust with respect to data error levels much larger than the tolerable level because we do not see the local minimum appearing inside the attraction basin.

The numerical experiments lead to the following recommendations, for implementing the MBTT reformulation of FWI to obtain large attraction basins for the background reconstruction.

First, the background space for $p$ has to be chosen smooth enough and the reflectivity level $β=‖r‖/‖p‖$ large enough so that the energy backscattered by $p$ is negligible compared to that backscattered by $r$⁠. On the other hand, the reflectivity level $β$ has to be chosen small enough so that the importance of multiples remains small compared to primary reflections in the synthetics. This can be achieved by using an adapted parameterization for the background $p$ (e.g., spline functions) and by tuning the reflectivity level $β$ to adjust the migration weight $W$ in equation 37. We are reminded that the relevance of the chosen settings can be checked before performing any iterative optimization algorithm, by computing the size of the $Θ$ and $RG$-attraction basins and the associated tolerable error. In addition, we have noted that the attraction basins tend to be larger in the direction of slowly increasing (with depth) backgrounds, which suggests underestimating the initial slope of the background velocity with respect to depth.

The dependency of the size of attraction basins with the reflectivity level suggests the following organization for MBTT waveform inversion. First, we need to scale down the given seismic data so that they can be resimulated by a reflectivity level $β$ of approximately 0.2%–1.0%, which amounts to bring the forward model close to the Born approximation. Second, we perform alternate minimization with respect to $s$ and $p$⁠, starting with $s$ first. Third, one can rescale (up) the processed seismic data to get closer to the original given data. Then, we pursue the alternate updates of $s$ and $p$⁠. Eventually, we continue until the processed and given original data coincide. Note that all of the difficulties related to the choice of the background and the reflectivity level disappear if one uses a linearized model (the Born approximation) as the forward model.

We have first studied the classic FWI problem for the acoustic wave equation in the time-frequency domain, and we computed directional attraction basins that give quantitative insight into optimizability of FWI. Unsurprisingly, the size of attraction basins becomes smaller when the frequency increases. In our tests, it is only of $±2%$ of the norm of the background at 4 Hz and $±1%$ at 7 Hz. Clearly, the use of complex frequencies including damping enlarges the attraction basin. However, the largest attraction basins obtained for low (or complex) frequency usually go with smaller tolerable error levels. Therefore, low and complex frequencies suffer more from noise.

We have observed that solving the problem for a sequence of increasing single frequency data is better in terms of the attraction basin than using an increasing sequence of frequency-bandwidth data. As expected, media with salt domes produce smaller attraction basins, and hence they are harder to recover.

We have then studied the MBTT approach, which uses a background/data-space reflectivity parameterization of the velocity model. In a direction of background perturbation, the largest attraction basin for MBTT shows an increase in size by a factor of 4–15 compared to FWI, and it is only barely reduced when the frequency increases. Correlatively, the tolerable error for MBTT is divided approximately by 10 compared to FWI. Nonetheless, the comparison of misfit functions has shown that the local minimum does not appear when the tolerable error is exceeded. In the direction of data-space reflectivity perturbations, in which the forward map is linear up to the multiple reflections, the attraction basin is very large, more than $±20$ times the norm of the nominal data-space reflectivity. Furthermore, we have identified the critical parameters for an efficient implementation of MBTT: the background smoothness and the reflectivity level.

In this paper, we provide an a priori analysis of the optimizability of nonlinear least-squares minimization problems. For this purpose, we compute the size of the attraction basin around nominal models and the associated tolerable error level, such that the misfit functional is guaranteed to have only one (hence global) minimum in the attraction basin, provided the noise on the data is below the tolerable error level. These estimates are defined a priori, and they depend only on the first- and second-order directional derivatives of the forward problem.

A first application is carried out with the seismic FWI problem, in which the misfit functional is known to have local minima in the directions associated with low spatial frequencies perturbation (background velocity). Yet, note that our methodology is totally applicable for other least-squares minimization problems. Our findings confirm the inherent difficulty of FWI, which requires unrealistically low-frequency data for the determination of the full velocity model and which we illustrate with experiments of reconstruction.

To see if these limitations could be overcome, the methodology is then applied to compare with the data-space reflectivity MBTT reformulation, originally designed to alleviate the phase shift problem. The results show a substantial increase of the attraction basins for MBTT, which can indeed overcome the requirement for low frequencies in the data.

As a side product of this analysis, we finally provide some guidelines for an efficient implementation of the MBTT reformulation of FWI.

The authors would like to thank the editor, associate editor, and the three anonymous referees for their valuable comments that have contributed to improve the quality of the paper. F. Faucher acknowledges funding from the Austrian Science Fund (FWF) under the Lise Meitner fellowship M 2791-N. The research of F. Faucher is also supported by the Inria-TOTAL strategic action DIP. The research of H. Barucq has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement no. 777778.

Data associated with this research are confidential and cannot be released.