Nonlinear regularization techniques for seismic tomography

https://doi.org/10.1016/j.jcp.2009.10.020Get rights and content

Abstract

The effects of several nonlinear regularization techniques are discussed in the framework of 3D seismic tomography. Traditional, linear, 2 penalties are compared to so-called sparsity promoting 1 and 0 penalties, and a total variation penalty. Which of these algorithms is judged optimal depends on the specific requirements of the scientific experiment. If the correct reproduction of model amplitudes is important, classical damping towards a smooth model using an 2 norm works almost as well as minimizing the total variation but is much more efficient. If gradients (edges of anomalies) should be resolved with a minimum of distortion, we prefer 1 damping of Daubechies-4 wavelet coefficients. It has the additional advantage of yielding a noiseless reconstruction, contrary to simple 2 minimization (‘Tikhonov regularization’) which should be avoided. In some of our examples, the 0 method produced notable artifacts. In addition we show how nonlinear 1 methods for finding sparse models can be competitive in speed with the widely used 2 methods, certainly under noisy conditions, so that there is no need to shun 1 penalizations.

Introduction

Since geophysical inverse problems are almost always underdetermined, regularization techniques are essential to obtain a meaningful solution. Two major classes of techniques exist. The first one, named ‘mollifying’ in the mathematical literature, or ‘optimally localized averaging’ (OLA) in helioseismology, can be traced back to the groundbreaking work of Backus and Gilbert [1], [2] in geophysics. In this approach one searches for the size of an averaging volume that can produce a local average of the model parameter with an acceptable variance. Since this method is computationally very expensive, it has found little application in large-scale geophysical inversions such as seismic tomography. To limit the computational effort, seismic tomographers instead search for a biased (‘damped’) solution. This has the disadvantage of introducing a systematic error – the bias – in lieu of the random error caused by the propagation of data errors. It can be turned into an advantage if the bias actually reflects a justified disposition of the scientist to prefer certain models over others, as long as the data are fit within their error bars.

Simple 2-norm damping, which biases model perturbations towards zero in the absence of information based on the data, is generally a bad choice to regularize the inverse problem for seismic tomography as it tends to introduce structures reflecting ray coverage into the images. For that reason, most tomographers prefer to bias the model towards ‘smooth’ anomalies, in effect trying to forge a compromise between Backus–Gilbert theory and the efficiency of damped inversions. The smoothness of the images has the advantage that large structures become easily visible. Sharp discontinuities, however, are blurred, and smaller structures, even when resolved, may be diluted beyond recognition. Recently, Loris et al. [3] — hereafter referred to as Paper I — introduced a third option for the bias in geophysical inversions: to minimize the 1 norm of the wavelet decomposition. This also biases the model towards zero, but it turns out that such reconstructions always have many (or most) wavelet coefficients exactly equal to zero (i.e. they are sparse). In a synthetic 2D experiment using surface wave data, we showed how structurally coherent features (in a geophysical sense), were more faithfully reproduced using this technique than with a simple 2-norm damping. In addition, as a result of their inherent sparsity, 1 reconstructions exhibit much less noise than their 2 counterparts.

Though Paper I clearly showed the feasibility of wavelet-based 1 regularization it left a number of questions unanswered, in particular which wavelet families work best, how they would perform against more sophisticated 2 norms (e.g. smoothness damping) and whether the computational feasibility as well as the positive conclusions for wavelet regularization do scale up to large, 3D models.

In this paper we therefore aim to refine the original conclusions of Paper I and to enlarge the scope of the investigation. We extend tests to 3D inversions of body wave travel times and investigate the use of different families of wavelets (Haar, D4, and dual tree). We include a comparison with smoothness damping, and with a fourth option named ‘total variation’ damping. Unlike Paper I, we also discuss so-called 0 constrained recovery. In contrast to the 1-norm technique which relies on iterative soft-thresholding, the 0 recovery method uses iterative hard-thresholding of wavelet coefficients [4], [5].

The (salt dome) model that we try to reconstruct here is more realistic than the model in Paper I and includes a wide range of length scales. The problem described in Paper I was also of a very limited size: there were only about 104 degrees of freedom in the reconstructed models. Here we perform 3D reconstructions, and increase the number of degrees of freedom by an order of magnitude to about ∼105. The number of data also increases accordingly to 24000. Our approach is complementary to that of [6] who expand the Fréchet kernels into wavelets to obtain a significant reduction in the memory requirements to store the kernel.

One disadvantage of the 1 norm, 0 and total variation penalties is that they lack the convenient linearity of the more conventional 2-norm minimizations. Making use of recent algorithmic improvements, we do demonstrate that finding a (nonlinear) sparse model reconstruction is not necessarily more expensive, computationwise, than 2-based (linear) reconstructions.

