share
1 Introduction
Full-waveform inversion can achieve high resolution of subsurface velocity reconstruction where the target is shallow and well illuminated by refracted seismic energy from wide-offset surveys Sirgue et al. (2010). However, illumination of deeper targets with refracted energy may require extra wide offset survey acquisitions, or otherwise suffer from poor constraining of deeper model sections. Mathematically, this is a manifestation of the fact that the full-waveform inversion is a mixed-determined problem, with shallow areas of the model overdetermined by the abundance of data, and deeper areas affected by poor resolution and spurious positioning errors. While geological priors such as well tie-ins may provide useful additional constraints in areas of poor illumination, regularization of the nonlinear inversion problems arising in full-waveform inversion is a well-established mathematical technique for dealing with underdetermined problems and noisy data. In particular, it is recognized that the total variation (TV) regularization promotes sparsity of model gradients, acting as an “edge preserving” constraint complementing or outweighing data fitting in problematic areas Anagaw and Sacchi (2011); Guitton (2012). However, and TV regularized optimization problems are difficult to solve, and the development of efficient numerical solution techniques is a subject of active ongoing research, see e.g. Boyd et al. (2011).
In this work we proposed a formulation of the full-waveform inversion as a problem of constrained optimization, and solve it using the iterative Bregman regularization technique, see e.g. Osher et al. (2005). We demonstrate advantages of the proposed method over unconstrained regularization. The paper concludes with an example of application to the Marmousi synthetic with added noise.
2 Method
We begin with the standard formulation of FWI as an unconstrained nonlinear least-squares fitting problem Virieux and Operto (2009):
| (1) |
where is the observed data, is the model (i.e., acoustic slowness) and is the nonlinear forward modeling operator. Problem (1) can be solved in either time or frequency domain, with either approach having its advantages Virieux and Operto (2009). Formulation (1) equally weights all data, resulting in better illuminated areas being better constrained. Weighted least-squares and priors may help improve recovery of deeper sections, however, the problem still remains underdetermined at greater depths or acquisition blind spots.
In this work we explore an alternative formulation of FWI as a constrained optimization problem,
| (2) | ||||
where
| (3) |
is the length of spatial slowness gradient, and is the total variation seminorm of Triebel (2006). Note that in (2) we use an equality constraint for data fitting, which is neither desirable nor realistic in applications to field data. Indeed, formulation (2) is suitable for strictly underdetermined problems, and the equality constraint can be enforced only for noise-free data. Since the full-waveform inversion problem is mixed-determined and field data are always noisy, we propose to solve the following constrained optimization problem:
| (4) | ||||
where is a weighting function that is above zero only in areas of poor resolution or illumination at depth. In a practical solution algorithm, solving (4) is equivalent to solving (2) with suitable stopping criteria when the desired data misfit is achieved.
We solve (4) using the iterative Bregman regularization technique proposed by Osher et al. (2005). Starting from
and given , we iteratively compute slowness as the solution of the following unconstrained TV-regularized problem
| (5) |
where an element of the subgradient of is computed as
| (6) |
Note that (5) describes a TV-regularized inversion Rudin et al. (1992) that may yield an edge-preserving or “blocky” approximation to the solution of the full-waveform inversion problem (1). The first two terms in (5) are known as the Bregman distance Bregman (1967), that is equivalent to extracting an approximately quadratic function from the regularization term. We use the regularization approach based on solving the single unconstrained problem similar to (5) with in our work on TV-regularized time-lapse FWI Maharramov and Biondi (2014). It is important to note, however, that in this work, solution of an unconstrained TV-regularized problem (5) is just a single iteration of an algorithm for solving the constrained problem (4). We solve problem (5) using nonlinear conjugate gradients (NCG) algorithm Nocedal and Wright (2006), with the smoothed TV regularization term
| (7) |
where is chosen as a threshold for realistic values of the slowness.
Choice of the regularization parameter in (5) can be based on achieving better conditioning of problem (5), and unlike the traditional penalty function/continuation methods for solving (2), does not increase between iterations Goldstein and Osher (2009); Cai et al. (2010). Furthermore, iterations (5,6) can be shown to converge to a solution of (2) (or a solution of (4) for some for noisy data) regardless of the value of , as can be demonstrated by a trivial extension of the technique of Cai et al. (2010). However, the rate of convergence does depend on the value of the regularization parameter , making the application of Bregman regularization to some nonlinear operators problematic. However, our experiments indicate that the value of the regularization parameter chosen to improve the convergence of NCG for (5) results in good overall convergence of Bregman iterations.
Iterations are stopped when the data misfit reaches a desired value of Osher et al. (2005); Cai et al. (2010). In practical applications where may not be known a-priori, iterations may continue until the effects of overfitting start exceeding the edge-preserving effects of the regularization term (7). Note that instead of using the NCG to solve (5) with the smoothed regularization term (7), problem (5) can be solved using split Bregman method that only requires iterative solution of nonlinear least squares problems and soft thresholding Goldstein and Osher (2009). However, our numerical experiments indicate that the NCG applied to the smoothed (7) has equivalent performance and accuracy.
3 Numerical examples
We apply the method to the synthetic dataset used in Maharramov and Biondi (2014), generated for the Marmousi velocity model over a 384122 grid with a 24 m grid spacing. The inversion is carried out in the frequency domain for 3.0, 3.6, 4.3, 5.1, 6.2, 7.5, 9.0, 10.8, 12.8, and 15.5 Hz with time-domain forward modeling Sirgue et al. (2010)
. The frequencies are chosen based on the estimated offset to depth range of the data
Sirgue and Pratt (2004). The acquisition has 192 shots at a depth of 16 m with a 48 m spacing, and 381 receivers at a depth of 15 m with a 24 m spacing. The minimum offset is 48 m. The source function is a Ricker wavelet centered at 10.1 Hz. Absorbing boundary conditions are applied along the entire model boundary, including the surface (thus suppressing multiples). A smoothed true model shown in Maharramov and Biondi (2014) is used as a starting model for the inversion. The smoothing is performed using a triangular filter with a 20-sample half-window in both vertical and horizontal directions. Random Gaussian noise is added to the noise-free synthetic data to produce a noisy dataset with 7 dB signal-to-noise ratio. The result of model inversion from the 7 dB SNR synthetic data is shown in Figure LABEL:fig:base. Up to 10 iterations of the nonlinear conjugate gradients algorithm Nocedal and Wright (2006) are performed for each frequency. Neither regularization nor model priors are used. basewidth= Inversion of 7dB SnR synthetic using the unregularized FWI, 10 iterations per frequency. Figure LABEL:fig:cfwi shows the results of solving the proposed constrained FWI (4). Only 5 NCG iterations were used for solving each problem (5), with only two outer (Bregman) iterations, resulting in roughly the same compute time as in our standard FWI experiment shown in Figure LABEL:fig:base (10 gradient evaluations using the adjoint state method). The weighting function was set to 1 below 2100 m and zero above 2000 m, thus the regularization is applied to less constrained areas. cfwiwidth= Inversion of 7dB SnR synthetic by solving constrained TV-regularized FWI (4) using Bregman iterative procedure, 10 iterations per frequency. Note that while the shallow parts are similar to Figure LABEL:fig:base, deeper sections below 2km are more focused, and the poorly illuminated and mispositioned intervals in the left part of the model have been improved. Our results in Figure LABEL:fig:cfwi indicate that CFWI has improved the deepest section of the model while matching the standard FWI in more shallow well-constrained areas. The result of Figure LABEL:fig:cfwi is closer to the clean synthetic inversion shown in Maharramov and Biondi (2014), and has better delineated interfaces.4 Conclusions
We have proposed a new formulation of FWI as a constrained optimization problem (CFWI), and demonstrated the CFWI to be a viable technique for improving depth resolution and accuracy of FWI. Application of Bregman iterative regularization provides a computationally efficient solution method for CFWI that can be easily built on top of the existing solvers. Application of CFWI to field data will be the subject of future work.
5 Acknowledgments
The authors would like to thank Jon Claerbout and Stewart Levin for a number of useful discussions, and the Stanford Center for Computational Earth and Environmental Sciences for providing computing resources.
References
- Anagaw and Sacchi (2011) Anagaw, A. Y., and M. D. Sacchi, 2011, Full waveform inversion with total variation regularization: Presented at the Recovery-CSPG CSEG CWLS Convention.
-
Boyd et al. (2011)
Boyd, S., N. Parikh, E. Chu, B. Peleato, and J. Eckstein, 2011, Distributed optimization and statistical learning via the alternating direction method of multipliers: Foundations and Trends in Machine Learning,
3, 1–122. - Bregman (1967) Bregman, L., 1967, The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming: {USSR} Computational Mathematics and Mathematical Physics, 7, 200 – 217.
- Cai et al. (2010) Cai, J., S. Osher, and Z. Shen, 2010, Split bregman methods and frame based image restoration: Multiscale Modeling & Simulation, 8, 337–369.
- Goldstein and Osher (2009) Goldstein, T., and S. Osher, 2009, The split bregman method for l1-regularized problems: SIAM Journal on Imaging Sciences, 2, 323–343.
- Guitton (2012) Guitton, A., 2012, Blocky regularization schemes for full-waveform inversion★: Geophysical Prospecting, 60, 870–884.
- Maharramov and Biondi (2014) Maharramov, M., and B. Biondi, 2014, Robust joint full-waveform inversion of time-lapse seismic datasets with total-variation regularization: SEP Report, 155, 199–208.
- Nocedal and Wright (2006) Nocedal, J., and S. J. Wright, 2006, Numerical optimization: Springer.
- Osher et al. (2005) Osher, S., M. Burger, D. Goldfarb, J. Xu, and W. Yin, 2005, An iterative regularization method for total variation-based image restoration: Multiscale Modeling & Simulation, 4, 460–489.
- Rudin et al. (1992) Rudin, L. I., S. Osher, and E. Fatemi, 1992, Nonlinear total variation based noise removal algorithms: Phys. D, 60, 259–268.
- Sirgue et al. (2010) Sirgue, L., O. Barkved, J. Dellinger, J. Etgen, U. Albertin, and J. Kommedal, 2010, Full waveform inversion: the next leap forward in imaging at Valhall: First Break, 28, 65–70.
- Sirgue and Pratt (2004) Sirgue, L., and R. Pratt, 2004, Efficient waveform inversion and imaging: A strategy for selecting temporal frequencies: Geophysics, 69, 231–248.
- Triebel (2006) Triebel, H., 2006, Theory of function spaces: Birkhauser.
- Virieux and Operto (2009) Virieux, J., and S. Operto, 2009, An overview of full-waveform inversion in exploration geophysics: Geophysics, 74, no. 6, WCC1–WCC26.
Comments
There are no comments yet.