Acoustic- and elastic-waveform inversion with total generalized p-variation regularization

Geophysical models usually contain both sharp interfaces and smooth variations, and it is difficult to accurately account for both of these two types of medium parameter variations using conventional full-waveform inversion methods. In addition, sparse geometry, noisy data and source encoding usually lead to strong inversion artifacts. We develop a novel full-waveform inversion method for acoustic and elastic waves using a total generalized p-variation regularization scheme to address these challenging problems. We decompose the full-waveform inversion into two subproblems and solve these two minimization subproblems using an alternating-direction minimization strategy. One important advantage of the total generalized p-variation regularization scheme is that it can simultaneously reconstruct sharp interfaces and smooth background variations of geophysical parameters. Such capability can also effectively suppress the noises in source-encoded inversion and sparse-data inversion, or inversion of noisy data. We demonstrate the advantages of our new full-waveform inversion algorithm using a checkerboard model, a modified elastic SEG/EAGE overthrust model, and a land field seismic dataset. Our results of synthetic and field seismic data demonstrate that our new method reconstructs both smooth background variations and sharp interfaces of subsurface geophysical properties accurately, reduces the inversion artifacts caused by source encoding, noisy data or insufficient data coverage effectively, and provides a useful tool for accurate and reliable inversion of field seismic data.


