Regularizing linearized seismic tomography with an ℓ1\ell_1 penalty on the model's wavelet coefficients preserves sharp structures (subducting slabs, rift zones) that classical smoothness constraints blur, and suppresses ray-path-aligned artifacts. Loris, Nolet, Daubechies and Dahlen (2007) demonstrated this on a synthetic two-dimensional surface-wave problem, but identified four coupled obstacles to routine use in large-scale, especially three-dimensional (3-D), tomography: slow convergence of the soft-thresholding iteration for ill-conditioned operators (their experiments needed 300 000 iterations); the memory cost of precomputing and storing the forward operator and its normal matrix ATA\mathbf{A}^T\mathbf{A}; the lack of suitable multiscale bases on the whole sphere and on irregular 3-D grids; and hand-tuned threshold selection with no automated procedure under realistic, correlated data noise. The problem is to remove these obstacles jointly so that the demonstrated advantages of ℓ1\ell_1 wavelet regularization become routinely usable at large 3-D scale: a fast, robust, provably convergent solver, memory-feasible (matrix-free) operator application, a suitable spherical multiscale parametrization or richer dictionaries, and automated, data-driven threshold selection under realistic data covariance. Later work by the same research program built accelerated matrix-free solvers with spherical wavelets (Simons et al. 2011) and applied them to real data from the USArray transportable seismic network (Charléty et al. 2013), but no global-scale, real-data, fully 3-D ℓ1\ell_1 wavelet-regularized inversion with automated parameter selection exists: the largest current linearized global inversion (UNICA25, 2026, 10.57 million arrival times) still uses ℓ2\ell_2 damping and smoothing.

Background

Seismic tomography reconstructs the three-dimensional (3-D) variations of seismic wave speed inside the Earth from measured seismic data, most commonly travel-time delays. In the linearized formulation one solves Am=d\mathbf{A}\mathbf{m}=\mathbf{d}, where m\mathbf{m} is the discretized model (for example, fractional wave-speed perturbations on a grid), d\mathbf{d} collects the delay-time data, and each row of the forward operator A\mathbf{A} is a Fréchet kernel, a sensitivity function; in finite-frequency tomography these kernels describe volumetric, off-path sensitivity of a travel time to Earth structure rather than an infinitely thin ray. Real inversions involve from 10510^5 up to more than 10710^7 data and far more model parameters than the data constrain, so A\mathbf{A} is severely ill-conditioned with a large null space, and a regularization — a prior preference among all models that fit the data to within their assigned errors — must be imposed. The conventional choice is quadratic (ℓ2\ell_2) regularization: damping of model parameters toward zero plus smoothing by penalizing the norm of the spatial gradient or Laplacian, or truncation to low-degree spherical harmonics, the classical global basis of functions on the sphere. Such priors suppress small-scale artifacts but blur genuinely sharp or localized features and can generate artifacts aligned with the ray coverage.

A sparsity-promoting alternative expands m\mathbf{m} in a multiscale wavelet basis — localized functions at a hierarchy of scales, so a sharp feature costs a few large fine-scale coefficients while a smooth background costs a few coarse-scale ones — and penalizes the ℓ1\ell_1 norm (sum of absolute values) of the wavelet coefficients w\mathbf{w}: minimize ∥ASw−d∥22+2τ∥w∥1\|\mathbf{A}\mathbf{S}\mathbf{w}-\mathbf{d}\|_2^2+2\tau\|\mathbf{w}\|_1, with S\mathbf{S} the wavelet synthesis operator and threshold τ\tau. Daubechies, Defrise and De Mol (2004) proved convergence of an iterative algorithm for such problems: each iteration applies the forward and adjoint operators and then a soft threshold, which shrinks every coefficient toward zero by τ\tau and sets small ones exactly to zero. Loris, Nolet, Daubechies and Dahlen (2007) demonstrated this on a synthetic two-dimensional regional surface-wave problem using dual-tree complex wavelets (a redundant, nearly shift-invariant wavelet frame whose 2-D elements distinguish six spatial orientations, avoiding directional bias), fitting the noisy data to the statistically expected misfit, χ2\chi^2 equal to the number of data, with noise modeled as uncorrelated and of constant variance (Σ=σ2I\boldsymbol{\Sigma}=\sigma^2\mathbf{I}) and τ\tau chosen by hand. Their reconstructions kept sharp rift structure while suppressing ray artifacts, outperforming ℓ2\ell_2 regularization of the model or of its wavelet coefficients.

