Tristan van Leeuwen

Seismic Waveform Inversion by Stochastic Optimization

We explore the use of stochastic optimization methods for seismic waveform inversion. The basic principle of such methods is to randomly draw a batch of realizations of a given misfit function and goes back to the 1950s. The ultimate goal of such an approach is to dramatically reduce the computational cost involved in evaluating the misfit. Following earlier work, we introduce the stochasticity in waveform inversion problem in a rigorous way via a technique called randomized trace estimation. We then review theoretical results that underlie recent developments in the use of stochastic methods for waveform inversion. We present numerical experiments to illustrate the behavior of different types of stochastic optimization methods and investigate the sensitivity to the batch size and the noise level in the data. We find that it is possible to reproduce results that are qualitatively similar to the solution of the full problem with modest batch sizes, even on noisy data. Each iteration of the corresponding stochastic methods requires an order of magnitude fewer PDE solves than a comparable deterministic method applied to the full problem, which may lead to an order of magnitude speedup for waveform inversion in practice.

Robust inversion, dimensionality reduction, and randomized sampling

We consider a class of inverse problems where it is possible to aggregate the results of multiple experiments. This class includes problems where the forward model is the solution operator to linear ODEs or PDEs. The tremendous size of such problems motivates dimensionality reduction techniques based on randomly mixing experiments. These techniques break down, however, when robust data-fitting formulations are used, which are essential in cases of missing data, unusually large errors, and systematic features in the data unexplained by the forward model. We survey robust methods within a statistical framework, and propose a semistochastic optimization approach that allows dimensionality reduction. The efficacy of the methods are demonstrated for a large-scale seismic inverse problem using the robust Students t-distribution, where a useful synthetic velocity model is recovered in the extreme scenario of 60% data missing at random. The semistochastic approach achieves this recovery using 20% of the effort required by a direct robust approach.

Estimating nuisance parameters in inverse problems

Many inverse problems include nuisance parameters which, while not of direct interest, are required to recover primary parameters. The structure of these problems allows efficient optimization strategies—a well-known example is variable projection, where nonlinear least-squares problems which are linear in some parameters can be very efficiently optimized. In this paper, we extend the idea of projecting out a subset over the variables to a broad class of maximum likelihood and maximum a posteriori likelihood problems with nuisance parameters, such as variance or degrees of freedom (d.o.f.). As a result, we are able to incorporate nuisance parameter estimation into large-scale constrained and unconstrained inverse problem formulations. We apply the approach to a variety of problems, including estimation of unknown variance parameters in the Gaussian model, d.o.f. parameter estimation in the context of robust inverse problems, and automatic calibration. Using numerical examples, we demonstrate improvement in recovery of primary parameters for several large-scale inverse problems. The proposed approach is compatible with a wide variety of algorithms and formulations, and its implementation requires only minor modifications to existing algorithms.

Resolution analysis by random probing

We develop and apply methods for resolution analysis in tomography, based on stochastic probing of the Hessian or resolution operators. Key properties of our methods are (i) low algorithmic complexity and easy implementation, (ii) applicability to any tomographic technique, including full-waveform inversion and linearized ray tomography, (iii) applicability in any spatial dimension and to inversions with a large number of model parameters, (iv) low computational costs that are mostly a fraction of those required for synthetic recovery tests, and (v) the ability to quantify both spatial resolution and interparameter trade-offs. Using synthetic full-waveform inversions as benchmarks, we demonstrate that autocorrelations of random-model applications to the Hessian yield various resolution measures, including direction- and position-dependent resolution lengths and the strength of interparameter mappings. We observe that the required number of random test models is around five in one, two, and three dimensions. This means that the proposed resolution analyses are not only more meaningful than recovery tests but also computationally less expensive. We demonstrate the applicability of our method in a 3-D real-data full-waveform inversion for the western Mediterranean. In addition to tomographic problems, resolution analysis by random probing may be used in other inverse methods that constrain continuously distributed properties, including electromagnetic and potential-field inversions, as well as recently emerging geodynamic data assimilation.

Robust Full-waveform Inversion Using the Student’s T-distribution

Full-waveform inversion (FWI) is a computational procedure to extract medium parameters from seismic data. Robust methods for FWI are needed to overcome sensitivity to noise and in cases where modeling is particularly poor or far from the real data generating process. We survey previous robust methods from a statistical perspective, and use this perspective to derive a new robust method by assuming the random errors in our model arise from the Student’s t-distribution. We show that in contrast to previous robust methods, the new method progressively down-weighs large outliers, effectively ignoring them once they are large enough. This suggests that the new method is more robust and suitable for situations with very poor data quality or modeling. Experiments show that the new method recovers as well or better than previous robust methods, and can recover models with quality comparable to standard methods on noise-free data when some of the data is completely corrupted, and even when a marine acquisition mask is entirely ignored in the modeling. The ability to ignore a marine acquisition mask via robust FWI methods offers an opportunity for stochastic optimization methods in marine acquisition.

