Efficient extended-search space full-waveform inversion with unknown source signatures
Hossein S. Aghamiry, Frichnel W. Mamfoumbi-Ozoumet, Ali Gholami, Stéphane Operto
EE FFICIENT EXTENDED - SEARCH SPACE FULL - WAVEFORMINVERSION WITH UNKNOWN SOURCE SIGNATURES
A P
REPRINT
Hossein S. Aghamiry
University Cote d’Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. [email protected]
Frichnel W. Mamfoumbi-Ozoumet
University Cote d’Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. [email protected]
Ali Gholami
Institute of Geophysics, University of Tehran, Tehran, Iran. [email protected]
St´ephane Operto
University Cote d’Azur - CNRS - IRD - OCA, Geoazur, Valbonne, France. [email protected]
February 9, 2021 A BSTRACT
Full waveform inversion (FWI) requires an accurate estimation of source signatures. Due to thecoupling between the source signatures and the subsurface model, small errors in the former cantranslate into large errors in the latter. When direct methods are used to solve the forward prob-lem, classical frequency-domain FWI efficiently processes multiple sources for source signatureand wavefield estimations once a single Lower-Upper (LU) decomposition of the wave-equationoperator has been performed. However, this efficient FWI formulation is based on the exact solu-tion of the wave equation and hence is highly sensitive to the inaccuracy of the velocity model dueto the cycle skipping pathology. Recent extended-space FWI variants tackle this sensitivity issuethrough a relaxation of the wave equation combined with data assimilation, allowing the wavefieldsto closely match the data from the first inversion iteration. Then, the subsurface parameters are up-dated by minimizing the wave-equation violations. When the wavefields and the source signaturesare jointly estimated with this approach, the extended wave equation operator becomes source de-pendent, hence making direct methods and, to a lesser extent, block iterative methods ineffective.In this paper, we propose a simple method to bypass this issue and estimate source signatures ef-ficiently during extended FWI. The proposed method replaces each source with a blended sourceduring each data-assimilated wavefield reconstruction to make the extended wave equation operatorsource independent. Besides computational efficiency, the additional degrees of freedom introducedby spatially distributing the sources allows for a better signature estimation at the physical locationwhen the velocity model is rough. We implement the source signature estimation with a variableprojection method in the recently proposed iteratively-refined wavefield reconstruction inversion(IR-WRI) method. Numerical tests on the Marmousi II and 2004 BP salt synthetic models confirmthe efficiency and the robustness against velocity model errors of the new method compared to thecase where source signatures are known.
Seismic wavefields carry information about subsurface and source, the latter being represented by its location andsignature. In controlled source seismic on which this study is focused, the source locations are generally known ac- a r X i v : . [ m a t h . O C ] F e b R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT curately, while the source signatures are usually unknown and need to be estimated to perform reliable full-waveforminversion (FWI) (Pratt et al., 1998; Tarantola, 1984; Virieux & Operto, 2009). Furthermore, it is well acknowledgedthat the estimation of the source signature is easier in the frequency domain than in the time domain since the time-harmonic wave equation can be solved for each frequency separately (Pratt, 1999; Song et al., 1995). In frequency-domain seismic modeling, direct methods are the most suitable ones to process a large number of sources efficientlyby forward/backward elimination, once a Lower-Upper (LU) decomposition of the so-called impedance matrix hasbeen performed once (Marfurt, 1984). When the size of the problem prevents using a direct solver, iterative methodsspeed-up the processing of multiple right-hand sides with block and recycling methods. Recently, Fang et al. (2018)tackled the source signature estimation problem in the framework of an extended formulation of frequency-domainFWI called wavefield reconstruction inversion, where the recorded data are assimilated during the wavefield recon-struction to closely match the data with inaccurate subsurface models and hence prevent cycle skipping (van Leeuwen& Herrmann, 2013). In the formulation of Fang et al. (2018), the data assimilation makes the extended wave-equationoperator source dependent, and hence the method is expensive since the LU decomposition needs to be performedfor each source when a direct solver is used. Also, block-processing of multiple right-hand sides is not possible any-more with iterative methods. The focus of this paper is to revisit the source signature estimation problem in extendedfrequency-domain FWI such that the forward operator remains source-independent, and hence multiple sources canbe processed efficiently as in the reduced-space FWI formulation.Source signatures may be estimated before FWI or updated jointly with subsurface parameters during FWI iterations.For a fixed velocity model, the source signature estimation can be formulated as a least-squares quadratic data fittingproblem (Pratt, 1999). The closed-form expression of the estimated source signature is given by the zero-lag cross-correlation between the calculated and recorded data, scaled by the auto-correlation of the calculated data. In thisframework, the source signature estimation can be implemented in the classical reduced-space FWI iterations withtwo different approaches: In the first, the source signatures and the subsurface parameters are updated in an alternatingmode, while the second approach enforces the closed-form expression of the estimated source signature as a functionof the subsurface parameters in the objective function (Aravkin & van Leeuwen, 2012; Aravkin et al., 2012; Li et al.,2013) through a variable projection approach (Golub & Pereyra, 2003). Plessix & Cao (2011) review these two formu-lations in the frame of the adjoint-state method and conclude that the variable projection method is more versatile toimplement the source signature estimation problem with specific data weighting, while Rickett (2013) showed that thevariable projection approach was more resilient to phase errors in the wavelet than the alternated optimization. Thissource signature estimation doesn’t introduce significant computational overhead in classical FWI since the gradientof the objective function with respect to the subsurface parameters is computed in the same way whether the sourcesignature is available or estimated on the fly during the FWI iterations (Aravkin & van Leeuwen, 2012; Rickett, 2013).In its more general form, FWI can be cast as a constrained optimization problem that aims to estimate the wavefieldsand the subsurface parameters by fitting the recorded data subject that the wave equation is satisfied (Haber et al.,2000). Regardless of the source signature estimation issue, it is well acknowledged that FWI is highly nonlinear whenthe full search space encompassed by the wavefields and the subsurface parameters is projected onto the parameterspace after elimination of the wavefield variables. This variable elimination, which is performed by forcing the wave-fields to satisfy exactly the wave equation at each FWI iteration, makes FWI prone to cycle skipping as soon as theinitial model is not accurate enough to predict recorded traveltimes with an error smaller than half a period (Virieux &Operto, 2009). To mitigate the cycle skipping issue, some approaches implement the wave equation as a soft constraintwith a penalty method such that the data can be closely matched with inaccurate subsurface models from the earlyFWI iterations (Abubakar et al., 2009; van Leeuwen & Herrmann, 2016, 2013). Then, the subsurface model is updatedby solving an overdetermined quadratic optimization problem, which consists of minimizing the source residuals gen-erated by the wave equation relaxation. In these extended approaches, the wavefields are reconstructed by solving ina least-squares sense an overdetermined linear system gathering the wave equation weighted by the penalty parameterand the observation equation relating the simulated wavefield to the data through a sampling operator. In other words,the wavefields are reconstructed with data assimilation. This approach was called Wavefield Reconstruction Inversion(WRI) by van Leeuwen & Herrmann (2013). A variant of WRI, based upon the method of multipliers or augmentedLagrangian method, was proposed by Aghamiry et al. (2019b) to increase the convergence rate and decrease the sen-sitivity of the algorithm to the relaxation (penalty) parameter choice. The augmented Lagrangian method combinesa penalty method and a Lagrangian method, where the penalty term is used to implement the initial relaxation of theconstraint while the Lagrangian term automatically tunes the sensitivity of the optimization to the constraint in itera-tions through the gradient ascent update of the Lagrange multipliers with the constraint violations. This method wascalled Iteratively-Refined(IR)-WRI, where the prefix IR refers to the iterative defect correction action of the Lagrangemultipliers.WRI was recently extended to jointly estimate the source signatures and the subsurface parameters (Fang et al., 2018).In this approach, the monochromatic data-assimilated wavefield and the signature of the source are gathered in anunknown vector and are estimated in a least-squares sense. Although the wave-equation relaxation increases the ro-bustness of the estimated source signatures against velocity model errors, the method is time-consuming because the2R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT normal system of the overdetermined wavefield-reconstruction problem is source-dependent, hence preventing effi-cient processing of multiple right-hand sides either with direct or iterative solvers. Another variant of WRI withunknown source signatures was proposed by Huang et al. (2018) where WRI is re-parametrized in terms of extendedsources and subsurface parameters. In this approach, the penalization (or annihilator) term is defined as the distancefunction from the real source position (Huang et al., 2018, Their equation 9), which means that the source signatureestimation is implicitly embedded in the extended source reconstruction. One issue with this approach is related to thepresence of the Green functions in the Hessian of the extended source reconstruction subproblem, which makes thenormal system very challenging to solve with a good accuracy (Huang et al., 2018, Their equation 11). Moreover,they update the subsurface parameters with a variable projection method, which precisely requires an accurate solutionof the normal system for the extended sources.In this paper, we implement a fast and robust (against model errors and recorded data inaccuracies) multi-source sig-nature estimation in IR-WRI with a variable projection method, namely, the closed-form expression of the sourcesignature is projected in the wavefield reconstruction subproblem. To achieve the computational efficiency of themulti-source signature and wavefield reconstructions, we reconstruct each individual wavefield with blended sources.This blending makes the normal operator of the wavefield reconstruction subproblems source independent and henceamenable to efficient multi-source processing. Although we use source blending, we stress that we estimate onewavefield per physical source thanks to the assimilation of the source-dependent recorded data in the right-hand sideof the normal system satisfied by the reconstructed wavefield. However, the source blending implies that, for eachreconstructed wavefield, the source signature is a vector of dimension equal to the number of individual sources inthe blended source. When the velocity model is accurate, each entry of the signature vector is zero except the onelocated at the position of the physical source. When the velocity model is inaccurate, the other entries also contributeto decrease the data misfit, although their contribution is much less than the component located at the physical sourceposition. Surprisingly, the additional degree of freedom provided by this spatially-distributed source signature helps toestimate a more accurate source signature when the velocity model is inaccurate, compared to the case where sourceblending is not used. The proposed algorithm solves IR-WRI with unknown source signatures in an alternating mode,first, jointly for data-assimilated wavefields and source signatures via variable projection, and then it solves a quadraticoptimization for updating the model parameters. Because the extension created by the blended source assumption isartificial, we should correct its effects during the iteration. We propose two different algorithms for this correction.Numerical tests from the Marmousi II and 2004 BP salt model show that the proposed method is more robust (againstinaccuracies in velocity model and data) and faster than traditional methods for source signature estimation and veloc-ity model inversion.This paper is organized as follows. In the method section, we first show how the source signature estimation can becombined with the extended wavefield reconstruction subproblem of IR-WRI by variable projection when each sourceis processed separately. We show that the extended wave equation operator becomes source dependent. The secondpart of the method section reviews the source blending approach that is used to make the extended wave equation op-erator source independent and hence amenable to efficient multi-source processing with direct methods. Two slightlydifferent algorithms are proposed to implement the method. The paper continues with a numerical example section.We first assess the sensitivity of the source signature estimation to several parameters, such as the accuracy of thevelocity model, noise, and the distance between sources and receivers. Then, we present applications of IR-WRI withthe proposed efficient source signature estimation on the Marmousi II model and the BP salt model and compare theresults when the source signatures are known and when they are estimated without the efficient blending strategy.
Frequency-domain FWI for multi-source acquisition with unknown source signatures can be formulated as the follow-ing constrained optimization problem (Aghamiry et al., 2019b) minimize u i ,s i , m ∈M R ( m ) subject to (cid:26) A ( m ) u i = s i φ i , i = 1 , , ..., n s , Pu i = d i , i = 1 , , ..., n s (1)where m ∈ R N × is the model parameter vector (squared slowness), N is the number of discretized points of themedium, n s is the number of sources, A ( m ) = ∆ + ω Diag ( m ) ∈ C N × N is the Helmholtz operator, ω is the angularfrequency, ∆ is the Laplacian operator, Diag( x ) denotes a diagonal matrix with the entries of the vector x on itsdiagonal, u i ∈ C N × and d i ∈ C M × denote the wavefield and the recorded data for the i ’th source, respectively, P ∈ R M × N is the observation operator and M is the number of receivers. Also, s i ∈ C is the source signature forthe i ’th source at frequency ω , and φ i ∈ R N × is a sparse vector defining the i ’th source location. Finally, R ( m ) is an appropriate regularization function on the model domain and M is a convex set defined according to our prior3R-WRI with unknown sources, Aghamiry et al. A P
REPRINT knowledge of m . For example, if we know the lower and upper bounds on m then M = { m | m min ≤ m ≤ m max } . (2)IR-WRI solves the constrained problem (1) with the augmented Lagrangian method (or the method of multipliers).The augmented Lagrangian method combines a penalty term to relax the constraints during the early iterations anda Lagrangian term to control how accurately the constraint is satisfied at the convergence point (Nocedal & Wright,2006). In this method, the primal variable and the Lagrange multiplier or dual variable are updated in alternatingmode using a primal descent/dual ascent approach. Moreover, to make the computational cost tractable, we updatethe primal variables u and m in an alternating mode in the framework of the alternating-direction method of multi-pliers (ADMM) (Boyd et al., 2010). The reader is referred to Aghamiry et al. (2019b), Aghamiry et al. (2019a) andAghamiry et al. (2020a) for more details about the ADMM-based IR-WRI algorithm. In the last two references, R ( m ) implements a total-variation (TV) regularization and an hybrid TV+Tikhonov regularization, respectively. Comparedto the above references, we extend IR-WRI to update the source signatures jointly with the wavefields during thewavefield reconstruction subproblem through a variable projection.Beginning with an initial model m and assume ˆd i = ˆb i = , ADMM solves iteratively the multivariate optimizationproblem, equation (1), with alternating directions as (see Aghamiry et al., 2019b; Boyd et al., 2010, for more details) ( u k +1 i , s k +1 i ) = arg min u i ,s i Ψ( u i , s i , m k , ˆb ki , ˆd ki ) , i = 1 , , ..., n s m k +1 = arg min m ∈M n s (cid:88) i =1 Ψ( u k +1 i , m , s k +1 i , ˆb ki , ˆd ki ) ˆb k +1 i = ˆb ki + s k +1 i φ i − A ( m k +1 ) u k +1 i , i = 1 , , ..., n s ˆd k +1 i = ˆd ki + d i − Pu k +1 i , i = 1 , , ..., n s (3a)(3b)(3c)(3d)where Ψ( u i , s i , m , ˆb ki , ˆb ki ) = R ( m ) + (cid:107) Pu i − d i − ˆd ki (cid:107) + λ (cid:107) A ( m ) u i − s i φ i − ˆb ki (cid:107) , (4)is the scaled form of the augmented Lagrangian (Boyd et al., 2010, section 3.1.1), • k is the value of • at iteration k ,the scalar λ > is the penalty parameter assigned to the wave equation constraint, and ˆb ki ∈ C N × , ˆd ki ∈ C M × are the scaled Lagrange multipliers, which are updated through a dual ascent scheme by the running sum of the con-straint violations (source and data residuals) as shown by equations (3c)-(3d). The penalty parameter λ can be tunedin equation (4) such that the estimated wavefields approximately fit the observed data from the first iteration at the ex-pense of the accuracy with which the wave equation is satisfied, while the iterative update of the Lagrange multipliersprogressively corrects the errors introduced by these penalizations such that both of the observation equation and thewave equation are satisfied at the convergence point with acceptable accuracies.Here, we focus on the optimization subproblem (3a). The readers are referred to Aghamiry et al. (2019a, 2020a, 2021)for the closed-form expression of the optimization subproblem (3b) with bound constraints and different regulariza-tions and to Aghamiry et al. (2020b) for a more robust implementation of this subproblem against velocity modelerrors with phase retrieval.Optimization problem (3a) is quadratic in u i and s i and can be written as arg min u i ,s i (cid:107) Pu i − d i (cid:107) + λ (cid:107) A k u i − s i φ i (cid:107) , (5)where A k ≡ A ( m k ) . Fang et al. (2018) solve Eq. (5) jointly for u i and s i by gathering them in a single vector andsolve an ( N + 1) × ( N + 1) linear system (Their Eq. 24) instead of the N × N system of the original WRI (vanLeeuwen & Herrmann, 2013). Here, we want to solve Eq. (5) with the variable projection method (Golub & Pereyra,2003). By taking derivative of Eq. (5) with respect to s i , we get that s i satisfies φ T i ( A k u i − s i φ i ) = 0 , where • T denotes the conjugate transpose of • , and since φ T i φ i = 1 , we get s i = φ T i A k u i . (6)Substituting the above expression of s i into Eq. (5) leads to a mono-variate optimization problem for the wavefield as arg min u i (cid:107) Pu i − d i (cid:107) + λ (cid:107) Q i A k u i (cid:107) , (7)where Q i = I − φ i φ T i , and I is the identity matrix. In equation (7), the matrix φ i φ T i is a diagonal matrix with onenonzero coefficient equal to 1 at the location of the source i , and Q i is another diagonal matrix complementary to φ i φ T i : its diagonal entries equal to 1 except at the source position where the coefficient is zero. The second term in4R-WRI with unknown sources, Aghamiry et al. A P
REPRINT equation (7) penalizes the predicted source A k u i at all spatial points except at the physical source location consistentlywith the elimination (or projection) of s i , equation (6), from the optimization variables. Minimization of Eq. (7) withrespect to u i gives (note that Q T i Q i = Q i = Q i ) u k +1 i = (cid:0) P T P + λ A T k Q i A k (cid:1) − P T d i . (8)The explicit relation between the estimated source signature and the data can be obtained as s k +1 i = φ T i A k u k +1 i = φ T i A k (cid:0) P T P + λ A T k Q i A k (cid:1) − P T d i = φ T i (cid:0) G T k G k + λ Q i (cid:1) − G T k d i , (9)where G k = PA − k is the rank-deficient forward operator sampling the Green function A − k at receiver positions.Equation (9) shows that the source signature is estimated by first back propagating the data in time from the receiverpositions, i.e. G T k d i , and then correct the blurring effects induced by the limited bandwidth of the data and the limitedspread of the receivers by applying the inverse of the Hessian, i.e. (cid:0) G T k G k + λ Q i (cid:1) .The optimization problem, Eq. (7) and its closed-form solution, Eq. (8), share some similarities with the extendedsource reconstruction method described in Huang et al. (2018, their Eqs. 10 and 11) as a source-independent variantof WRI, although there are two differences in their formulation: first, their state variables are the extended sourcesinstead of the extended wavefields, i.e., b ei = A k u i where b ei is the i ’th extended source; Second, they used anotheranihilator function for Q i (Their equation 9), which is zero at the source location and linearly increases away from it.The difficulty with the method of Huang et al. (2018) is related to the presence of the Green’s functions A − k in theHessian (cid:0) G T k G k + λ Q i (cid:1) which makes the system very difficult, if not impossible, to solve exactly. The proposed method for joint estimation of source signature and data-assimilated wavefield, Eq. (7), is robust againstvelocity model error but it is computationally-expensive because the normal operator (cid:0) P T P + λ A T k Q i A k (cid:1) is source-dependent, hence preventing efficient multi-source processing with either direct or iterative solvers. This difficultyalso exists in the method of Fang et al. (2018) for joint estimation of wavefield and source signature (Their Eq. 24).Although they proposed a block matrix formulation to overcome the computational challenges (Their Eq. 33), it seemsinefficient and not applicable for large scale problems.To overcome this computational issue and design a fast and accurate multi-source signature and wavefield estimation,we assume that a virtual blended source generates each wavefield. The true signature of this blended source is thephysical source signature at the physical source location and is zero elsewhere. By doing so, each source can bewritten as Φs i where Φ = [ φ φ ... φ n s ] ∈ R N × n s is a tall matrix including the shifted delta functions ( φ i ) inits columns and the true s i ∈ C n s × is a vector that contains the i ’th physical source signature ( s i ) and it is zeroelsewhere. We stress at this stage that this reformulation of the source, Φs i , is equivalent to the original one s i Φ i given in equation (1).Plugging the new source expression in the objective function (5) and taking the derivative with respect to s i gives s i = Φ T A k u i . (10)Projecting this expression into the cost function and remembering that Φ T Φ = I give arg min u i (cid:107) Pu i − d i (cid:107) + λ (cid:107) QA k u i (cid:107) , (11)where Q = I − ΦΦ T = n s (cid:80) n s i =1 Q i .The closed-form expression of the wavefield u i at iteration k + 1 is now given by u k +1 i = (cid:0) P T P + λ A T k QA k (cid:1) − P T d i . (12)Comparing Eq. (12) and Eq. (8) shows that the Hessian (cid:0) P T P + λ A T k QA k (cid:1) is source independent, hence preservingthe benefit of direct solver method to process efficiently multiple sources once one LU factorization has been per-formed. This solves the computational issue. However, the new parametrization of the source makes the optimizationproblem blind to the fact that the vector s i ∈ C n s × should have only one non zero entry at index i . It gives equalprobability to all source positions to reconstruct the wavefield u i , equation (11), hence leading to a blended wave-field. This is highlighted by the fact that A k u i is potentially dense in the closed-form expression of the reconstructedsignature s i , equation (10), when the velocity model is inaccurate. Conversely, when the velocity model is the trueone, A true u i = s i Φ i , is sparse. The blended source assumption gives an extra degree of freedom to the optimizationproblem to decrease the cost function. Surprisingly, we will show in the Numerical results section that the estimatedsource signature obtained with the source blending approach, Eq. (11), is more accurate than the counterpart obtained5R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT without blending, Eq. (7), when we start IR-WRI from a rough initial velocity model.The source blending is artificial and its effects (extra non zero coefficients in s i vectors) must be removed duringiterations. In the following, we propose two algorithms to achieve this goal. For the sake of compact notations, werecast from now the optimization problem, eq. (1), in matrix form.Multi sources can be processed efficiently in frequency-domain modeling by gathering them in the right-hand side(rhs) of the Helmholtz system in a matrix format. Considering n s sources, the multi-rhs Helmholtz system is writtenas AU = ΦS where U = [ u u ... u n s ] ∈ C N × n s and S ∈ C n s × n s is a square diagonal matrix with the sourcesignatures on its main diagonal, S ii = s i . Introducing the data matrix D = [ d d ... d n s ] ∈ C M × n s gives thefollowing optimization problem for multi-source signature and wavefield estimations arg min U , S (cid:107) PU − D (cid:107) F + λ (cid:107) A k U − ΦS (cid:107) F , (13)where (cid:107) • (cid:107) F denotes the Frobenius norm. Solving Eq. (13) with a variable projection method gives S = Φ T A k U , (14)where the diagonal components of S are dominant, and the off-diagonal coefficients represent the non-physical sourcecomponents associated with each of the n s blended sources. Plugging the explicit expression of S into Eq. (13) leadsto a mono-variate optimization problem for U , the closed-form expression of which is given by U k +1 = (cid:0) P T P + λ A T k QA k (cid:1) − P T D . (15)Eq. (15) is the same as Eq. (12) but for multi-data D . At this stage, we didn’t impose any constraint on the structure of S which is potentially dense when the velocity model is inaccurate. Plugging the expression of U k +1 in equation (14)gives the explicit expression of SS k +1 = Φ T A k U k +1 = Φ T A k (cid:0) P T P + λ A T k QA k (cid:1) − P T D . (16)Before proceeding with the subsurface parameter updating, we extract the approximate signature of the physicalsources as s = diag ( S ) , where diag( X ) denotes a vector that contains the diagonal elements of matrix X . Thedifferent steps of IR-WRI with source signature estimation are reviewed in Algorithm (1), which begins with an initialmodel m and initial dual variables ˆB = and ˆD = .The U -subproblem in Algorithm 1 introduce errors in the extended wavefield reconstruction due to the source blend- Algorithm 1
IR-WRI with unknown sources
Require:
Initial velocity model m set ˆB = and ˆD = repeat U k +1 = arg min U (cid:107) PU − D − ˆD k (cid:107) F + λ (cid:107) QA k U − ˆB k (cid:107) F S k +1 = Diag ( diag ( Φ T A k U k +1 )) m k +1 = arg min m R ( m ) + λ (cid:107) A ( m ) U k +1 − S k +1 − ˆB k (cid:107) F ˆB k +1 = ˆB k + S k +1 − A k +1 U k +1 ˆD k +1 = ˆD k + D − PU k +1 until stopping conditions are satisfied.ing. These errors can be corrected iteratively by the action of the Lagrange multipliers, which are formed by the sourceresiduals computed without source blending (Line 6 of algorithm 1). This is shown by the fact that the diagonal com-ponents of the source signature matrix are used (Line 4 of algorithm 1) instead of the whole matrix, equation (16). Theright-hand side correction term ˆB k in the objective function of the U -subproblem in Algorithm 1 gathers the runningsum of the source residuals of previous iterations. This iterative refinement leads to the error forgetting property dis-cussed by Yin & Osher (2013) in the frame of Bregman iterations, which means that the error correction performed atthe current iteration is made independent of the error corrections performed at previous iterations. Here, this iterativesolution refinement by right-hand side updating is necessary to correct three sources of errors: the first results fromthe fact that each primal subproblem is solved keeping fixed the other primal variable, the second from the fact thatwe solve a constrained problem with a penalty method keeping the penalty parameter fixed and the third from thenon-physical source blending. Another application of iterative refinement in AVO inversion is presented in Gholamiet al. (2018), where the linearized Zoeppritz equations are used to simplify the primal problem, while the dual problem6R-WRI with unknown sources, Aghamiry et al. A P
REPRINT compensates the linearization-related errors by computing the residuals with the exact Zoeppritz equations.When we seek to reconstruct a complicated velocity model starting from a rough initial model, the Algorithm 1 wasn’table to fully remove the detrimental effects of the source blending and failed to reach the same minimizer as IR-WRIwith a known source. This prompts us to propose Algorithm 2 that includes one extra step compared to Algorithm 1 byre-estimating the wavefields (Line 5 of Algorithm 2) with the diagonal restriction of the source signature matrix (Line4 of Algorithm 2). By doing so, the pollution effects of the non-physical sources are removed from the reconstructedwavefields. The improvement provided by this wavefield refinement is illustrated in the next
Numerical results section.
Algorithm 2
IR-WRI with unknown sources with wavefield correction
Require: starting point m set ˆB = and ˆD = repeat U k + = arg min U (cid:107) PU − D − ˆD k (cid:107) F + λ (cid:107) QA k U − ˆB k (cid:107) F S k +1 = Diag ( diag ( Φ T A k U k + )) U k +1 = arg min U (cid:107) PU − D − ˆD k (cid:107) F + λ (cid:107) A k U − S k +1 − ˆB k (cid:107) F m k +1 = arg min m R ( m ) + λ (cid:107) A ( m ) U k +1 − S k +1 − ˆB k (cid:107) F ˆB k +1 = ˆB k + S k +1 − A k +1 U k +1 ˆD k +1 = ˆD k + D − PU k +1 until stopping conditions are satisfied.Algorithm 2 requires two LU decompositions at each IR-WRI iteration (one for U k + and the other one for U k +1 ),but still remains much cheaper than the algorithm performing one LU decomposition per source. We first investigate different aspects of the proposed method for efficient source signature estimation in IR-WRI(referred to as joint approach in the following) with the Marmousi II model (Fig. 1a) and compare its performancewhen each source is processed separately in IR-WRI (referred to as separate approach in the following), equation (3),and when the source signatures are estimated with the method of Pratt (1999) (referred to as conventional approach),as S = (cid:0) Φ T G T GΦ (cid:1) − Φ T G T D . (17)Then, we compare the performances of the proposed Algorithms 1 and 2 with those of separate approaches when thesource signatures are known (classical IR-WRI) and unknown, equation (3). We use the Marmousi II model and ascaled version of the left target of the challenging 2004 BP salt model (Billette & Brandsberg-Dahl, 2004) for thiscomparison.For all the numerical tests, we use a 9-point finite-difference staggered-grid stencil with PML boundary condition(along the model’s edges except for the top where the free-surface boundary condition is used) and anti-lumped massto solve the Helmholtz equation. We first illustrate the performance of the separate, joint and conventional methods for source signature estimation andinvestigate the robustness of these methods against the accuracy of the initial velocity model, noise in the recordeddata, the vertical distance between the source and the receiver profiles, and the number of sources. We use fourinitial models for these tests (Fig. 1b-1e), which are referred to as models 1 -4. The fixed-spread acquisition contains114 point sources, the source signatures of which are Ricker wavelets of different central frequency and initial phase(Figure 2), and a line of receivers spaced 50 m apart at the surface. The central frequency for each source signature isselected randomly between [7-15] Hz, and the peak of each wavelet is centered randomly between [0-0.4] s.
First, we put the line of sources at 75 m depth, generate the data with the true velocity model (Fig. 1a) and use the 1Dgradient velocity model (model 4, Fig. 1e) to estimate the source signatures with the conventional method (Eq. (17)),7R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT
Figure 1: Marmousi II model. (a) True model. (b-e) Smoothed versions of the true model, referred to as models 1 (b),2 (c), 3 (d) and 4 (f) in the paper.Figure 2: Amplitude spectrum of source signatures used in Marmousi II test.the separate method (Eq. (9)), and the joint method (Eq. (16)). We show the magnitude and phase of the estimatedsource signatures for a couple of sources with source numbers 18 and 33 (Fig. 2) in Fig. 3. First, the results clearlyshow the improvement achieved by the relaxed wave-equation methods (the separate and joint methods) (Fig.3b,c,e,f)compared to the conventional method (Fig.3a,d). Second, both separate and joint methods estimate accurate sourcesignatures but with a different computational burden (one LU decomposition for the joint method against 114 LUdecompositions for the separate method).To gain more quantitative insights on the accuracy of methods, we plot the relative error (RE) of the estimated sourcesignatures as a function of the source number for each method in Fig. 4. The RE for the estimated source signature is8R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT
Figure 3: Estimated source signature (magnitude and unwrapped phase) with different methods for a couple of sourceswith source numbers 18 and 33 in Fig. 2 when model 4 (Fig. 1e) is used as the initial velocity model. (a) Theconventional method of Pratt (1999), (b) the separate method, and (c) the joint method. The blue and red points showthe estimated source signatures for the first and second sources, respectively, and the black dashed lines show the trueones. (d-e) same as (a-c) but for unwrapped phase.defined as RE = (cid:113)(cid:80) f n j = f [ x ∗ j − ˆ x j ] (cid:113)(cid:80) f n j = f x ∗ j , (18)where x ∗ j and ˆ x j are the true and estimated source signatures at frequency j , respectively, and f and f n are theminimum and maximum frequencies, respectively. In this figure, we show the RE of the estimated source signaturesfor the separate and joint methods when the rough velocity model 4 (Fig. 1e) and the kinematically accurate model2 (Fig. 1c) are used as background velocity model. We don’t show the RE of the conventional method in this figurebecause it is much higher than those obtained with the separate and joint approaches. It is shown that when the initialvelocity model is rough, the joint method performs better than the separate method (the blue and red curves in Fig.4). This probably results from the extra degrees of freedom available in the joint method compared to the separatecounterpart. More precisely, the error in the estimated source signature generated by the inaccuracy of the initialvelocity model is entirely mapped at the physical location of the source in the separate method. In contrast, this erroris distributed across the different components of the blended source in the joint method, hence, providing a betterestimation of the source signature at the location of the physical source. On the other hand, as the velocity modelbecomes more accurate, the separate and joint methods reach the same accuracy for the source signature estimation(the green and orange curves in Fig. 4). We continue by assessing the robustness of the methods against the noise in the recorded data and the error in the initialvelocity model. We repeat the same test as before with different initial velocity models and different levels of randomGaussian noises in the data. The average RE over all the sources for the conventional, separate, and joint methods areshown in Figs. 5a-5c, respectively, as functions of the signal to noise ratio (SNR) of the data and the initial velocitymodel. In this paper, the SNR of data is defined asSNR = 20 log (cid:18)
Clean data RM S amplitudeN oise RM S amplitude (cid:19) . (19)First of all, the conventional method is robust against the noise in the data, but it is sensitive to the errors in the initialvelocity model, as illustrated in the previous section. In contrast, the separate and joint methods are robust against9R-WRI with unknown sources, Aghamiry et al. A P
REPRINT
Figure 4: RE of the estimated source signatures for separate and joint method when a rough initial model 4 (Fig. 1e)and the kinematically accurate initial model 2 (Fig. 1c) are used.Figure 5: Average RE over all the sources for the estimated source signatures of (a) conventional method, (b) separatemethod and (c) joint method as a function of SNR of recorded data and initial velocity model.inaccuracy of the initial velocity model due to the extended search space allowing for data fitting with an inaccuratemodel but are sensitive to the noise in recorded data due to the risk of noise overfitting. Although the RE of the separateand joint methods increases with the amount of noise, it remains however far less than the conventional method. It canbe seen that even in the worst scenario (with the lowest SNR data and a rough initial model), such methods based onwave-equation relaxation still can estimate an acceptable source signature.
In IR-WRI, the accuracy of the estimated source signature is directly controlled by the distance between the sourceand the closest receiver. This results because the extended wavefield, from which the source signature is estimated,equation (14), matches well the recorded data only near the receivers when the background velocity model is inaccu-rate. This implies that when the source is close to a receiver, it will be estimated from an accurate estimation of thewavefield at this receiver. Moreover, the impedance matrix in equation (14), which is built from a potentially inaccu-rate velocity model, will not generate significant errors when applied to the wavefield to generate the source when thelatter and the receiver are close to each other. In this section, we show the robustness of the different methods againstthis distance. We do the same test as before, but we change the depth of the source line while keeping the receiverline at the surface. The average RE of the estimated source signatures, summed over all the sources, are plotted inFig. 6 as a function of the depth of the source profile for the conventional, separate, and joint methods and for therough model 4 (Fig. 6a) and the kinematically accurate model 2 (Fig. 6b). First, it can be seen that for the verticaldistances up to 1500 m, the joint method has a better performance for both initial velocity models, but it becomesunstable beyond this vertical distance where the subsurface becomes more complex (green curves in Fig. 6). Thissuggests that the additional degrees of freedom in the joint method drives the least-squares problem, Eq. (13), towardan inaccurate local minimizer when the extended wavefield becomes too inaccurate at the source location. For surfaceacquisitions, this should however not be an issue in practice. For towed-streamer acquisitions, the sources are close tothe nearest receiver, and both of them are in the water. In seabed acquisition, the reciprocal sources are on the seabedand potentially far away from the receivers. However, the medium between the receiver and the source layouts (i.e.,10R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT
Figure 6: Average RE over all the sources for the estimated source signatures as a function of the vertical distancebetween the source and receiver lines when (a) the rough velocity model 4 (Fig. 1e), (b) the kinematically accuratevelocity model 2 (Fig. 1c) are used.Figure 7: Average RE over all the sources for the estimated source signatures as a function of the number of thesources when the rough initial velocity model 4 (Fig. 1e) and a line of 320 receivers at the surface are used.the water) is known. On land, areal acquisitions are classically designed with sources and receivers at the surface withshort nearest offset.
The other aspect that we need to investigate is the number of sources and receivers. We repeat the same test as beforeseveral times with the rough initial velocity model 4 and a line of 320 receivers with 50 m spacing at the surface.For each of them, we use a line of sources with a different number of sources ranging from 1 to 360 with an interval20 (sources) deployed at 75 m depth. The average RE over all the sources for the estimated source signatures as afunction of the number of sources for the conventional, separate, and joint methods is shown in Fig. 7. Both theconventional method and the separate method have a stable behavior (blue and red curves), while the joint methodbecomes unstable as soon as the number of sources exceeds the number of receivers. In this case, the source signatureestimation problem becomes under-determined in the case of the joint method, and equation Eq. (13) converges to theleast-norm solution. In practice, this issue can be easily bypassed by subdividing the sources into patches of suitabledimension, gathering possible closely-spaced sources.
We continue by assessing Algorithms 1 and 2 as well as the separate method for FWI on the Marmousi II modelwhen the inversion is started using the rough initial model 4 (Fig. 1e) and a 3 Hz frequency. The fixed-spread surfaceacquisition consists of 114 sources spaced 150 m apart with the source signatures depicted in Fig. 2 at the surface, and340 hydrophone receivers spaced 50 m apart at 75 m depth.We first show the effects of the source blending in the matrix S , Eq. (16), for the rough initial model 4 (Fig. 8a).The column of the matrix at X = 34 is plotted in black to give more precise insights on the relative magnitude ofthe diagonal and off-diagonal elements. Also, the off-diagonal elements of this matrix are plotted separately in Fig.11R-WRI with unknown sources, Aghamiry et al. A P
REPRINT
Figure 8: The estimated real components of source signature matrix, i.e. Φ T A k U k , (a-b) with the rough initialmodel 4(Fig. 1e) for 3 Hz, and (c-d) with the kinematically accurate velocity model 2(Fig. 1c) for 12 Hz. (a,c)shows Φ T A k U k where the column at X = 34 of these matrices is plotted in black to have a visual insight aboutthe magnitude of diagonal and off-diagonal elements. (b,d) show the non-diagonal components of (a,c) [ Φ T A k U k − Diag ( diag ( Φ T A k U k )) ].Figure 9: (a) The reconstructed wavefield using Eq. (15) using the rough initial model 4 (Fig. 1e) for 3 Hz. (b) Thedifference between (a) and the reconstructed wavefield with known true source signature. (c-d) Same as (a-b), butusing the kinematically accurate velocity model 2 (Fig. 1c) for 12 Hz.8b. It is shown that the diagonal elements of this matrix are dominant, the maximum amplitude of the off-diagonalcomponents being less than 1 percent of the maximum diagonal element. We remind that such off-diagonal elementspartially absorb the errors in the estimated physical source signatures when the initial velocity model is rough, but weneed to remove these effects during the inversion, which is the goal of Algorithms 1 and 2.The effects of the source blending are also seen in the reconstructed wavefields. The reconstructed monochromaticwavefield associated with the source located in Fig. 8a is shown in Fig. 9a. To assess the effects of the source blending,we show the difference between this wavefield and the extended wavefield reconstructed with the true source signaturein Fig. 9b. It can be seen that the differences are not significant. This is the worst scenario because of the significantinaccuracy of the velocity model 4. However, the estimated source signature matrix becomes close to a diagonal matrixas soon as the accuracy of the velocity model improves. This statement is verified in Figures 8c-8d and 9c-9d, whichare similar to Figures 8a-8b and 9a-9b, except that the source signature and the wavefield are now estimated from thekinematically accurate velocity model 2 and a frequency of 12 Hz. It can be seen that the source matrix, Fig. 8c, tendsto a diagonal matrix, and the differences between the estimated and the true wavefields tend to zero in Fig. 9d.We continue by performing the frequency-domain IR-WRI in the 3 Hz - 12 Hz frequency band with a frequencyinterval of 0.5 Hz. Mono-frequency batches are successively inverted following a classical frequency continuationstrategy. We perform three paths through the frequency batches to improve the inversion results, using one path’s finalmodel as the initial model of the next one. The starting and finishing frequencies of the paths are [3, 6], [3, 7], [3,12R-WRI with unknown sources, Aghamiry et al. A P
REPRINT
Figure 10: Marmousi II inversion results when the inversion starts from the crude initial velocity model 4 (Fig. 1e).Estimated velocity model using IR-WRI (a) with known sources, (b-d) with unknown sources using (b) separatemethod, (c) joint method with Algorithm 1, and (d) joint method with Algorithm 2.Figure 11: Marmousi II inversion results. (a-d) The model error (difference between true model and estimated model)for estimated models Figs. 10a-10d.12] Hz. We compare the results of IR-WRI with the known sources (Aghamiry et al., 2019b) (Fig. 10a), IR-WRI withunknown sources using the separate method (Fig. 10b), the joint method with Algorithm 1 (Fig. 10c), and finally, thejoint method with Algorithm 2 (Fig. 10d), to assess the effectiveness of the proposed methods. Also, the model error(the difference between the estimated and true velocity model) for the different estimated models (Figs. 10a-10d) areshown in Fig. 11a-11d. It can be seen that all of these methods work well but with a different computational cost.Also, IR-WRI with known sources (Fig. 10a) and IR-WRI using Algorithm 2 (Fig. 10d) are really close together andhave a better performance at the reservoir level (compare Figs. 11a and 11d).Finally, the wave-equation residual, data residual, and RE for IR-WRI with known and unknown sources are shownin Fig. 12. In summary, we see that the separate method and both of the joint-approach algorithms work well andcan reconstruct velocity models close to what we get from IR-WRI with known source signatures but with a differentcomputational burden. The test in the next section shows that this conclusion is not valid for a more complicatedvelocity model. 13R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT
Figure 12: Marmousi II test. A comparison between the convergence history of IR-WRI with known source (green),unknown sources with separate method (blue), joint method with Algorithm 1 (red) and joint method with Algorithm2 (orange) as a function of iteration number. (a) PDE misfit, (b) Data misfit, and (c) Model RE.Figure 13: BP Salt test. RE of the estimated source signatures for separate and the joint method when a rough initialmodel (Fig. 16b) and the accurate velocity model (Fig. 16f) are used.
We continue by assessing the performance of the methods on a re-scaled version of left-target of challenging 2004BP salt model (Billette & Brandsberg-Dahl, 2004) (Fig. 16a) when a 1D gradient initial model (Fig. 16b) and the3 Hz frequency are used to start the inversion. We use 65 point sources with 250 m spacing and a line of receiverswith 50 m spacing at 75 m depth. Like the previous test, the source signatures are random Ricker wavelets withdifferent central frequencies between [8-12] Hz and the initial phases between [0-0.4] s. We apply the inversion inthe 3 Hz-13 Hz frequency band with a frequency interval of 0.5 Hz. We perform three paths through the frequencybatches, using the final model of one path as the initial model of the next one, and each batch contains two frequencieswith one frequency overlap. The starting and finishing frequencies of the three paths are [3, 6], [4, 8.5], [6, 13] Hz,respectively. Bound-constrained Tikhonov + Total variation (BTT)-regularization (Aghamiry et al., 2020a) is appliedon IR-WRI for all the cases to decrease the ill-posedness of the problem. We first plot the RE for the estimated sourcesignatures of the separate and joint method as a function of the source number in Fig. 13 for the first (3 Hz, using initialvelocity model Fig. 16b) and the final iteration of the inversion (13 Hz, using updated velocity model 16f). We see thesame effects as those revealed by Marmousi II in Fig. 4 in the sense that the joint approach provides more accuratesource signature estimation. Also, the estimated source signature matrices, Eq. (16), at the first iteration and the finaliteration of IR-WRI, are shown in Fig. 14 and the related reconstructed monochromatic wavefields are shown in Fig.15. Again, we see similar effects as those revealed by the Marmousi II test, except that the off-diagonal elements arelarger relative to the diagonal counterparts, hence revealing the more complex structure of the BP salt model. Weshow in Fig. 16 the final results of BTT-regularized IR-WRI with the known sources (Fig. 16c), the unknown sourcesusing the separate method (Fig. 16d), the joint method with Algorithm 1 (Fig. 16e), and finally, the joint method withAlgorithm 2 (Fig. 16f). Also, the model errors for different estimated models are shown in Fig. 17. In contrast to theMarmousi II test, the different methods don’t converge to the same minimizer. Let’s consider the BTT-regularized IR-WRI with known sources as the benchmark model (Fig. 16c). Only IR-WRI with unknown sources using Algorithm14R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT
Figure 14: BP Salt test. The estimated real components of source signature matrix, i.e. Φ T A k U k , for the (a-b) firstand (c-d) final iteration of IR-WRI. (a,c) shows Φ T A k U k where the column at X = 34 of these matrices is plotted inblack to have a visual insight about the magnitude of diagonal and off-diagonal elements. (b,d) show the non-diagonalcomponents of (a,c) [ Φ T A k U k − Diag ( diag ( Φ T A k U k )) ].Figure 15: BP Salt test. (a) The reconstructed wavefield using Eq. (15) at the first iteration of IR-WRI using jointmethod with the rough initial model (Fig. 16b) for 3 Hz. (b) The difference between (a) and the reconstructed wavefieldwith true source. (c-d) Same as (a-b), but for 12 Hz at the final iteration of IR-WRI with updated velocity model Fig.16f.2 reaches approximately the same results. The failure of the separate method (Fig. 16d and 17b) probably results fromthe limited quality of the estimated wavelets at the early iterations when the initial velocity model is inaccurate (Fig.13). Also, the failure of Algorithm 1 (Fig. 16e and 17c) may result from the significant amplitudes of the off-diagonalelements of the estimated source signature matrix (Fig. 14), and it seems that the iterative refinement implemented inAlgorithm 1 is not enough to correct all of these effects. For such complicated velocity models, we need to recomputethe wavefields from the diagonalized source signature with Algorithm 2 (lines 4-5). We extended the recently proposed iteratively-refined wavefield reconstruction inversion (IR-WRI) to estimate theunknown source signatures. The source signatures and wavefields are jointly estimated with a variable projectionmethod during the extended wavefield reconstruction subproblem. We first show that the source signature estimationgenerates computational overhead when each source is processed separately because the extended wave equation15R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT
Figure 16: 2004 BP salt test: (a) True velocity model. (b) Initial velocity model. (c-f) BTT-regularized IR-WRI (c)with known sources, (d) with unknown sources using separate method, (e) with unknown sources using Algorithm 1,and (f) using Algorithm 2.Figure 17: 2004 BP salt test: (a-d) The model error (difference between true model and estimated model) for estimatedmodels Figs. 16c-16f.operator becomes source dependent. This computational overhead becomes prohibitive when the augmented waveequation system is solved with the direct method as a LU factorization needs to be performed for each source. Tobypass this issue and make the operator source independent, we proposed a method that blends the sources duringeach wavefield reconstruction. Accordingly, for each source of the experiment, the proposed method searches fora virtual blended source that best fits each single-source dataset during the wavefield reconstruction. Regardless of16R-WRI with unknown sources, Aghamiry et al.
A P
REPRINT the computational efficiency issue, we also show that when the background velocity model is inaccurate, this sourceblending provides a more accurate source signature estimation at the physical source location than the case wheresource blending is not used. This probably results from the additional degrees of freedom provided by the extravirtual sources. Once the spatially distributed source signatures have been estimated, we restrict them at the positionof the physical source to mimic the true source signatures, and we re-estimate the extended wavefields with theselocalized source signatures. Although we assume in this study that the source locations match the grid points of thecomputational domain, the method can be readily used for arbitrary source positions.
References
Abubakar, A., Hu, W., Habashy, T. M., & van den Berg, P. M., 2009. Application of the finite-difference contrast-source inversion algorithm to seismic full-waveform data,
Geophysics , (6), WCC47–WCC58.Aghamiry, H., Gholami, A., & Operto, S., 2019a. Implementing bound constraints and total-variation regularization inextended full waveform inversion with the alternating direction method of multiplier: application to large contrastmedia, Geophysical Journal International , (2), 855–872.Aghamiry, H., Gholami, A., & Operto, S., 2019b. Improving full-waveform inversion by wavefield reconstructionwith alternating direction method of multipliers, Geophysics , , R139–R162.Aghamiry, H., Gholami, A., & Operto, S., 2020a. Compound regularization of full-waveform inversion for imagingpiecewise media, IEEE Transactions on Geoscience and Remote Sensing , (2), 1192–1204.Aghamiry, H., Gholami, A., & Operto, S., 2020b. Robust wavefield inversion with phase retrieval, GeophysicalJournal International , (2), 1327–1340.Aghamiry, H., Gholami, A., & Operto, S., 2021. Full Waveform Inversion by Proximal Newton Methods usingAdaptive Regularization, Geophysical Journal International , (1), 169–180.Aravkin, A. & van Leeuwen, T., 2012. Estimating nuisance parameters in inverse problems, Inverse Problems , (11).Aravkin, A. Y., van Leeuwen, T., Calandra, H., & Herrmann, F. J., 2012. Source estimation for frequency-domainFWI with robust penalties, in Expanded Abstracts, 74 th Annual EAGE Meeting .Billette, F. J. & Brandsberg-Dahl, S., 2004. The 2004 BP velocity benchmark, in
Extended Abstracts, 67 th AnnualEAGE Conference & Exhibition, Madrid, Spain , p. B035.Boyd, S., Parikh, N., Chu, E., Peleato, B., & Eckstein, J., 2010. Distributed optimization and statistical learning viathe alternating direction of multipliers,
Foundations and trends in machine learning , (1), 1–122.Fang, Z., Wang, R., & Herrmann, F. J., 2018. Source estimation for wavefield-reconstruction inversion, Geophysics , (4), R345–R359.Gholami, A., Aghamiry, H., & Abbasi, M., 2018. Constrained nonlinear AVO inversion using Zoeppritz equations, Geophysics , , R245–R255.Golub, G. & Pereyra, V., 2003. Separable nonlinear least squares: the variable projection method and its applications, Inverse problems , (2), R1.Haber, E., Ascher, U. M., & Oldenburg, D., 2000. On optimization techniques for solving nonlinear inverse problems, Inverse problems , (5), 1263.Huang, G., Nammour, R., & Symes, W. W., 2018. Source-independent extended waveform inversion based on space-time source extension: Frequency-domain implementation, Geophysics , (5), R449–R461.Li, M., Rickett, J., & Abubakar, A., 2013. Application of the variable projection scheme for frequency-domain full-waveform inversion, Geophysics , , R249–R257.Marfurt, K., 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations, Geophysics , , 533–549. 17R-WRI with unknown sources, Aghamiry et al. A P
REPRINT
Nocedal, J. & Wright, S. J., 2006.
Numerical Optimization , Springer, 2nd edn.Plessix, R. E. & Cao, Q., 2011. A parametrization study for surface seismic full waveform inversion in an acousticvertical transversely isotropic medium,
Geophysical Journal International , , 539–556.Pratt, R. G., 1999. Seismic waveform inversion in the frequency domain, part I: theory and verification in a physicalscale model, Geophysics , , 888–901.Pratt, R. G., Shin, C., & Hicks, G. J., 1998. Gauss-Newton and full Newton methods in frequency-space seismicwaveform inversion, Geophysical Journal International , , 341–362.Rickett, J., 2013. The variable projection method for waveform inversion with an unknown source function, Geophys-ical Prospecting , (4), 874–881.Song, Z., Williamson, P. R., & Pratt, R. G., 1995. Frequency-domain acoustic-wave modeling and inversion ofcrosshole data. Part 2 : Inversion method, synthetic experiments and real-data results, Geophysics , (3), 786–809.Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation, Geophysics , (8), 1259–1266.van Leeuwen, T. & Herrmann, F., 2016. A penalty method for PDE-constrained optimization in inverse problems, Inverse Problems , , 1–26.van Leeuwen, T. & Herrmann, F. J., 2013. Mitigating local minima in full-waveform inversion by expanding thesearch space, Geophysical Journal International , , 661–667.Virieux, J. & Operto, S., 2009. An overview of full waveform inversion in exploration geophysics, Geophysics , (6),WCC1–WCC26.Yin, W. & Osher, S., 2013. Error forgetting of Bregman iteration, Journal of Scientific Computing ,54