The 2007 paper itself identified what blocked routine adoption at large, especially 3-D, scale. First, the provably convergent soft-threshold iteration is slow for ill-conditioned A\mathbf{A}: constructing trade-off curves required 150 000+150 000 iterations, which the authors called prohibitive for any larger-scale, realistic 3-D application, and their two-step data-shift trick did not remove the problem; they suggested combining an efficient linear method, such as conjugate gradient (a Krylov-type iterative solver for least-squares systems), with an adaptive thresholding scheme, but did not design, analyze, or test such a hybrid. Second, each iteration needs products with ATA\mathbf{A}^T\mathbf{A}, which they precomputed and stored, and doubling the resolution in every direction multiplies the number of stored matrix elements by sixteen. Third, the machinery was 2-D and Cartesian, whereas global 3-D tomography needs multiscale bases on the whole sphere and on irregular grids, including directionally sensitive representations; and if the real Earth is not parsimoniously representable in the chosen wavelet family, richer building blocks must be explored. Fourth, the threshold τ\tau and its component-dependent variants (for example, a lower threshold for the few coarse scaling-function coefficients than for fine-scale wavelet coefficients) were hand-tuned, with no automated, data-driven selection procedure, and with correlated or heterogeneous noise and a non-ideal data covariance Σ\boldsymbol{\Sigma} untreated.

Problem Statement

Make ℓ1\ell_1 wavelet-regularized seismic tomography routinely usable for large-scale, real-data, three-dimensional inversions by jointly resolving all four of the obstacles identified by Loris, Nolet, Daubechies and Dahlen (2007) — none of the four is optional. Here ℓ1\ell_1 wavelet-regularized seismic tomography means linearized travel-time tomography in which the discretized Earth model is expanded in a multiscale wavelet basis and the sum of absolute values of the wavelet coefficients is penalized against data misfit:

  1. Fast, robust solvers: design and analyze algorithms for the ill-conditioned, sparsity-promoting problem that retain the convexity and sparsity-promoting character of the formulation, come with practical convergence behavior and provable convergence guarantees for finite-frequency tomographic operators (forward operators expressing volumetric, off-ray sensitivity of travel times to Earth structure), and remove the need for iteration counts of the hundreds of thousands observed in 2007.

  2. Memory- and CPU-feasible operators: avoid precomputing and storing the forward operator A\mathbf{A} and its normal matrix ATA\mathbf{A}^T\mathbf{A} at realistic 3-D scale, for instance by on-the-fly (matrix-free) application of the forward and adjoint kernel actions and of the wavelet transforms, including their parallel and distributed implementation, so that the approach scales to problems with millions of data and 10610^6–10810^8 model unknowns.

  3. Multiscale parametrization: choose or construct suitable multiscale bases on the whole sphere and on irregular 3-D grids — including directionally sensitive frames and the trade-off between redundancy and compactness, or richer dictionaries if the Earth's true structure is not parsimoniously representable in a wavelet family — integrated with the discretization and with geophysical priors, and quantify which representations actually sparsify realistic Earth structure.

  4. Automated threshold selection: develop principled, data-driven procedures for choosing the threshold τ\tau and component-dependent thresholds (e.g. scaling versus fine-scale wavelet coefficients) that balance data fit against sparsity in practice, under correlated, heterogeneous noise and realistic, non-ideal data covariance matrix Σ\boldsymbol{\Sigma}, without ground truth.

The claimed advantages — sharp-feature recovery, artifact suppression, parsimony (fewer significant wavelet coefficients) at equal data misfit — must be shown to survive at regional-to-global 3-D scale (millions of data and 10610^6–10810^8 model unknowns), demonstrated on real data with resolution and recovery tests and comparison against established damping-and-smoothing inversions.

