{'fraction_non_alphanumeric': 0.030383424925778536, 'fraction_numerical': 0.020339839555302885, 'mean_word_length': 4.560765718299965, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extended full waveform inversion in the time domain by the augmented Lagrangian method

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


Introduction
Full-waveform inversion (FWI) is the state-of-the-art data-fitting based imaging method for subsurface parameters estimation. The parameters of interest appear as coefficients in system of hyperbolic partial differential equations (PDE) which serves as the forward modeling engine for simulating seismic waves. Sampling these waves at receiver locations gives us synthetic seismograms. FWI tries to find the set of parameters from which we can compute synthetic seismograms which best fit (in the least-squares sense) the observed seismograms (data) (Pica et al., 1990;Virieux and Operto, 2009). The forward problem is nonlinear thus iterative linearization methods are natural choice to solve the problem (Tarantola, 2005).
Traditional formulation of FWI eliminates the intermediate variable (wavefields) from the equations and represents observed data directly in terms of the Green's functions which are nonlinearly related to the model parameters. The standard gradient (steepest descent) method solves this nonlinear inverse problem iteratively. At each iteration, the parameters are updated along the steepest-descent direction. This direction is efficiently determined by using adjointstate method (Plessix, 2006) that is a simple zero-lag cross-correlation between two calculated wavefields: one is the A PREPRINT forward wavefield in the current model and the other is a (backward) wavefield obtained by back-propagation of the data residuals acting as if they were sources at receiver locations. A more accurate search direction is obtained by using Newton or quasi-Newton methods by including curvature information via the Hessian matrix or an approximation of it (Métivier et al., 2017). A simple while effective approximation of the Hessian matrix can also be obtained by using only its diagonal terms that are calculated by zero-lag correlation of the forward wavefield (Shin et al., 2001). Calculated step direction is then scaled by an appropriate step length that can be chosen using a line search method, such as the backtracking line search, and added to the previous model to get the next iterate. This process is repeated until the convergence is achieved. In this way, the quality of the search direction is the main factor governing the algorithm convergence and it is a directly influenced by the quality of forward and backward wavefields.
In contrast to the FWI, extended FWI brings flexibilities to fit the data by softening the wave equation constraint at least at early iterations (Aghamiry et al., 2019;Huang et al., 2018;van den Berg and Abubakar, 2001;van Leeuwen and Herrmann, 2013). The associated objective function is bivariate and besides model parameters we need to also manipulate data-assimilated wavefields (Aghamiry et al., 2020), specific wavefields that are obtained by solving the data equation and the wave equation simultaneously. The objective function is biconvex and the augmented Lagrangian (AL) method provides an excellent framework for its minimization with a minimal sensitivity to the penalty parameter choice (Aghamiry et al., 2019;Gabay and Mercier, 1976;Nocedal and Wright, 2006). However, there are two main obstacles that limit the usage of this method in the time domain (Aghamiry et al., 2020;Rizzuti et al., 2019;van Leeuwen, 2019;Wang et al., 2017). The first is reconstruction of data-assimilated wavefields and the second is the need to store the Lagrange multipliers, which are of the size of the wavefields and thus infeasible for the field-scale problems.
In this paper, we provide a new formulation for solving the AL-based extended FWI. An equivalent formulation of the iteratively-refined wavefield reconstruction inversion (Aghamiry et al., 2019) is derived in the data space which require manipulating Lagrange multipliers of the data size. Furthermore, it replaces the wavefield reconstruction problem in the full space with a residual data reconstruction problem in the data space. The resulting reconstruction system is of the data size and thus more tractable compared with the wavefield reconstruction in the original formulation. Based on the fact that the AL algorithm is convergent for a wide range of penalty parameters (Aghamiry et al., 2019;Gabay and Mercier, 1976;Nocedal and Wright, 2006), we derive a simplified algorithm which only requires four wave equation solves at each iteration which can be carried out by usual time-stepping. An interesting observation is that the AL-based extended FWI updates the model parameters via a normalized adjoint-state formulation similar to the quasi-Newton algorithm for the classical FWI. The only difference is in the forward and backward wavefields which are cross-correlated. This can be thought of as a generalized algorithm which can be applied to both the extended and reduced FWI and in both the time and frequency domains.
The contents of this paper are organized as follows: We begin with the general form of the FWI problem as a partialdifferential equation (PDE) constrained data-fitting problem. We then review the two approaches of classical FWI and the extended FWI for solving this optimization problem. A brief review of the augmented Lagrangian (AL) method is presented with its advantages and limitations for time domain FWI. Then we derive an algorithm for ALbased extended FWI and show how it can be efficiently applied in the time domain. Finally, we will demonstrate the efficiency of the proposed algorithm by using Marmousi II and BP 2004 benchmark velocity models.

Theory
FWI seeks the subsurface parameters by solving the following PDE constrained optimization problem (Aghamiry et al., 2019;van Leeuwen and Herrmann, 2013): where · 2 denotes the Euclidean norm, u ∈ R N ×1 is the wavefield, d ∈ R M ×1 is the recorded seismic data generated by the sampling operator P ∈ R M ×N (with M N ), A(m) ∈ R N ×N is the PDE operator, and b ∈ R N ×1 is the source term. More specifically, N = N x × N z × N t and M = N t × N r with N x and N z the number of spatial grid samples of the model parameter m (square of slowness) in horizontal and vertical dimensions, N t the number of time samples, and N r is the number of receivers.

FWI
Eliminating the wavefield from equation 1 leads to the conventional (reduced) formulation of FWI as (Pratt et al., 1998) where G(m) = PA(m) −1 denotes the forward modeling operator. This objective function is very hard to solve due to the presence of the inverse operator A(m) −1 (which includes the Green's functions as its columns). The standard gradient-based iterative methods solve this nonlinear inverse problem iteratively (Pica et al., 1990;Pratt et al., 1998;Tarantola, 1988;Virieux and Operto, 2009), in which at iteration kth the parameters are updated as where α k is a step length that can be chosen using a line search method, such as the backtracking line search, and δm k is a step direction. This step direction is selected as opposite gradient of the objective function (g k , in the steepest descent method) or determined by solving a Newton or quasi-Newton system: where H k is the Hessian matrix or an approximation of it. Both the gradient and the Hessian are evaluated at the current model m k . The gradient can be computed efficiently by using the adjoint state method (Pica et al., 1990;Plessix, 2006).
whereü r k (t) is the second time derivative of the forward wavefield induced by the source and v r k (t) is the adjoint wavefield obtained by back-propagation of the data residuals δd r k (which acts as a source emitting in a time-reversed manner).
(6) In equation 5, both wavefields are propagated in the current model m k and they are recorded over the time interval (0, τ ). Also, • denotes the element-wise product operator.
A simple while effective approximation of the Hessian or the pseudo-Hessian includes only the diagonal terms (Shin et al., 2001), which gives where the division and the square-power are element-wise operations. This indeed is the zero-lag crosscorrelation image between the forward and adjoint wavefields (imaging condition) normalized by the square of the source illumination strength (Claerbout, 1971). A more accurate Hessian matrix can also be implement to improve the model update by solving equation 4 via an iterative algorithm while computing Hessian-vector products by second-order adjoint state formulas (Métivier et al., 2017).
However, FWI requires us to provided an initial model of sufficient quality for convergence to the desired solution (Virieux and Operto, 2009). The next section provides a remedy for this issue using extended FWI. In the formulas above, superscript r shows that the variable corresponds to the reduce FWI. We will use superscript e to refer to the extended FWI.

Extended FWI
Extended FWI methods have been developed to decrease the nonlinearity of the objective function and hence increase the robustness of local optimization methods to the initial model (Aghamiry et al., 2019;Huang et al., 2018;van den Berg and Abubakar, 2001;van Leeuwen and Herrmann, 2013). The simple form of the extended formulation can be obtained by the penalty formulation of the original equation 1: where µ > 0 is the penalty parameter. This extended formulation can be thought of as an "error-in-variable" method which accounts for errors in the wave equation, and it reduces to the FWI objective in equation 2 as µ → ∞. Equation 8 is biconvex and it is usually solved by alternating minimization (Beck, 2017). This algorithm however suffers from two major drawbacks which limit its application in a time domain practical setting. The first is manipulating the A PREPRINT u-subproblem (computation and storage of the wavefield for a given model) and the other is the issue related to the penalty parameter choice.
The issue of penalty parameter can be overcome by adding a Lagrangian term of the form λ T (A(m)u e − b) to the penalty function (Aghamiry et al., 2019;Nocedal and Wright, 2006): where superscript T denotes adjoint operator, Penalty denotes the penalty objective function in equation 8 and λ is the Lagrange multiplier for the PDE constraint. The alternating direction method of multipliers (Aghamiry et al., 2019;Boyd et al., 2010;Gabay and Mercier, 1976;Goldstein and Osher, 2009) solves the augmented Lagrangian in equation 9 iteratively by choosing u and λ and then minimizing the objective as a function of m. The resulting m is used to choose a new u and λ, and the process repeats until the convergent is achieved.
beginning with initial u 0 (which is built from an initial model) and λ 0 = 0, where A k ≡ A(m k ) and k is the iteration count. Note that Aghamiry et al. (2019) first optimize in the wavefield u e and then in the model parameters m but for the method to develop in this paper the permuted order in equation 10 is more convenient. This does not affect the final solution because the order of the optimization in biconvex problems can be permuted. The iteration in equation 10 without the "add-back-the residual" term reduces to the wavefield reconstruction inversion (WRI) (van Leeuwen and Herrmann, 2013). This term accounts for the Lagrangian and besides increasing the convergence rate, reduces the sensitivity of the algorithm to the penalty parameter choice by "adding-back the source residual". The advantage of the augmented Lagrangian method is thus to use a small µ (to make the optimization well posed) and "add-backthe-residual" (Goldstein and Osher, 2009). Aghamiry et al. (2019) consider this iteration while solving the feasibility form of equation 1 in which the data are also updated by adding back the data residual.
beginning with b 0 = b and d 0 = d which have the role of scaled Lagrange multipliers. This algorithm has been successfully applied in the frequency-domain FWI where the most challenging subproblem, equation 10b, has a size of N x × N z and thus can be solved directly by sparse solvers. The Lagrangian multiplier of size N x × N z can also be conveniently stored during the iterations. In the time domain, however, solving equation 10b is more challenging because the size of the system would be prohibitively large (N x × N z × N t ) due to the extra temporal dimension. Furthermore, b k of the size N x × N z × N t would be dense; which causes serious storage problems. Several attempts have been make by different authors to solve these issues (Aghamiry et al., 2020;Rizzuti et al., 2019;van Leeuwen, 2019;Wang et al., 2017). In the next section, we propose an algorithm for efficient treatment of these issues which performs in the data space instead of the full wavefield space.

Data Reconstruction Inversion
The general form of data-assimilated wavefield reconstruction, equation 10b, is Minimization of equation 12 gives u e as the solution of the following system: This system can be re-expressed by decoupling Au e from other terms, making it suitable for time-stepping: However, since the unknown wavefield u e appears in both sides of this equation, researchers tried different approaches to solve it. For example, Wang et al. (2016) replaced the wavefield at the right-hand-side with the wavefield obtained by the reduced FWI. Later, Aghamiry et al. (2020) developed an algorithm which solves equation 14 iteratively, beginning with the wavefield constructed by Wang et al. (2016).
In this paper, we solve equation 14 in a two step procedure. Noting that u e appears in the right-hand side only as (WRI predicted data) Pu e thus if we have an estimate of Pu e then u e can be determined by solving equation 14. Accordingly, we rewrite this equation as are respectively the source residual and data residual (of WRI). Note that the source residual is zero for the reduced FWI. Up to now we have reduced the optimization dimension from N (directly searching for u e in equation 13 or 14) to M (searching for the residual vector δd e in equation 15). However, the question is now "what is the optimization problem for δd e ?" By eliminating u e from equation 16 we get the following relation between the two residual vectors δb e and δd e : Gδb e + δd e = δd r .
where δd r is the data residual of the FWI, equation 6.
Furthermore, comparing equations 14 and 16 shows that Substituting δb e from equation 18 into equation 17 immediately gives an explicit M × M system for δd e in the data space: This equation clearly shows that the FWI data residual, δd r , is a blurred version of the WRI data residual, δd e , in which the blurring kernel consists of the correlation of receiver-side Green's functions. The solution of equation 19 is also the minimizer of the following quadratic objective function: Plugging δb e and δd e from equations 18 and 19 into the primal objective in equation 12 gives an equivalent objective for the extended FWI which is a weighted norm of the FWI data residual: where y 2 Q(m) = y T Q(m)y and Q(m) is the coefficient matrix in equation 19. This equivalence has also been confirmed differently by van Leeuwen (2019) and Symes (2020).

IR-WRI Iteration in Data Space
Using the data space formulation given above the IR-WRI iteration given in equation 11 can be performed more efficiently in the data space. From equations 11c, 11d and 16 we have that Furthermore, from equations 18 and 19 we have Accordingly the wavefield u e k in equation 10a can be written as (see equation 15) (while replacing G k−1 with its update G k ) based on which the model update in equation 10a simplifies to From equation 20, we see that δd e k is the minimizer of the following objective function Differentiating this objective with respect to δd e gives Substituting the values of d k−1 and b k−1 from equations 22 and 23 into equation 31 gives Note that, in equation 34, we replaced G k−1 with its update G k . Setting the gradient, equation 35, equal to zero and solving for δd e gives the exact minimizer, δd e k , of the objective It is seen that in the augmented Lagrangian method the data residual vector is updated via a Newton like algorithm, equation 36, in which the search direction is defined by the FWI data residual δd r k . Q k has the role of Hessian matrix. After introducing the auxiliary variable y k = δd e k /µ and δy k = y k − y k−1 , equations 28 and 29 simplify to the following iteration where G T k = A −T k P T . As seen, the wavefield reconstruction step in equation 11 is replaced by a residual data reconstruction step, δy k in equation 37, and thus we call this new iteration data reconstruction inversion (DRI). DRI only needs to store the dual variable y k lying in the data space.

DRI with Gradient Descent
A fundamental property of all augmented Lagrangian based algorithms is that they are convergent for an approximate solution of the subproblems (Aghamiry et al., 2019;Boyd et al., 2010;Gabay and Mercier, 1976;Gholami, 2015;Gholami et al., 2018;Yin and Osher, 2013). This motivated us to perform the minimization of objective f in equation 30 approximately via a gradient step rather than full minimization given in equation 36.
where α k is a step length and ∇f is the gradient of the objective function determined in equation 35, evaluated at δd e k−1 . From equation 35 we get that ∇f (δd e k−1 ) = −δd r k .
(39) Accordingly, δy k in equation 37 becomes δy k ≈ α k δd r k , (40) which (after replacing y k and y k−1 by y k /α k and y k−1 /α k ) simplifies the DRI iteration as

Step length determination
In the DRI algorithm, equation 41, the model parameters are updated by assuming a large penalty parameter and thus more weight given to minimization of the source residuals rather than the data residual. Accordingly, we determine the value of α k by minimizing the data residuals: giving that where q k = G k G T k δd r k .

Update of the Model Parameters
We consider constant density acoustic wave equation as the forward modeling operator, given by where m is the square of slowness, u(t) is the discretized wavefield andü(t) denotes its second time derivative, ∇ 2 is the Laplacian operator with respect to the spatial coordinates, and b(t) is the source term. Then, minimization of the source residual, equation 41, leads to where u e k (t) obeys the following equation (see equation Using the expression for ∇ 2 u e k (t), equation 45 can also be written as where in which v e k (t) is the backward-wavefield obtained by back-propagation of the adjoint source y k + δd r k . As seen, the new algorithm simply adds back the residual to improve the estimate.

Computational Cost
The computational burden of DRI algorithm lies in the simulation of propagating wavefields using a finite difference method. More precisely, two fields u e k and v e k must be simulated and need to be available at the same time level to compute the integral in equation 49. One can compute these integrations by first computing the forward wavefield u r k (solution of one forward problem) and storing it in the computer memory. This first step also delivers the data residual δd r k . Then the residual field δu k is computed by back-propagation of the data residual followed by forward propagation of the resulting backward field (solution of one forward plus one backward problems), during which the step length α k , equation 43, and the wavefield u e k = u r k + α k δu k are computed. Finally, the backward field v e k is computed by back-propagating the adjoint source y k + δd r k (solution of one backward problem) and the imaging condition in equation 49 is computed simultaneously. Therefore, this approach requires one to solve 2 + 2 wave equations at each iteration. In practice, storing the wavefield is memory intensive. To mitigate this issue, one can use the compression and reconstruction methods to decrease demanded memory (e.g., Boehm et al., 2016;Gholami et al., 2019;Zand et al., 2019). An alternative method, which is less memory demanding but adds the cost of extra wave equation solves, is to store only the spatial boundary conditions and final time slice and then recomputing the wavefields (Gauthier et al., 1986;Symes, 2007).

Camembert model test
In order to show the robustness of the proposed method with respect to the initial model, we consider the simple while challenging "Camembert" model (Gauthier et al., 1986). The subsurface model contains a circular anomaly of velocity 4.6 km/s embedded in a homogeneous background of velocity 4.0 km/s (Figure 1a). The dimensions of the model are 4.8 km in distance and 6 km in depth, and the grid spacing is 35.5 m. The crosshole acquisition consists of 14 equally spaced sources of 10 Hz Ricker wavelet in the left side of the model and 170 equally spaced receivers deployed vertically on the opposite side. This model was previously used by Engquist and Yang (2020) for an application of FWI based upon an optimal-transport distance. Due to the large diameter of the anomaly the inverse problem is nonlinear and thus the initial model plays an important role in the convergence of the local optimization algorithms. We start the inversion from the homogeneous background model to meet the nonlinear regime of the FWI and perform one cycle of inversion with the 10 Hz Ricker wavelet. The results obtained using FWI and the DRI are presented in Figures 1c and 1d, respectively. Figure 2 shows the trajectory of the misfit function versus iteration for both methods. The result of FWI is quite different from the true one, which indicates that the method suffers from severe cycle skipping due to the rough initial model and 10 Hz frequency inversion. In contrast, the DRI mitigates efficiently the cycle skipping. This is clearly seen from the reconstructed model which is quite close to the ground truth.

Checkerboard model test
We design a checkerboard model composed of a homogeneous background of velocity 1.5 km/s and checkerboard velocity perturbations of velocity 4.0 km/s and dimension 200 m by 200 m (Figure 3a). We use a full acquisition with 92 equally-spaced sources placed all around the domain with 50 m distance from the model boundary. Each source is associated with a total of 100 equally spaced receivers placed on the model boundary. The model is discretized over a 101 by 101 grid with a spatial step of 20 m. The source wavelet is a band-pass filtered Ricker wavelet between 2.5 Hz and 5 Hz with a peak frequency of 3.75 Hz (Figure 4). The initial model is the background velocity. Figures  5a and 5b show two data-sets computed with the exact velocity model and with the initial homogeneous-velocity model. We compare the inversion results obtained performing 300 iterations ( Figure 6) for both FWI and DRI and the corresponding model estimates are shown in Figures 3b and 3c. The final data residuals corresponding to both methods are also shown in Figures 5c and 5d. For this test, again, unlike FWI, the DRI reconstructed a model which is quite close to the ground truth.

Marmousi II Model test
We now consider a more realistic multisource experiment using the Marmousi II velocity model (Figure 7). The model is 17 km long and 3.5 km deep. The fixed spread acquisition consists of equally-spaced 91 sources and equally-spaced 273 receivers on top side of the model. We build reference data using a highpass filter minimum-phase Ricker wavelet with no signal below 2.5 Hz. The total recording time is 7.5 seconds with a sample internal 3 ms. The initial model is a 1D velocity model linearly increasing in depth from 1.5 km/s to 4.5 km/s. Figure 8 shows the residual gather for the source at distance 8.5 km. Obviously classical FWI will stuck in a local optimum due to severe cycle skipping caused by the 1D initial model and the absence of low-frequencies in the data (see Aghamiry et al., 2019, Their Figure  6) thus we don't show the corresponding results here. The inversion result of the DRI is shown in Figure 9 which clearly shows that the method mitigated efficiently the cycle skipping. The reconstruction is quite close to the true model in most parts of the model. In order to compare the DRI algorithm with the IR-WRI algorithm we performed both algorithm in the frequency domain to invert the data corresponding to the Marmousi II. The results obtained by both algorithms are shown in Figures 9b-9c and as seen the estimated velocity models are similar which confirm the accuracy of the approximation we considered in DRI algorithm, given in equation 40.

2004 BP Model test
Finally, we assess the DRI method against a target of the challenging 2004 BP salt model when a crude starting model is used (Figure 10). The selected target corresponds to the left part of the 2004 BP salt model and was previously used by Métivier et al. (2016) for an application of FWI based upon an optimal-transport distance and by Aghamiry et al. (2019) for an application of frequency-domain WRI based upon the augmented Lagrangian formulation (IR-WRI). The subsurface model is 16.25 km wide and 5.825 km deep. The fixed spread acquisition consists of equally-spaced 109 sources and equally-spaced 326 receivers on top side of the model. We build reference data using a highpass filter minimum-phase Ricker wavelet with no signal below 2.5 Hz. The total recording time is 7.5 seconds with a sample internal 2 ms. We used a smoothed version of the true velocity model without the salt as initial model ( Figure  10). A comparison between a common-shot gather computed in the true and initial models shows severe traveltime mismatches ( Figure 11) that drive the classical FWI to a local minimum (see Aghamiry et al., 2019, Their Figure  12c). The DRI, however, mitigated efficiently the cycle skipping as it can be seen from the reconstructed model shown in Figure 12. As seen from this figure, the precise geometry of the salt body and the subsalt low-velocity anomalies are mostly recovered with a quite high resolution.

Conclusions
Time domain full waveform inversion (FWI) is appealing since the wave equation can be solved efficiently by explicit time-stepping but this method require sufficiently accurate initial models which may not be available in practice. Distance (km) Extended FWI with wavefield reconstruction is robust with respect to the initial model and it can be solved with a minimal sensitivity to the penalty parameter choice when using the augmented Lagrangian solver. However, its application in the time domain is a challenging task due to the absent of explicit time-stepping and because computation and memory burdens grow quickly as the size of model increases due to the inherent high-dimensional nature of the state variable. In this paper, we proposed a new augmented Lagrangian algorithm for this task which performs in the data space, a significantly lower dimensional subspace than the dimensionality of the original full space. In this algorithm, the wavefield reconstruction appears as a residual data reconstruction and the model parameters are updated by a quasi-Newton like algorithm with a search direction computed by a generalized adjoint state method. This latter is computed by zero-lag cross-correlation between the data-assimilated wavefields and adjoint fields coming from backpropagation of all previous data residuals added back to the current residual vector. Approximate update of the residual reconstruction by a gradient descent method resulted to an efficient algorithm which only requires four wave equation solves (two forward solves plus two backward solves) at each iteration. Different numerical tests including Marmousi II and 2004 BP velocity models confirm the efficiency of this algorithm for time-domain applications similar to the classical FWI while benefiting from the excellent robustness of the extended FWI.

Acknowledgments
This study was partially funded by the WIND consortium (https://www.geoazur.fr/WIND), sponsored by Chevron, Shell, and Total. This study was granted access to the HPC resources of SIGAMM infrastructure (http://crimson.oca.eu), hosted by Observatoire de la Côte d'Azur and which is supported by the Provence-Alpes Côte d'Azur region, and the HPC resources of CINES/IDRIS/TGCC under the allocation A0050410596 made by GENCI."

A Appendix A
Consider the acoustic wave equation in a homogeneous 3D medium 1 c 2 ∂ 2 ∂t 2 u(t, x) − ∇ 2 u(t, x) = b(t, x), in which the acoustic velocity c is a constant. The solution of equation 50 can be written in the following integral form: where are the Green's functions and δ denotes the delta function. Writing equation 51 in an operator notation as u = A −1 b then the goal is to build Q = P A −1 A −T P T .
(53) where P is the restriction matrix at the receiver locations. In order to build Q we need to compute the adjoint wavefield for each column of P T followed by computing the forward wavefield while using the ajoint field as the source and then sampling the result at receiver location. Each column of P T can be written as in which x j denotes the location of the jth receiver and τ denotes the time lag. The adjoint wavefield for the delta function in equation 54 is Then applying the forward operator on this adjoint wavefield gives Finally, sampling this wavefield at the ith received location, thus counting for the matrix P , gives Substituting for the Green's functions from equation 51 and performing the integration over t gives In order to compute this integral, we note that for each receiver pair i and j and a point x the Q matrix includes a linear event parallel to the main diagonal described by in which the second term at right hand side determines the amount of shift from the main diagonal. From the backward triangle inequality we have that which implies that the maximum shift is 1 c x i − x j 2 . Summing over all points x thus leads to a symmetric banded matrix with a band centred on the main diagonal. Figure 13 shows the Q matrix for a homogeneous model of velocity 2000 m/s for 5 equally spaced receivers each of 251 time samples. The distance between receivers is 200 m. The band diagonal structure of each block is clearly seen and the bandwidth follows the theoretical finding.