Introduction
Full-waveform inversion (FWI) attempts to reconstruct subsurface medium properties by iteratively minimizing the difference between synthetic and observed data (Tarantola, 1984(Tarantola, , 1986Mora, 1987Mora, , 1988. As self-explained by its name, FWI uses full wavefield information, including both the amplitude and the traveltime of seismic signals, to invert for subsurface medium properties. Theoretically, FWI should be the most accurate inversion method for subsurface model building. However, numerous studies have shown that FWI is a highly nonlinear, ill-posed inverse problem (e.g., Luo and Schuster, 1991;Virieux and Operto, 2009). The misfit function of FWI may have numerous local minima. This characteristics of FWI usually leads to unsuccessful inversions and unreliable results when the initial model used in FWI is far away from the true model. This problem is often true in practical applications, particularly when the subsurface geology is complicated and/or the acquired seismic data are sparse and noisy.
Given the difficulties in applying the conventional FWI to practical problems, numerous studies have been conducted to seek a practically applicable and reliable inversion scheme, and continuous efforts is one of the most active research areas in the geophysical community. These studies are mostly in three categories.
The first category is employing a more convex function as the misfit function of FWI, rather than using the traditional 2 -norm waveform difference. The misfit function in the conventional FWI relies on the absolute waveform matching between synthetic data and observed data, and cycle skipping issue can easily occur if there is a significant phase and/or amplitude difference. Studies of seeking more convex function include the wave-equation traveltime based misfit function (Luo and Schuster, 1991;Luo et al., 2016), the envelope-based misfit function (Wu et al., 2014;Chi et al., 2014), the instantaneous-phase-based misfit function (Bozdag et al., 2011;Jiao et al., 2015), the correlation-based misfit function (Van Leeuwen and Mulder, 2010;Luo and Sava, 2011;Chi et al., 2015;Choi and Alkhalifah, 2016), the deconvolution-based misfit function (Luo and Sava, 2011;Warner and Guasch, 2016), the dynamic-time-warping-based misfit function (Ma and Hale, 2013), the Huber-norm misfit function (Guitton and Symes, 2003;Ha et al., 2009), and the optimal transport approach (Métivier et al., 2016a,b;Yang et al., 2017Yang et al., , 2018, etc. These misfit functions are generally more convex with respect to model perturbations, and therefore less prone to the cycle skipping issue. The second category is conducting FWI in a Gauss-Newton or quasi-Newton minimization framework. Conventional FWIs are essentially based on the first-order perturbation theory. The second-order terms (or equivalently the Hessian) can be important in FWI to obtain accurate and high-resolution model estimations (Pratt et al., 1998;Tang and Lee, 2010;Fichtner and Trampert, 2011). The Hessian can also be important in multi-parameter inversion (Pan et al., 2016). A direct computation of the inverse Hessian can be prohibitively expensive using current computational architecture (Fichtner and Trampert, 2011), and therefore several approximation methods were developed, including the limited-memory BFGS (L-BFGS) scheme (e.g., Nocedal and Wright, 2006), the truncated Newton method (Métivier et al., 2013(Métivier et al., , 2017, the quasi-Newton method with projected Hessian (Ma and Hale, 2013), combined Newton and conjugate gradient schemes (Epanomeritakis et al., 2008), and pseudo-Gauss-Newton scheme (Pan et al., 2015), etc. Several preconditioners were also developed to approximated the Hessian (e.g., Zhang et al., 2012).
The third category is employing regularization to accelerate the convergence. Applying regularization to geophysical inverse problems has a fairly long history (Zhdanov, 2002), and the Tikhonov regularization is probably the most frequently used regularization scheme (Tikhonov et al., 1995;Asnaashari et al., 2013).
The Tikhonov regularization tends to produce smooth models. In image analysis and processing, the total variation (TV) method (Rudin et al., 1992) was developed to promote sharp interfaces in the image. The method has been applied in geophysical inverse problems to promote the sharp interfaces of layers (Anagaw, 2011). Guitton (2012) designed a blocky regularization scheme for FWI. He employed a 1 -norm and a Cauchy function to enforce blockiness of the model. The blocky regularization is essentially edge-promoting: it can enforce sharp interfaces even at locations where interfaces probably do not exist, a feature in accordance with minimizing the total variation of a model. Therefore, the inversion results tend to be piecewise constant. This feature is also an inherent limitation of the first-order TV regularization.
The first-order TV regularization can be less efficient when the data is noisy or the model is complex. Lin and Huang (2014) developed a novel regularized FWI based on a modified TV (MTV) regularization scheme. They decomposed FWI into two interlacing inversion problems: one is the conventional FWI with a Tikhonov regularization term, and the other is an image denoising problem using the first-order TV method. The second problem is solved efficiently with the split-Bregman iteration, a technique that is proven to be especially suitable to produce clean and accurate TV denoising results (Goldstein and Osher, 2009). The model parameter and the auxiliary model parameter are updated in an alternating-direction approach, resulting in an efficient first-order TV regularized FWI. However, this regularization scheme is based on the first-order TV.
Although it avoids spike noises and provides more reliable inversion results, it still cannot avoid the inherent disadvantage of the first-order TV, i.e., the staircase artifacts in inversion results. Esser et al. (2016) developed an asymmetric TV regularized FWI. The asymmetric TV, or hinge-loss constraint TV as they called, penalizes the model discontinuities only in the vertical direction. They show that this asymmetric TV regularized FWI can facilitate the automatic delineation of high-contrast medium parameter anomalies such as salt bodies, as well as depth structures.
The Tikhonov and TV regularizations are not the only regularization schemes used in fullwaveform inversion. Guitton et al. (2012) developed a preconditioned FWI using the geological information derived from the migration image. The preconditioner is solved via estimating the local dip information from the migration image followed by a directional Laplacian filter (Hale, 2007). Lewis et al. (2014) developed a similar approach based on the anisotropic diffusion filtering, where the structural tensor is estimated from a migration image. In their approach, the gradient is preconditioned based on the geological structures derived from structural images, thus becomes geologically meaningful. In the case of sparse sources and therefore insufficient subsurface model coverage, such preconditioning can facilitate faster convergence towards meaningful inversion results. Xue et al. (2017) developed a similar regularization scheme using sparsity promotion in the seislet domain.
The total generalized variation (Bredies et al., 2010;Knoll et al., 2011;Zhang et al., 2016) is a technique that incorporates the higher-order total variations in image reconstruction. It is usually applied to various problems in its second-order form with a 1 -norm framework. The total generalized variation can reconstruct both the sharp interfaces and smooth variations of an image, leading to fewer artifacts compared with the first-order total-variation approach.
We develop a novel full-waveform inversion method for acoustic and elastic waves using a total generalized p-variation regularization scheme (TGPV-FWI). We combine the advantages of pnorm compressive sensing techniques and the second-order total variation to reconstruct both sharp interfaces and smooth background variations of geophysical parameters. We formulate our TGPV-FWI in an alternating-direction minimization framework. We decompose the TGPV-FWI into two interlacing minimization problems. The first minimization problem is a conventional FWI problem with the Tikhonov regularization, whereas the second minimization problem is a second-order compressive sensing model denoising problem. Once the second subproblem is solved correctly, we only need to add the difference between the FWI model in the last inversion and its "denoising" result from the second minimization problem to the gradient in the current inversion step. This strategy forces the inversion to evolve towards the TGPV model that is less prone to artifacts when seismic data are sparse, noisy or when using dynamic source encoding. With this alternating minimization framework, the FWI eventually converges to a more accurate and reliable result.
Our paper is organized as follows. In the Methodology section, we present the formulation of our new TGPV-FWI. We present the detailed algorithms for solving the second subproblem and our new TGPV-FWI in two appendices. In the Numerical Results section, we give three numerical examples, including two synthetic data examples and one field data example, to verify the advantages of our new TGPV-FWI over FWI with the conventional regularization schemes including the Tikhonov and TV regularizations. We give our findings in the Conclusions section.

Methodology
We formulate our new full-waveform inversion with the total generalized p-variation regularization (TGPV-FWI) as the following minimization problem: where d is an observed dataset, f (m) is a synthetic dataset, m is the medium parameter model to be inverted, and T p (m) is the TGPV regularization term.
Unlike the conventional regularization schemes such as the Tikhonov and TV regularizations that are defined using the 2 -norm and 1 -norm, respectively, the TGPV regularization term T p (m) is defined through a p -norm minimization problem: where w = (w x , w y ) is an auxiliary vector variable. The p -norm (or more precisely, the p quasinorm) with 0 < p < 1 of a number x ∈ R n is defined as In the 2D case, the gradient of model m is and the symmetric gradient ε(·) reads The minimization problem defining the TGPV regularization term is based on the total generalized variation in image processing (Bredies et al., 2010;Knoll et al., 2011;Zhang et al., 2016).
The advantage of the total generalized variation is that it penalizes both the first-order gradient and high-order gradients of an image or a model, and thus it can effectively avoid the staircase effect compared with the first-order total variation. The advantage of the p -norm penalty over the traditional 1 -norm penalty is observed from the results in compressive-sensing medical imaging (Chartrand, 2009), in which it can handle even sparse data than the 1 -norm for accurate image reconstruction. The p -norm regularization term leads to a nonconvex minimization problem. This introduces difficulties in solving the FWI problem. In the context of compressive sensing, Chartrand (2009) shows that there exist techniques to efficiently solve the p -norm nonconvex problem.
We solve our TGPV-FWI using the similar approaches.
Unlike the conventional TV-regularized FWI, it is impossible to solve the minimization problem in eq. (1) directly. To solve our TGPV-FWI efficiently, we reformulate the minimization in eq. (1) into an alternating-direction minimization problem as by Lin and Huang (2014): which can be decomposed into two interlacing minimization problems: where l is the iteration number in the FWI. That is, the medium parameter m and the auxiliary variable u are updated alternatively in the inversion procedure. The auxiliary variable serves as the prior information in the TGPV-FWI, and it is also updated through iterations.
The first minimization problem in eq. (7a) is a conventional FWI problem with a zeroth-order Tikhonov regularization term. Numerous methods can be adopted to solve this nonlinear minimization problem, such as the conjugate-gradient (CG) method and the limited-memory BFGS (L-BFGS) method (Nocedal and Wright, 2006), etc. Misfit functions other than the simple 2norm-squared waveform difference in eq. (1), as discussed in the Introduction, can also be applied to the first subproblem. We adopt the split-Bregman iteration methodology (Goldstein and Osher, 2009) to solve the second minimization problem.
We show detailed algorithms for solving the second minimization problem in the 2D and 3D cases in Appendices A and B, respectively.
Our TGPV-FWI contains several parameters to be adjusted. We describe how to select these parameters in the following.
The positive parameter λ 1 controls the "strength" of the regularization term in the first minimization problem, and a larger λ 1 leads to a stronger TGPV regularization. In applications, we use the following simple rule to determine λ 1 : where γ is a scaling factor that can be tuned during iterations, and ∇χ(m) represents the gradient of the misfit function w.r.t. some model parameter at certain iteration. · 2 represents the 2norm of a quantity. For different applications, the value of γ can vary. Normally, a value between 0.05 to 0.5 should be suitable for most applications, and larger values of γ are not encouraged. In our numerical tests, we determine different λ 1 s for different model parameters (e.g., V p or V s in elastic-waveform inversion) in each iteration using eq. (8). We find γ = 0.1 can serve as a suitable value.
The method to determine the regularization parameter λ 1 is slightly different from that in such as Lin and Huang (2014). The general principle, however, is in common: the regularization term and the data misfit term should be in decent balance to avoid excessive or insufficient regularization.
There are six parameters in the second minimization problem in eq. (A-11): α 0 , α 1 , η 0 , η 1 , µ and the norm p. The positive parameter λ 2 (or equivalently 1/µ) controls the "strength" of denoising in the second minimization problem, and a larger λ 2 (or equivalently a smaller µ) leads to a stronger smoothing. Note that the smoothing is not a spatial smoothing in a usual sense, such as that of Gaussian spatial filtering; it is a smoothing in the total-variation sense. In extreme cases, µ → 0 can smooth out all the features of a model, while µ → +∞ leaves the model unchanged.
The most important feature of this smoothing is that it simultaneously preserves sharp interfaces and smooth variations of the model in the framework of the high-order total variation.
Therefore, for the second subproblem, only the parameter µ should be tuned for different FWI problems. This indicates that in our TGPV-FWI, in most cases we only need to tune two parameters, the Tikhonov regularization coefficient λ 1 and the regularization parameter µ (or equivalently λ 2 ) as described above, leading to an efficient inversion system.
We summarize our TGPV-FWI algorithm in Appendix C.

Numerical Results
In the following, we use Tikhonov-TV to denote the FWI with the Tikhonov regularization, TV-FWI to denote the FWI with the TV regularization, and TGPV-FWI to denote the FWI with our TGPV regularization.

Synthetic data example I: Checkerboard model
Our TGPV-FWI for acoustic and elastic waves can be used for either large-scale tomography of the Earth using earthquake data or small-scale seismic reflection inversion. For the former case, transmission seismic waves are mostly used. To verify the improved inversion accuracy of our new TGPV-FWI algorithm using mostly transmission signals, we design a checkerboard model with randomly distributed sources within the model. The goal of this numerical test is to verify the applicability of our new method using mostly transmission signals, rather than using Earthquake data.
The model is defined in a region of 970 m×970 m with a grid interval of 10 m in both directions. This checkerboard model is composed of a smoothly-varying background and checkerboard velocity perturbations, as shown in Fig. 1a. The velocity perturbation is approximately 10 percent throughout the model, and the perturbation can be as high as approximately 20 percent in some regions. We use a smooth velocity model without the large-contrast velocity perturbations as displayed in Fig. 1b as the initial velocity model for FWI.
We place a total of 24 sources randomly distributed throughout the model (black stars in The source wavelet is a Ricker wavelet with a center frequency of 20 Hz. We employ the L-BFGS inversion framework and terminate the inversion after 150 iterations. We conduct three tests using the Tikhonov-FWI, TV-FWI and TGPV-FWI.

Synthetic data example II: Modified elastic SEG/EAGE overthrust model
We use a modified elastic SEG/EAGE overthrust model to verify the efficacy of our TGPV-FWI method for elastic-waveform inversion using surface reflection data. In the following numerical tests, we verify our TGPV-FWI method for: (1) noise-free data with adequate numbers of sources and receivers (a regular source/receiver geometry); (2) noise-free, sparse data; (3) noisy data with a regular source/receiver geometry; and (4) noise-free data with dynamic source encoding. We employ a total of 80 equally-spaced sources and 399 equally-spaced receivers at a depth of 50 m in the model. The source interval is 250 m and the receiver interval is 50 m. We use a Ricker wavelet with a center frequency of 8 Hz as the source wavelet. The source is vertical source force acted on the particle velocity wavefield. We employ the CG inversion framework for inversion. The initial P-and S-wave velocity models are shown in Figs. 6a and b, respectively. We also show the region of interest from the true velocity models in Figs. 7a and b for comparison of inversion results. We compare FWI inversion results obtained using the Tikhonov-FWI, TV-FWI and TGPV-FWI.
In the first numerical test, we compare the results obtained using FWI with three different regularization schemes with all data of the 80 sources and the 399 receivers. The Tikhonov-FWI reconstructs smooth P-and S-wave velocity models as shown in Figs. 8a and b, respectively. The TV-FWI improves reconstruction of the interfaces and inversion resolution as depicted in Figs. 8c and d. By contrast, our TGPV-FWI produces best results among the three methods, as displayed in Figs. 8e and f. The method not only improves reconstruction of the interfaces, but also inversion accuracy. In addition, Figs. 8e and f contains fewer inversion artifacts than the other results in Fig. 8. Comparing between the inversion results obtained using the TV-FWI and our TGPV-FWI, we find that, the TV-FWI results contain staircase artifacts in regions between thin layers where seismic velocities vary smoothly in space. By contrast, our TGPV-FWI accurately reconstructs these smoothly varying velocities.
The comparison among the inversion results in Fig. 8 in the region of interest manifests the significant improvements of our TGPV-FWI compared with the other two FWI methods, as shown in Fig. 9.
In the second numerical test, we verify the efficacy of our TGPV-FWI for sparse seismic data acquired using inadequate numbers of sources and receivers. We use seismic data for one third of We conduct the third numerical test using noisy data. Fig. 12a shows the vertical component data of the 40th common-shot gather for the modified elastic SEG/EAGE overthrust model. Fig. 12b depicts the same vertical component data with random noise. The random noise obviously deteriorates the quality of the synthetic data, making many reflections indiscernible from the random noises. We also add the same level of random noise to the horizontal component.
We show the inversion results obtained using the Tikhonov-FWI, the TV-FWI and the TGPV-FWI in Fig. 13 and Fig. 14. In this case, the noises in the data leads to degraded inversion results compared with those in the previous two numerical tests. Nevertheless, our TGPV-FWI results shown in Figs. 14e and f are better than those produced with Tikhonov-FWI (Figs. 14a and b) and TV-FWI (Figs. 14c and d).
In the fourth numerical test, we adopt dynamic source encoding to verify the capability of our TGPV-FWI method. We form six encoded super gathers using the data for 80 sources, and dynamically encode the phase and amplitude of the super gathers with random time delays and random polarity reversal, respectively, over iterations. Fig. 15 and Fig. 16 show the inversion results for the aforementioned three FWI methods. The results show that our TGPV-FWI produces the most accurate inversion results among the three FWI methods.
Finally, we compare the convergences of the data and model misfits for the three FWI methods in Fig. 17. Fig. 17a shows the relative data misfits in the first numerical test using noise-free data with a regular source/receiver geometry. The data misfits of the Tikhonov-FWI and TV-FWI differ from each other insignificantly. By contrast, the data misfit for our TGPV-FWI decreases to a much smaller value than those of the other two methods. The corresponding data misfits in Fig. 17b for the noise-free, sparse data resemble the results in Fig. 17a.
Model misfit is another important indicator for FWI convergence. The model misfits in Fig. 17c and d for FWI with the Tikhonov, TV and TGPV regularizations and the sparse seismic data demonstrate again that our TGPV regularization is the most effective among the three regularization schemes. The P-and S-wave velocity model misfits for the source-encoding FWI test displayed in Fig. 17c and d, respectively, indicate that the model misfits for the source-encoding TGPV-FWI are always smaller than those for the other two inversion methods.

Field data example: Soda Lake geothermal field
We apply our TGPV-FWI method to surface seismic data acquired at the Soda Lake geothermal field in Nevada, USA, and compare the result with those obtained using the Tikhonov-FWI and TV-FWI.  and TGPV-FWI. Comparing these three results, we find that Tikhonov-FWI produces only a lowresolution velocity model without many structural details, particularly in the region beneath the basalt body.
The TV-FWI result in Fig. 18b provides more details compared with that of Tikhonov-FWI in Fig. 18a. The TGPV-FWI result in Fig. 18c has the highest resolution among the three results.
The accuracy of FWI with field seismic data can be validated using the convergence curve of the data misfit. We plot the convergence curves of the data misfits for the Tikhonov-FWI, TV-FWI, and our TGPV-FWI of the field seismic data in Fig. 19, showing that our TGPV-FWI converges fastest and to the smallest value after 100 iterations among the three methods.

Conclusions
We have developed a novel full-waveform inversion method for acoustic and elastic waves using a total generalized p-variation regularization scheme. We decompose the regularized fullwaveform inversion problem into two interlacing minimization problems. The first minimization problem is a conventional full-waveform inversion with a Tikhonov regularization term, and the second minimization problem is a denoising problem using the p -norm total generalized variation.
We have developed an efficient algorithm to solve the second minimization problem based on the split-Bregman iteration method. We have used two synthetic data examples and one field data example to verify the improved accuracy of our new method for both acoustic and elastic full-waveform inversion. Our numerical results demonstrate that our new acoustic and elastic full-waveform inversion with the total generalized p-variation regularization produces accurate inversion results that preserve both sharp interfaces and smooth background velocity variations in the model. In addition, our new inversion method produces reliable inversion results when using sparse seismic data, noisy data, and dynamic source encoding. Our results of synthetic and field seismic data demonstrate that our new inversion method can be used as a robust and accurate tool for subsurface velocity model building, subsurface reservoir characterization, and large-scale tomographic reconstruction using Earthquake data. Appendix A: Solution to the second minimization problem in the 2D case Eq. (7b) is equivalent to the following problem:

Acknowledgments
where h and s are dual variables, µ = 1/λ 2 is a parameter chosen for convenience, and w = [w x , w y ] T .
Using the split-Bregman iteration technique (Goldstein and Osher, 2009;Chartrand, 2009), we reformulate the constrained minimization problem as the following optimization system: whereh ands are auxiliary split-Bregman variables. Note that in the above formulation, l is the index of FWI inversion iteration, while k is the index of the split-Bregman iteration in the second minimization problem. In addition, η 0 and η 1 are two regularization parameters for the constrained minimization problem in eq. (A-12).
We formulate the following four subproblems associated with the variables h, s, u and w from eq. (A-12).
The first subproblem is the minimization of u: The first-order optimality condition of eq. (A-15) is Because the linear system in eq. (A-16) is strictly diagonally dominant, it can be efficiently solved with the Gauss-Seidel iteration: Goldstein and Osher (2009) developed a similar solution for the total-variation image denoising problem.
The second subproblem is the minimization of w: The first-order optimality condition of eq. (A-18) leads to two linear systems: which can also be solved via the Gauss-Seidel iteration because of their strict diagonal dominance.
We use eq. (A-19) to explain the procedure, and the solution to eq. (A-20) can be obtained analogously.
For convenience, we first define and compute an intermediate variable We then obtain The third subproblem is the minimization of h: This problem can be efficiently solved with where the generalized p-shrinkage S p reads (Chartrand et al., 2013) S p ξ , 1 The shrinkage operation is an element-wise operation and therefore can be solved with fairly small computational cost.
The fourth subproblem is the minimization of s: which can also be solved with the generalized p-shrinkage: Lastly, the update of the auxiliary split-Bregman variables is trivial. For example, forh, we And the update of the auxiliarys is similar.

Appendix B: Solution to the second minimization problem in the 3D case
In the 3D case, the TGPV regularization term in eq. (1) reads: where the gradient is and the symmetric gradient ε(·) is We still use a dual-variable alternating-direction minimization strategy to solve the regularized FWI. We transform the second minimization problem to arg min u,w,h,s which also contains four subproblems as is in the 2D case.
The first-order optimality condition for the minimization of variable u is The Gauss-Seidel solution to this linear system is x|i, j,k ) y|i, j,k ) The first-order optimality condition for the minimization of variable w is a system composed of three equations: Solutions to the above three minimization equations can be obtained with the Gauss-Seidel iteration analogously to those for the 2D case. For instance, for the first equation (B-9), we first define and compute two intermediate variables v y,x for ∇ x w y and v z, We then obtain Solutions to the other two equations can be obtained similarly.
The minimization problems for the other variables, including the split-Bregman variables, are trivial to solve.

Appendix C: Algorithm implementation of the TGPV-FWI
We summarize the algorithm of our TGPV-FWI in the following Algorithm 1. Step (3) of this algorithm, n max is the maximum number of the outer split-Bregman iterations, which can also be replaced by a criterion such as ||u (k) − u (k−1) || 2 < ε 0 where ε 0 is a small number; m max is the maximum number of inner split-Bregman iterations, and we find that m max = 2 suffices to produce an accurate result for the second minimization problem in various FWI problems. The computation time of the second minimization problem is trivially small compared to the entire FWI inversion process. Inverted checkerboard velocity models after 150 iterations using ( Figure 17: Data misfits of Tikhonov-FWI, TV-FWI and TGPV-FWI for (a) the first numerical test using all data and for (b) the second numerical test using sparse data. Panels in (c) and (d) show the P-and S-wave velocity model misfits in the second numerical test using sparse data. Panels in (e) and (f) depict the P-and S-wave velocity model misfits in the fourth numerical test using source-encoding inversions.  Figure 19: Comparison of the convergence curves of data misfits for Tikhonov-FWI, TV-FWI and TGPV-FWI of surface seismic data from the Soda Lake geothermal field.