Mauricio D. Sacchi

High‐resolution velocity gathers and offset space reconstruction

We present a high-resolution procedure to reconstruct common-midpoint (CMP) gathers. First, we describe the forward and inverse transformations between offset and velocity space. Then, we formulate an underdetermined linear inverse problem in which the target is the artifacts-free, aperture-compensated velocity gather. We show that a sparse inversion leads to a solution that resembles the infinite-aperture velocity gather. The latter is the velocity gather that should have been estimated with a simple conjugate operator designed from an infinite-aperture seismic array. This high-resolution velocity gather is then used to reconstruct the offset space. The algorithm is formally derived using two basic principles. First, we use the principle of maximum entropy to translate prior information about the unknown parameters into a probabilistic framework, in other words, to assign a probability density function to our model. Second, we apply Bayes’s rule to relate the a priori probability density function (pdf) with the pdf corresponding to the experimental uncertainties (likelihood function) to construct the a posteriori distribution of the unknown parameters. Finally the model is evaluated by maximizing the a posteriori distribution. When the problem is correctly regularized, the algorithm converges to a solution characterized by different degrees of sparseness depending on the required resolution. The solutions exhibit minimum entropy when the entropy is measured in terms of Burg’s definition. We emphasize two crucial differences in our approach with the familiar Burg method of maximum entropy spectral analysis. First, Burg’s entropy is minimized rather than maximized, which is equivalent to inferring as much as possible about the model from the data. Second, our approach uses the data as constraints in contrast with the classic maximum entropy spectral analysis approach where the autocorrelation function is the constraint. This implies that we recover not only amplitude information but also phase information, which serves to extrapolate the data outside the original aperture of the array. The tradeoff is controlled by a single parameter that under asymptotic conditions reduces the method to a damped least-squares solution. Finally, the high-resolution or aperture-compensated velocity gather is used to extrapolate near- and far-offset traces.

Minimum weighted norm interpolation of seismic records

In seismic data processing, we often need to interpolate and extrapolate data at missing spatial locations. The reconstruction problem can be posed as an inverse problem where, from inadequate and incomplete data, we attempt to reconstruct the seismic wavefield at locations where measurements were not acquired.We propose a wavefield reconstruction scheme for spatially band‐limited signals. The method entails solving an inverse problem where a wavenumber‐domain regularization term is included. The regularization term constrains the solution to be spatially band‐limited and imposes a prior spectral shape. The numerical algorithm is quite efficient since the method of conjugate gradients in conjunction with fast matrix–vector multiplications, implemented via the fast Fourier transform (FFT), is adopted. The algorithm can be used to perform multidimensional reconstruction in any spatial domain.

Latest views of the sparse Radon transform

The Radon transform (RT) suffers from the typical problems of loss of resolution and aliasing that arise as a consequence of incomplete information, including limited aperture and discretization. Sparseness in the Radon domain is a valid and useful criterion for supplying this missing information, equivalent somehow to assuming smooth amplitude variation in the transition between known and unknown (missing) data. Applying this constraint while honoring the data can become a serious challenge for routine seismic processing because of the very limited processing time available, in general, per common midpoint. To develop methods that are robust, easy to use and flexible to adapt to different problems we have to pay attention to a variety of algorithms, operator design, and estimation of the hyperparameters that are responsible for the regularization of the solution. In this paper, we discuss fast implementations for several varieties of RT in the time and frequency domains. An iterative conjugate gradient algorithm with fast Fourier transform multiplication is used in all cases. To preserve the important property of iterative subspace methods of regularizing the solution by the number of iterations, the model weights are incorporated into the operators. This turns out to be of particular importance, and it can be understood in terms of the singular vectors of the weighted transform. The iterative algorithm is stopped according to a general cross validation criterion for subspaces. We apply this idea to several known implementations and compare results in order to better understand differences between, and merits of, these algorithms.

Simultaneous seismic data denoising and reconstruction via multichannel singular spectrum analysis

We present a rank reduction algorithm that permits simultaneous reconstruction and random noise attenuation of seismic records. We based our technique on multichannel singular spectrum analysis (MSSA). The technique entails organizing spatial data at a given temporal frequency into a block Hankel matrix that in ideal conditions is a matrix of rank k , where k is the number of plane waves in the window of analysis. Additive noise and missing samples will increase the rank of the block Hankel matrix of the data. Consequently, rank reduction is proposed as a means to attenuate noise and recover missing traces. We present an iterative algorithm that resembles seismic data reconstruction with the method of projection onto convex sets. In addition, we propose to adopt a randomized singular value decomposition to accelerate the rank reduction stage of the algorithm. We apply MSSA reconstruction to synthetic examples and a field data set. Synthetic examples were used to assess the performance of the method in two...