3D FREQUENCY-DOMAIN SEISMIC INVERSION WITH CONTROLLED SLOPPINESS ∗

Seismic waveform inversion aims at obtaining detailed estimates of subsurface medium parameters, such as the spatial distribution of soundspeed, from multiexperiment seismic data. A formulation of this inverse problem in the frequency domain leads to an optimization problem constrained by a Helmholtz equation with many right-hand sides. Application of this technique to industry-scale problems faces several challenges: First, we need to solve the Helmholtz equation for high wave numbers over large computational domains. Second, the data consist of many independent experiments, leading to a large number of PDE solves. This results in high computational complexity both in terms of memory and CPU time as well as input/output costs. Finally, the inverse problem is highly nonlinear and a lot of art goes into preprocessing and regularization. Ideally, an inversion needs to be run several times with different initial guesses and/or tuning parameters. In this paper, we discuss the requirements of the various compo...

Automatic Migration Velocity Analysis And Multiples

Migration velocity analysis attempts to find a velocity model by focusing a migration image with respect to a redundant coordinate, here a shift at depth between forward and timereversed wavefields. An automatic algorithm requires a wave propagation algorithm and an optimization scheme. A 2D frequency-domain finite-difference code for the two-way wave equation provides the wavefields. A gradient-based minimization scheme should find the velocity model best focuses the migration image. In the presence of multiples, the method may provide incorrect results because there is no mechanism to distinguish between primary and multiple reflections. Two modifications are proposed that may improve the result in the presence of multiples.

Easy implementation of advanced tomography algorithms using the ASTRA toolbox with Spot operators

Mathematical scripting languages are commonly used to develop new tomographic reconstruction algorithms. For large experimental datasets, high performance parallel (GPU) implementations are essential, requiring a re-implementation of the algorithm using a language that is closer to the computing hardware. In this paper, we introduce a new MATLAB interface to the ASTRA toolbox, a high performance toolbox for building tomographic reconstruction algorithms. By exposing the ASTRA linear tomography operators through a standard MATLAB matrix syntax, existing and new reconstruction algorithms implemented in MATLAB can now be applied directly to large experimental datasets. This is achieved by using the Spot toolbox, which wraps external code for linear operations into MATLAB objects that can be used as matrices. We provide a series of examples that demonstrate how this Spot operator can be used in combination with existing algorithms implemented in MATLAB and how it can be used for rapid development of new algorithms, resulting in direct applicability to large-scale experimental datasets.

A modified, sparsity-promoting, Gauss-Newton algorithm for seismic waveform inversion

Images obtained from seismic data are used by the oil and gas industry for geophysical exploration. Cutting-edge methods for transforming the data into interpretable images are moving away from linear approximations and high-frequency asymptotics towards Full Waveform Inversion (FWI), a nonlinear data-fitting procedure based on full data modeling using the wave-equation. The size of the problem, the nonlinearity of the forward model, and ill-posedness of the formulation all contribute to a pressing need for fast algorithms and novel regularization techniques to speed up and improve inversion results. In this paper, we design a modified Gauss-Newton algorithm to solve the PDE-constrained optimization problem using ideas from stochastic optimization and compressive sensing. More specifically, we replace the Gauss-Newton subproblems by randomly subsampled, ℓ1 regularized subproblems. This allows us us significantly reduce the computational cost of calculating the updates and exploit the compressibility of wavefields in Curvelets. We explain the relationships and connections between the new method and stochastic optimization and compressive sensing (CS), and demonstrate the efficacy of the new method on a large-scale synthetic seismic example.

Total Variation Regularization Strategies in Full-Waveform Inversion

We propose an extended full-waveform inversion formulation that includes general convex constraints on the model. Though the full problem is highly nonconvex, the overarching optimization scheme arrives at geologically plausible results by solving a sequence of relaxed and warm-started constrained convex subproblems. The combination of box, total variation, and successively relaxed asymmetric total variation constraints allows us to steer free from parasitic local minima while keeping the estimated physical parameters laterally continuous and in a physically realistic range. For accurate starting models, numerical experiments carried out on the challenging 2004 BP velocity benchmark demonstrate that bound and total variation constraints improve the inversion result significantly by removing inversion artifacts, related to source encoding, and by clearly improved delineation of top, bottom, and flanks of a high-velocity high-contrast salt inclusion. The experiments also show that for poor starting models t...