Nonlinear regularization techniques for seismic tomography
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 -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 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 -norm damping. In addition, as a result of their inherent sparsity, reconstructions exhibit much less noise than their counterparts.
Though Paper I clearly showed the feasibility of wavelet-based regularization it left a number of questions unanswered, in particular which wavelet families work best, how they would perform against more sophisticated 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 constrained recovery. In contrast to the -norm technique which relies on iterative soft-thresholding, the 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 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 norm, and total variation penalties is that they lack the convenient linearity of the more conventional -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 -based (linear) reconstructions.
The use of 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 norms have found application in other geophysical contexts such as deconvolution and interpolation [9], [10], [11], [12]. The use of 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 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 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: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 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 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)
- et al.
Iterative hard thresholding for compressed sensing
Applied and Computational Harmonic Analysis
(2009) - et al.
Fast solutions of large sparse linear systems: application to seismic travel time tomography
Journal of Computational Physics
(1988) - et al.
Adaptive total variation image deblurring: a majorization–minimization approach
Signal Processing
(2009) - et al.
Wavelets on the interval and fast wavelet transforms
Journal for Applied and Computational Harmonic Analysis
(1993) Complex wavelets for shift invariant analysis and filtering of signals
Applied and Computational Harmonic Analysis
(2001)- et al.
Coordinate and subspace optimization methods for linear least squares with non-quadratic regularization
Applied and Computational Harmonic Analysis
(2007) - et al.
Numerical applications of a formalism for geophysical inverse problems
Geophysical Journal of the Royal Astronomical Society
(1967) - et al.
Uniqueness on the inversion of inaccurate gross Earth data
Philosophical Transactions of the Royal Society of London
(1970) - et al.
Tomographic inversion using -norm regularization of wavelet coefficients
Geophysical Journal International
(2007) - et al.
Iterative thresholding for sparse approximations
Journal of Fourier Analysis and Applications
(2008)
Multiscale finite frequency Rayleigh wave tomography of the Kaapvaal craton
Geophysical Journal International
Seismic wave propagation and seismic tomography
Deconvolution with the norm
Geophysics
Recovery of the acoustic impedance from reflection seismograms
Geophysics
Linear inversion of band-limited reflection seismograms
SIAM Journal on Scientific and Statistical Computing
Interpolation and extrapolation using a high-resolution discrete Fourier transform
IEEE Transactions on Signal Processing
Stable signal recovery from incomplete and inaccurate measurements
Communications on Pure and Applied Mathematics
Compressed sensing
IEEE Transactions on Information Theory
For most large underdetermined systems of linear equations the minimal -norm solution is also the sparsest solution
Communications on Pure and Applied Mathematics
Non-parametric seismic data recovery with curvelet frames
Geophysical Journal International
New insights into one-norm solvers from the Pareto curve
Geophysics
Probing the Pareto frontier for basis pursuit solutions
SIAM Journal on Scientific Computing
Sensitivity kernels for shear wave splitting in transverse isotropic media
Geophysical Journal International
Cited by (54)
Parallel three-dimensional electrical capacitance data imaging using a nonlinear inversion algorithm and L<sup>p</sup> norm-based model regularization
2018, Measurement: Journal of the International Measurement ConfederationParameter identification for continuous point emission source based on Tikhonov regularization method coupled with particle swarm optimization algorithm
2017, Journal of Hazardous MaterialsCitation 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.
Adaptive Spectral Inversion for inverse medium problems
2023, Inverse ProblemsPHYSICS-GUIDED UNSUPERVISED DEEP-LEARNING SEISMIC INVERSION WITH UNCERTAINTY QUANTIFICATION
2023, Journal of Seismic Exploration