{'fraction_non_alphanumeric': 0.031042228212039533, 'fraction_numerical': 0.009838274932614556, 'mean_word_length': 4.6962640736949846, 'pattern_counts': {'":': 1, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Full Waveform Inversion by Source Extension: Why it works

An extremely simple single-trace transmission example shows how an extended source formulation of full waveform inversion can produce an optimization problem without spurious local minima ("cycle skipping"). The data consist of a single trace recorded at a given distance from a point source. The velocity or slowness is presumed homogeneous, and the target source wavelet is presumed quasi-impulsive or focused at zero time lag. The source is extended by permitting energy to spread in time, and the spread is controlled by adding a weighted mean square of the extended source wavelet to the data misfit, to produce the extended inversion objective. The objective function and its gradient can be computed explicitly, and it is easily seen that all local minimizers must be within a wavelength of the correct slowness. The derivation shows several important features of all similar extended source algorithms. For example, nested optimization, with the source estimation in the inner optimization (variable projection method), is essential. The choice of the weight operator, controlling the extended source degrees of freedom, is critical: the choice presented here is a differential operator, and that property is crucial for production of an objective immune from cycle-skipping.


INTRODUCTION
Full Waveform Inversion (FWI), or estimation of earth structure by model-driven least squares data fitting, is now well-established as a useful tool for probing the earth's subsurface (Virieux and Operto, 2009;Fichtner, 2010). However, so-called "cycle-skipping", the tendency of iterative FWI algorithms to stagnate at suboptimal and geologically uninformative earth models, still impedes its use. Because the computational size of field inversion tasks is very large, only iterative local (descent) minimization of the data misfit function is computationally feasible. However local descent methods avoid suboptimal stagnation only if initial models are already quite close to optimal, in the sense of predicting the arrival times of seismic events to within a small multiple of a dominant wavelength (Gauthier et al., 1986;Plessix et al., 2010). This paper concerns one of the several ideas that have been advanced to overcome cycle-skipping, namely so-called extended inversion (Symes, 2008). "Extended" signifies that addional degrees of freedom are provided to the modeling process, in the hope of opening up more effective routes to geologically informative models with acceptable data fit. Since these extended degrees of freedom are not part of the basic physics chosen to model the data acquisition process, they should be suppressed in the eventual solution. Extended inversion methods differ by the choice of additional degrees of freedom, and by choice of penalty applied to eliminate them in the final result.
Many of these extended inversion concepts sound plausible, and appear to work at least to some extent as one might hope from their heuristic justifications. However very few of these approaches have been underwritten by mathematical argument: in essence, they are mostly justified only "in the rear-view mirror", with no assurance that failure is not just around the corner, at the next example. On top of that, some of these approaches, for example those based on the computationally attractive Variable Projection Method ("VPM") of Golub and Pereyra (2003), are cast in such form that the reasons for success are not readily apparent.
This note shows exactly how VPM leads to successful velocity updates for a particular extended inversion approach to a very simple inverse problem, which asks that a homogeneous velocity field be deduced from one trace at known offset. I put forward this inverse problem and extension-based solution not because there are not simpler ways of answering the question it poses -there certainly are -but because the formal ingredients of waveformbased velocity estimation in this very simple setting are common to many similar extended inversion algorithms, and because in this case every computation can be done analytically, nearly to completion. In particular, it becomes clear why the VPM gradient formula produces a constructive update, with no possibility of stagnation away from the global minimum.
The extended inversion approach developed here uses a source extension, in which source parameters form the additional degrees of freedom. This type of extension presumes that the actual or target source is constrained in some way; the extended source is allowed to violate the constraint. For recent overview of source extension methods, see Huang et al. (2019). Source extension methods have computational complexity approximately the same as that of FWI, a signal advantage over the alterantive medium extension class described for instance in (Symes, 2008).
For the problem considered here, the source model amounts to a wavelet, and the target wavelet is assumed to be non-zero only in a short time interval (an approximate impulse, perhaps as the result of signature deconvolution). The extension consists in permitting energy to spread in time at intermediate iterations of the inversion. A simple penalty for energy spread (second moment of square amplitude) drives the extended source towards a focused source approximately satisfying the assumed constraint. Not all penalties are created equal: the penalty used here has the critical (pseudo-)differential attribute necessary for avoidance of cycle-skipping, as will be explained in the Discussion section.
I begin with a quick sketch of constant density acoustics, and describe the single-trace transmission inverse problem. To make the role of data frequency content clear, I introduce a family of noise-free data parametrized by wavelength. For completeness, I show how the standard FWI approach to this problem generates multiple local minima that will be found by any descent method unless the initial estimate predicts travel time from source to receiver with an error on the order of a wavelength. The next section describes the source extension objective, and the reduced objective produced by VPM. As VPM eliminates the extended source, this function depends only on the velocity, just as does the FWI objective. A nearly-explicit calculation of the VPM gradient shows that the only stationary points are "within a wavelength" of the correct velocity, used to build the data: that is, cycle-skipping cannot occur. The paper ends with a discussion of the parallels between the calculations presented here and the structure of other extended inversion methods applicable to fieldscale velocity estimation, and the critical role that the differential nature of the extension penalty plays in the success of this and other extension methods.

PRELIMINARIES
Assume small amplitude (linearized) acoustic propagation, constant density, and isotropic point source and receiver. Denote by m(x) the slowness (reciprocal velocity) at spatial position x, f (t) the time dependence of the point source ("wavelet") at location x = x s . Then the (excess) pressure field p(x, t) obeys a scalar wave equation: Suppose that a single trace is recorded, at distance r > 0 from the source position x s . The dominant information in a single trace is the transient signal time of arrival, constraining only the mean slowness in the region between source and receiver, so assume that the slowness is constant, that is, independent of position x. The pressure field is simply the the source wavelet f (t) convolved with the acoustic Green's function, for which an analytic expression is available in the constant m case (Courant and Hilbert, 1962): The receiver location x r lies at distance r from the source location x s , that is, |x r −x s | = r. The predicted signal at p(x r , t) depends nonlinearly on the slowness m and linearly on the source wavelet f . Therefore it is naturally represented as the action of a m-dependent linear operator S[m] on f : Ignoring amplitude, this map implements a m-dependent time shift. This time shift operator is the basis of many descriptions of the cycle-skipping phenomenon (for example, Virieux and Operto (2009) Figure 7), so it is unsurprising that an analysis of cycle-skipping can be based on the simple modeling operator described above, which amounts essentially to a time shift. To make the link with wavelet frequency content manifest, I introduce a family {f λ } of wavelets indexed by λ, a parameter having dimensions of time, The argument s of the "mother wavelet" f 1 is nondimensional. The only constraints placed on f 1 are that (i) f 1 (s) = 0 for |s| ≥ 1, and (ii) f 1 has positive mean-square, that is, does not vanish identically. Note that the scaling is such that the mean-square is independent of λ.
I shall refer to λ as "wavelength": if f 1 has a dominant period of oscillation, then so does f λ , and it is proportional to λ.
To this family of wavelets and a choice of target slowness m * corresponds a family of noise-free data This family of data in turn defines a family of inverse problems, to which I now turn.

FULL WAVEFORM INVERSION
The preceding section provided all of the raw ingredients to define full waveform inversion for estimation of m from a single trace. It is only m that is to be determined: the λ−dependent family of wavelets {f λ } is regarded as known, along with the data family {d λ }. The aim is to chose m to minimize for all values of λ > 0.
Written out in detail, this objective function is using the previously observed independence of f λ from λ.
That is, for slowness m in error by more than 2λ/r from the target slowness m * , the FWI objective J FWI is perfectly flat: all nearby values of m are local minima. Therefore the local exploration of the objective gives no useful information whatever about constructive search directions towards the global minimizer m = m * , for which the objective value is of course = 0, if the initial estimate of m in error by an amount greater than a fixed multiple of λ. This is precisely the behaviour of FWI noted many times in the literature: the model must be known to "within a wavelength".

EXTENDED SOURCE INVERSION
As mentioned in the introduction, the modeling operator introduced in the last section, m → S[m]f λ , may be extended simply by including the source wavelet as one of the model parameters: that is, the model vector becomes (m, f ), and the modeling operator, (m, f ) → S[m]f . The reader will have no trouble seeing that the data misfit using this extended modeling operator can always be made to vanish entirely by proper choice of wavelet f , unless f must satisfy some additional constraints. Huang et al. (2019) describe a plethora of possible constraints for this and similar source extensions. Many (but not all) take the form of a quadratic penalty, that is, the mean square of Af , A being a suitable operator, commonly dubbed an annihilator: in many examples the ideal output for A applied to a source obeying the target constraints is the zero vector (Symes, 2008). With an annihilator A, to be chosen below, the penalty form of extended inversion is: minimize over The choice of penalty weight α has a profound influence on the character of this optimization problem. For the moment, I will mandate only that α > 0. A detailed discussion of methods for setting α is beyond the scope of this brief paper. In the Discussion section, I mention an effective method for choosing α to optimize convergence of iterative minimization, that applies (and has been applied) to much larger-scale inversion problems with features similar to those of this simple example.
While minimization of J ESI might be tackled directly -by alternately minimizations between m and f , or by computing updates for m and f simultaneously -such joint mimization performs poorly, as Huang (2016) has shown. The reason for this poor performance is that J ESI has dramatically different sensitivity to m versus f , especially for high frequency f , as the reader will see below. Instead, a nested approach, in which f is eliminated in an inner optimization, generally gives far better numerical performance. This Variable Projection Method (VPM) (Golub and Pereyra, 2003) takes advantage of J ESI being quadratic in f to solve for f given m, thus producing a reduced objective of m alone: where f [m] is the minimizer of J ESI over f for given m. For the problem considered here, J VPM is explicitly computable. First observe that apart from amplitude, S[m] is unitary: Therefore the normal equation for the minimizer on the RHS of equation 10 is At this point I have to come clean about the actual choice of A. Recall that A is called an "annihilator", as it vanishes when applied to the target source field, at least in an idealized limit. Thus the choice of A depends on modeling assumptions, not fundamental physics, as is characteristic of extended source inversion. As λ → 0, the target source wavelet f λ focuses at time t = 0, in the sense that it vanishes for |t| > λ. Therefore I choose A to penalize energy away from t = 0: This particular annihilator has been employed in earlier papers on extended source inversion (Plessix et al., 2000;Luo and Sava, 2011;Warner and Guasch, 2014;Huang and Symes, 2015;Warner and Guasch, 2016;Huang et al., 2017).

Symes
With this choice of A, the normal operator on the LHS of 13 is simply multiplication by a positive function of time, and can be inverted by inspection. Using the identity 11 for the adjoint operator, obtain thanks to the definition of the data d λ (equation 5). After a little algebra, and So The gradient of J VPM can be extracted by elementary means from the identity 18, but instead we will derive it in using some important features that this inverse problem shares with other inverse wave problems with more complex physics. To begin with, the gradient of a VPM objective of the form 10 is given by the formula This easily derived result is in some sense the main content of (Golub and Pereyra, 2003 and is NOT the same as the adjoint denoted by S[m] T , which is the source-space-to-dataspace adjoint. That is, this adjoint make this relation true, for all m, δm, f, and d: On the left side is the dot product in model space -since model space is just 1D in this example, that's just the numerical product. On the right is the dot product in data space, in the idealized continuum limit.
Note that the VPM gradient formula is remarkable in two ways. First, it is exactly the same as the FWI gradient, that is, the gradient of J FWI if you happen to insert the solution f The key to unlocking the meaning of the VPM gradient formula for this and similar problems is a remarkable feature of the derivative operator DS[m]: from the definition 3, where That is, Q[m]δm is a skew-adjoint operator depending linearly on δm -more on this below.
A calculation, detailed in Appendix A, yields an expression for the gradient 19 in terms of Q: Here Note that the annihilator A is explicitly present in 23. The structure displayed in this expression is common to many other extended inversion methods. The extreme simplicity of the factorization 21, 22 and the gradient expression 23 is modified for more complex inversion problems by asymptotically negligible corrections (Symes, 2014;ten Kroode, 2014;Symes, 2015).
Remember that A T A amounts to multiplying by t 2 , and Q is the scaled time derivative Combine this identity with the formula 15 for the solution of the inner problem, and divide out the common factor δm, to obtain Therefore the integral on the RHS of equation 26 can be re-written If m > m * + λ/r, then the entire interval of integration is a proper subset of the negative half-axis. Consequently the first factor in the integrand satisfies t 1 + (4πr) 2 α 2 t 2 ≤ −(m − m * )r + λ < 0 over the interval of integration. Consequently, the integral is negative (since f λ has a positive mean square) and so the gradient is positive. Similar reasoning applies to the case m < m * − λ/r.

Symes
• if m > m * + λ/r, then ∇J VPM [m] > 0, and That is, J VPM has no local minima further than O(λ) from the global minimum: the gradient has the correct sign and slowness updates computed from it will be constructive, unless the slowness estimate is already "within a wavelength" of being correct.
On the other hand, careful examination of equation 26 shows that J VPM is not convex: there is an inflection point O(1/α) from the global minimizer m * , and in fact as α → ∞ for fixed λ, the J VPM approximates J FWI .

DISCUSSION
There remain several important points to be made.
First, the formal computations centering on the operator Q (equations 21 through 23) depend only on the relation 21 and the skew-symmetry of Q, which hold with minor modifications for other more complex waveform inversion problems. These other more complex problems do not submit to such a simple treatment as is shown in the equations following 23, but for example it is possible in some cases to use a relation analogous to 21 to extract a relation between extended waveform inversion and traveltime tomography, via analysis of the Hessian at a zero-residual global minimizer. See for instance ten Kroode (2014); Symes (2014Symes ( , 2015. Up to that point the reasoning is quite general, and central to the understanding gained so far of extended inversion methods. Another issue mentioned in the last section is that the penalty weight α must be chosen somehow, and its choice influences heavily the convergence rate of iterative optimization algorithms applied to J VPM (or J ESI , though (as also mentioned above) it is a poor candidate for local optimization). Two solutions to the weight assignment problem has emerged in the last few years, both quite effective. One applies directly to VPM problem, via use of the Discrepancy Principle (Morozov, 1984). This algorithm controls the size of the data misfit term in J ESI to lie in an interval chosen to contain the assumed level of data noise, by adjusting α sporadically as the optimization proceeds. While this approach requires an estimate of data misfit to be achieved at the optimal model, it is very effective in maintaining convergence (Fu and Symes, 2017). Another approach modifies the VPM problem by means of the Augmented Lagrangian method (Nocedal and Wright, 1999). This reformulation appears to suppress much of the sensitivity to α of VPM optimization, and has been successfully used in extended inversion (Aghamiry et al., 2019).
Perhaps the most important general message implicit in the example presented here is that the choice of the annihilator A determines whether the extended inversion algorithm achieves global or semiglobal convexity, as is accomplished in the present example. The annihilator used here has a property whose importance can be guessed by examining the gradient formula 23. The operator Q is a first order differential operator. The gradient is a quadratic form whose Hessian is the communator [Q, A T A], and whose argument is f [m]. In order that the VPM objective be continuous for any model m and finite energy data d, this form should admit any finite-energy source as argument: in technical terms, it should be a bounded (or continuous) operator on the Hilbert space of finite energy traces.
If one asks for a bit more, namely that VPM objective function have derivatives of arbitrary order, then it is not too hard to see that the iterated commutators [Q, ...[Q, A T A]...] must all be bounded operators. This is a very strong restriction on A T A: this operator must be pseudodifferential, that is, a combination of differential operators and powers of the Laplace operator (Taylor, 1981). For the present operator, A T A is multiplication by t 2 , a very simple pseudodifferential operator. For more discussion of this requirement in the context of annihilators in extended inverseion, see Symes (2008), where you can also find references to the technical backstory.
This constraint actually rules out some popular approaches to FWI. To begin with, basic least-squares FWI as formulated in the third section above can be reformulated as a quadratic form whose Hessian turns out not to be pseudodifferential. Therefore the FWI objective is not smooth in data and model jointly, a fact that is linked to the cycle-skipping behaviour demonstrated above. More surprising, perhaps, the same turns out to be true for Wavefield Reconstruction Inversion (WRI), an extended source inversion algorithm introduced by van Leeuwen and Herrmann (2013), and further developed by van Leeuwen and Herrmann (2016), Wang et al. (2016), and Aghamiry et al. (2019), amongst others. This approach turns out to be closely linked to basic FWI, and can be formulated as minimization of a similar quadratic form which once again has a nonpseudodifferential Hessian. Not coincidentally, it also exhibits cycle-skipping behaviour. Just as for FWI, in simple cases such as the problem studied in this paper, WRI can be shown explicitly to have local minima far from the global minimum, and a region of attraction for the global minimum on the order of a wavelength in diameter, just as does FWI (Symes, 2020).

CONCLUSION
Desipite its simplicity, the single-trace transmission inversion problem proves typical of many more complex waveform inversion problems. The structure of the derivative is similar in many of these problems, and for the particularly simple one explained here, can be analysed on paper to the point of showing explicitly why a simple extended source approach to waveform inversion works -that is, generates an objective all of whose local minima are "within a wavelength" of the global minimizer. Otherwise stated, this particular extended source inversion is genuinely immune to cycle-skipping. The simple structure of this problem showcases the importance of the variable projection reduction (elimination of the extended source) and a proper choice of annihilator in the formulation of the basic objective. It also makes clear the central role played by a factorization of the linearized modeling operator, a feature shared with many more complex extended source methods applicable at field scale.