The use of 1 norms in seismic tomography was to the best of our knowledge first proposed in [7] in the form of an iteratively reweighted least squares method (IRLS, see also [8]), but has never found much favor, possibly because the convergence of IRLS was not guaranteed. Besides its use in seismic tomography 1 norms have found application in other geophysical contexts such as deconvolution and interpolation [9], [10], [11], [12]. The use of 1 norms in combination with a carefully chosen basis (such as wavelets) is, however, more recent and is largely inspired by the recent development of compressed sensing [13], [14], [15], that shows that under certain conditions an exact reconstruction can be obtained by solving an 1 regularized inverse problem, provided there is an underlying basis in which the desired model is sparse. This emphasizes the importance of studying different “dictionaries” as we do here. Recently [16] has successfully applied the compressed sensing idea to wavefield reconstruction, albeit on small-scale problems only. The success and promise of compressed sensing has therefore also increased interest in the speed-up of such 1 problems to be able to handle practical geophysical applications (e.g. [17], [18]).

Section snippets

Forward problem formulation

We plan to test the regularization methods on a synthetic data set generated for a salt dome model. Since the main goal of this paper is to evaluate and compare a number of algorithms, numerical efficiency is more important than the wish to have a tomographic problem at hand that is fully realistic. We have thus taken a few shortcuts to be able to run inversions quickly in Matlab on a single processor. However, we took pains to ensure that we would invert for a model that has a large range of

Reconstruction methods

Reconstructing the model m from the data d is done, in principle, by solving the linear system Am=d, where A denotes, as before, the matrix containing the kernels discretized in the voxel basis. This system may contain incompatible equations (due to noise), and at the same time be underdetermined (not enough data to reconstruct all of m).

The problem of incompatible data can be solved by replacing the original problem with the minimization of a data fidelity term:m¯=argminmAm-d2.Here and in

Reconstructions

In this section we present some sample reconstructions using the algorithms mentioned in Section 3, applied to the finite frequency tomography problem described in Section 2. First we will consider a simple checkerboard input model that we will try to reconstruct from incomplete and noisy data. We will also look at a more complicated salt dome model which we obtained from BP America, Inc.

In each case, our procedure will be the following. We start from a known input model minput from which we

Conclusions

In Paper I we showed how a large scale anomaly could be reconstituted even where it was ill resolved because of the selective nature of the wavelet coefficients and the 1 criterion: one wavelet coefficient reconstituting a large, circular, anomaly gave a better optimization than a couple of coefficients reconstituting only the resolved part. With the results of the much more complex salt dome model at hand, we must now conclude that this probably represents more a (lucky) exception than a

Acknowledgments

The authors thank BP America Inc., for kindly providing the 3D salt dome velocity model and several referees for their suggestions that helped improve the manuscript. Part of this research has been supported by the Francqui Foundation (I.L.), the VUB-GOA 062 grant (I.D. and I.L.), the FWO-Vlaanderen grant G.0564.09N (I.D. and I.L.) and NSF grant DMS-0530865 (I.D.).

References (46)

  • S. Chevrot et al.

    Multiscale finite frequency Rayleigh wave tomography of the Kaapvaal craton

    Geophysical Journal International

    (2007)
  • G. Nolet

    Seismic wave propagation and seismic tomography

  • H.L. Taylor et al.

    Deconvolution with the 1 norm

    Geophysics

    (1979)
  • D. Oldenburg et al.

    Recovery of the acoustic impedance from reflection seismograms

    Geophysics

    (1983)
  • F. Santosa et al.

    Linear inversion of band-limited reflection seismograms

    SIAM Journal on Scientific and Statistical Computing

    (1986)
  • M.D. Sacchi et al.

    Interpolation and extrapolation using a high-resolution discrete Fourier transform

    IEEE Transactions on Signal Processing

    (1998)
  • E.J. Candès et al.

    Stable signal recovery from incomplete and inaccurate measurements

    Communications on Pure and Applied Mathematics

    (2006)
  • D. Donoho

    Compressed sensing

    IEEE Transactions on Information Theory

    (2006)
  • D.L. Donoho

    For most large underdetermined systems of linear equations the minimal 1-norm solution is also the sparsest solution

    Communications on Pure and Applied Mathematics

    (2006)
  • F.J. Herrmann et al.

    Non-parametric seismic data recovery with curvelet frames

    Geophysical Journal International

    (2008)
  • G. Hennenfent et al.

    New insights into one-norm solvers from the Pareto curve

    Geophysics

    (2008)
  • E. van den Berg et al.

    Probing the Pareto frontier for basis pursuit solutions

    SIAM Journal on Scientific Computing

    (2008)
  • N. Favier et al.

    Sensitivity kernels for shear wave splitting in transverse isotropic media

    Geophysical Journal International

    (2003)
  • Cited by (54)

    • Parameter identification for continuous point emission source based on Tikhonov regularization method coupled with particle swarm optimization algorithm

      2017, Journal of Hazardous Materials
      Citation Excerpt :

      Neupauer et al. applied Tikhonov regularization to recover the history of water contaminant [18]. Loris et al. applied nonlinear regularization techniques for 3D seismic tomography [22]. They reported that nonlinear L1 method was much more efficient than classical L2 minimization.

    View all citing articles on Scopus
    View full text