References

  1. Loris, I., Nolet, G., Daubechies, I., Dahlen, F. A. (2007). Tomographic inversion using ℓ1-norm regularization of wavelet coefficients. Geophysical Journal International 170(1), 359–370. DOI: 10.1111/j.1365-246x.2007.03409.x, https://doi.org/10.1111/j.1365-246x.2007.03409.x
  2. Daubechies, I., Defrise, M., De Mol, C. (2004). An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics 57(11), 1413–1457. DOI: 10.1002/cpa.20042, https://doi.org/10.1002/cpa.20042
  3. Loris, I. (2009). On the performance of algorithms for the minimization of ℓ1-penalized functionals. Inverse Problems 25(3), 035008. DOI: 10.1088/0266-5611/25/3/035008, https://doi.org/10.1088/0266-5611/25/3/035008
  4. Loris, I., Douma, H., Nolet, G., Daubechies, I., Regone, C. (2010). Nonlinear regularization techniques for seismic tomography. Journal of Computational Physics 229(3), 890–905. DOI: 10.1016/j.jcp.2009.10.020, https://doi.org/10.1016/j.jcp.2009.10.020
  5. Simons, F. J., Loris, I., Nolet, G., Daubechies, I. C., Voronin, S., Judd, J. S., Vetter, P. A., Charléty, J., Vonesch, C. (2011). Solving or resolving global tomographic models with spherical wavelets, and the scale and sparsity of seismic heterogeneity. Geophysical Journal International 187(2), 969–988. DOI: 10.1111/j.1365-246x.2011.05190.x, https://doi.org/10.1111/j.1365-246x.2011.05190.x
  6. Charléty, J., Voronin, S., Nolet, G., Loris, I., Simons, F. J., Sigloch, K., Daubechies, I. (2013). Global seismic tomography with sparsity constraints: Comparison with smoothing and damping regularization. Journal of Geophysical Research: Solid Earth 118(9), 4887–4899. DOI: 10.1002/jgrb.50326, https://doi.org/10.1002/jgrb.50326
  7. Chevrot, S., Martin, R., Komatitsch, D. (2012). Optimized discrete wavelet transforms in the cubed sphere with the lifting scheme — implications for global finite-frequency tomography. Geophysical Journal International 191, 1391–1402. DOI: 10.1111/j.1365-246X.2012.05686.X, https://doi.org/10.1111/j.1365-246X.2012.05686.X
  8. Fang, H., Zhang, H. (2014). Wavelet-based double-difference seismic tomography with sparsity regularization. Geophysical Journal International 199(2), 944–955. DOI: 10.1093/gji/ggu305, https://doi.org/10.1093/gji/ggu305
  9. Beck, A., Teboulle, M. (2009). A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems. SIAM Journal on Imaging Sciences 2(1), 183–202. DOI: 10.1137/080716542, https://doi.org/10.1137/080716542
  10. Bianco, M. J., Gerstoft, P. (2018). Travel Time Tomography With Adaptive Dictionaries. IEEE Transactions on Computational Imaging 4(4), 499–511. DOI: 10.1109/TCI.2018.2862644, https://doi.org/10.1109/TCI.2018.2862644
  11. Hung, S.-H., Chen, W.-P., Chiao, L.-Y. (2011). A data-adaptive, multiscale approach of finite-frequency, traveltime tomography with special reference to P and S wave data from central Tibet. Journal of Geophysical Research 116(B6). DOI: 10.1029/2010JB008190, https://doi.org/10.1029/2010JB008190
  12. Chen, S., Wang, G. (2022). Pre-stack Seismic Inversion with L1-2-norm Regularization Via a Proximal DC Algorithm and Adaptive Strategy (preprint). DOI: 10.21203/rs.3.rs-1361049/v1, https://www.researchsquare.com/article/rs-1361049/v1
  13. Nolet, G., Bonnieux, S., Chen, Y. J., Hello, Y., Hosseini, K., Kondo, Y., Namjesnik, D., Obayashi, M., Sigloch, K., Simon, J., Simons, F., Slezak-Oreshina, I., Sugioka, H., Yoshimitsu, J., Yu, Y. (2026). New Global P-wave Tomographic Model of the Earth's Mantle. Seismica 5(1). DOI: 10.26443/seismica.v5i1.1818, https://doi.org/10.26443/seismica.v5i1.1818