A Bayes tour of inversion: A tutorial

It is unclear whether one can (or should) write a tutorial about Bayes. It is a little like writing a tutorial about the sense of humor. However, this tutorial is about the Bayesian approach to the solution of the ubiquitous inverse problem. Inasmuch as it is a tutorial, it has its own special ingredients. The first is that it is an overview; details are omitted for the sake of the grand picture. In fractal language, it is the progenitor of the complex pattern. As such, it is a vision of the whole. The second is that it does, of necessity, assume some ill‐defined knowledge on the part of the reader. Finally, this tutorial presents our view. It may not appeal to, let alone be agreed to, by all.

Least‐squares wave‐equation migration for AVP/AVA inversion

We present an acoustic migration/inversion algorithm that uses extended double‐square‐root wave‐equation migration and modeling operators to minimize a constrained least‐squares data misfit function (least‐squares migration). We employ an imaging principle that allows for the extraction of ray‐parameter‐domain common image gathers (CIGs) from the propagated seismic wavefield. The CIGs exhibit amplitude variations as a function of half‐offset ray parameter (AVP) closely related to the amplitude variation with reflection angle (AVA). Our least‐squares wave‐equation migration/inversion is constrained by a smoothing regularization along the ray parameter. This approach is based on the idea that rapid amplitude changes or discontinuities along the ray parameter axis result from noise, incomplete wavefield sampling, and numerical operator artifacts. These discontinuities should therefore be penalized in the inversion.The performance of the proposed algorithm is examined with two synthetic examples. In the first...

Accurate interpolation with high‐resolution time‐variant Radon transforms

It is well known that a sparse hyperbolic Radon transform (RT) can be used to extend the aperture of aperture limited data, filter noise, and fill gaps. In the same manner, an elliptical RT can achieve similar results when applied to slant stack sections. A problem with these transformations is that they have a time-variant kernel that results in slow implementation. By defining a model space in terms of an irregularly sampled velocity space to minimize the number of unknowns during the inversion and using sparse matrices, however, the computation time can be kept low enough for practical application. We implement hyperbolic and elliptical time domain RTs by using inversion via weighted conjugate gradient methods with a sparseness constraint. The hyperbolic RT performs accurate interpolation in common-midpoint (CMP) gathers, while the elliptical RT attenuates sampling artifacts in slant stack sections obtained from CMP gathers with poor sampling and gaps.

Fourier reconstruction of nonuniformly sampled, aliased seismic data

There are numerous methods for interpolating uniformly sampled, aliased seismic data, but few can handle the combination of nonuniform sampling and aliasing. We combine the principles of Fourier reconstruction of nonaliased, nonuniformly sampled data with the ideas of frequency-wavenumber (f-k) interpolation of aliased, uniformly sampled data in a new two-step algorithm. In the first step, we estimate the Fourier coefficients at the lower nonaliased temporal frequencies from the nonuniformly sampled data. The coefficients are then used in the second step as an a priori model to distinguish between aliased and nonaliased energy at the higher, aliased temporal frequencies. By using a nonquadratic model penalty in the inversion, both the artifacts in the Fourier domain from nonuniform sampling and the aliased energy are suppressed. The underlying assumption is that events are planar; therefore, the algorithm is applied to seismic data in overlapping spatiotemporal windows.

Multistep autoregressive reconstruction of seismic records

Linear prediction filters in the f-x domain are widely used to interpolate regularly sampled data. We study the problem of reconstructing irregularly missing data on a regular grid using linear prediction filters. We propose a two-stage algorithm. First, we reconstruct the unaliased part of the data spectrum using a Fourier method (minimum-weighted norm interpolation). Then, prediction filters for all the frequencies are extracted from the reconstructed low frequencies. The latter is implemented via a multistep autoregressive (MSAR) algorithm. Finally, these prediction filters are used to reconstruct the complete data in the f-x domain. The applicability of the proposed method is examined using synthetic and field data examples.

High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution

Three-term AVO inversion can be used to estimate P-wave velocity, S-wave velocity, and density perturbations from reflection seismic data. The density term, however, exhibits little sensitivity to amplitudes and, therefore, its inversion is unstable. One way to stabilize the density term is by including a scale matrix that provides correlation information between the three unknown AVO parameters. We investigate a Bayesian procedure to include sparsity and a scale matrix in the three-term AVO inversion problem. To this end, we model the prior distribution of the AVO parameters via a Trivariate Cauchy distribution. We found an iterative algorithm to solve the Bayesian inversion and, in addition, comparisons are provided with the classical inversion approach that uses a Multivariate Gaussian prior. It is important to point out that the Multivariate Gaussian prior allows us to include the correlation of the AVO parameters in the solution of the inverse problem. The Trivariate Cauchy prior not only permits us ...