We propose the use of 1 regularization in a wavelet basis for the solution of linearized seismic tomography problems Am = d, allowing for the possibility of sharp discontinuities superimposed on a smoothly varying background. An iterative method is used to find a sparse solution m that contains no more fine-scale structure than is necessary to fit the data d to within its assigned errors. Key words: inverse problem, one-norm, sparsity, tomography, wavelets. Like most geophysical inverse problems, the linearized problem Am = d in seismic tomography is underdetermined, or at best offers a mix of over and underdetermined parameters. It has, therefore, long been recognized that it is important to suppress artefacts that could be falsely interpreted as 'structure' in the earth's interior. Not surprisingly, strategies that yield the smoothest solution m have been dominant in most global or regional tomographic applications; these strategies include seeking global models represented as a low-degree spherical harmonic expansion (Dziewonski et al. 1975; Dziewonski & Woodhouse 1987; Masters et al. 1996) as well as regularization via minimization of the gradient (∇m) or second derivative (∇ 2 m) norm of a dense local parametrization (Constable et al. 1987; Nolet 1987; Spakman & Nolet 1988; VanDecar & Snieder 1994; Trampert & Snieder 1996) . Smooth solutions, however, while not introducing small-scale artefacts, produce a distorted image of the earth through the strong averaging over large areas, thereby making small-scale detail difficult to see, or even hiding it. Sharp discontinuities are blurred into gradual transitions. For example, the inability of global, sphericalharmonic, tomographic models to yield as clear an image of uppermantle subduction zones as produced by more localized studies has long been held against them. and optimize images of upper-mantle slabs to fit physical models of heat diffusion, in an effort to suppress small-scale imaging artefacts while retaining sharp boundaries. Portniaguine & Zhdanov (1999) use a conjugate-gradient method to seek the smallest possible anomalous domain by minimizing a norm based on a renormalized gradient ∇m/(∇m · ∇m + γ 2 ) 1 2 , where γ is a small constant. Like all methods that deviate from a least-squares type of solution, both these methods are non-linear and pose their own problems of practical implementation. From a mathematical point of view, the linear inverse problem Am = d requires regularization whenever A has either very small singular values or a null space. Even small variations in d can then cause large differences in m; regularization restores stability to the problem (Bertero & Boccacci 1998). The choice among the many different regularization techniques depends on the problem, and can be governed, to some extent, by the desire to preserve or emphasize particular features of the (unknown) m that are expected. For instance, if the small singular values of A are mostly associated to highly oscillatory singular functions (as is often the case), and if one is mostly interested in large-scale features of m, then it is natural to introduce a regularizing constraint or penalty term involving ∇m or ∇ 2 m. Or, if one seeks to find solutions m that are 'sparse' with respect to a given basis {ϕ ; = 1, . . . , L}, that is, solutions of the form m = L =1 c ϕ in which only a small fraction of the coefficients c 1 , c 2 , . . . , c L are non-zero, then this can be achieved by a constraint or a penalty term involving c 1 = L =1 | m, ϕ |. This regularization procedure selects solutions that are particularly 'simple' when expressed in terms of the basis {ϕ ; = 1, . . . , L}. The notion that we should seek the 'simplest' model m that fits a measured set of data d to within the assigned errors is intuitively equivalent to the notion that the model should be describable with a small number of parameters. However, clearly, restricting the model to a few low-degree spherical-harmonic or Fourier coefficients, or a few large-scale blocks or tetrahedra, does not necessarily lead to a geophysically plausible solution. In this paper, we investigate whether a multiscale representation based upon wavelets (Daubechies 1992) has enough flexibility to represent the class of models we seek. We propose an 1 -norm regularization method which yields a model m that has a strong tendency to be sparse in a wavelet basis, meaning that it can be faithfully represented by a relatively small number of non-zero wavelet coefficients. This allows for models that vary smoothly in regions of limited coverage without sacrificing any sharp or small-scale features in well-covered regions that are required to fit the data. Our approach is different from an approach briefly suggested by de Hoop & van der Hilst (2005) , in which the mapping between data and model is decomposed in
For example, REF regularized the model in the wavelet domain with 1 -norm constraints, allowing sharp velocity discontinuities (model perturbations) to be superimposed on a smooth background model.
Tomographic inversion using 1-norm regularization of wavelet coefficients
Geophysical Journal International 170 (2007) 359-370
Physics Mathematics