Michael P. Lamoureux

Gabor deconvolution: Estimating reflectivity by nonstationary deconvolution of seismic data

We have extended the method of stationary spiking deconvolution of seismic data to the context of nonstationary signals in which the nonstationarity is due to attenuation processes. As in the stationary case, we have assumed a statistically white reflectivity and a minimum-phase source and attenuation process. This extension is based on a nonstationary convolutional model, which we have developed and related to the stationary convolutional model. To facilitate our method, we have devised a simple numerical approach to calculate the discrete Gabor transform, or complex-valued time-frequency decomposition, of any signal. Although the Fourier transform renders stationary convolution into exact, multiplicative factors, the Gabor transform, or windowed Fourier transform, induces only an approximate factorization of the nonstationary convolutional model. This factorization serves as a guide to develop a smoothing process that, when applied to the Gabor transform of the nonstationary seismic trace, estimates the magnitude of the time-frequency attenuation function and the source wavelet. By assuming that both are minimum-phase processes, their phases can be determined. Gabor deconvolution is accomplished by spectral division in the time-frequency domain. The complex-valued Gabor transform of the seismic trace is divided by the complex-valued estimates of attenuation and source wavelet to estimate the Gabor transform of the reflectivity. An inverse Gabor transform recovers the time-domain reflectivity. The technique has applications to synthetic data and real data.

Gabor deconvolution of seismic data for source waveform and Q correction

Summary We present a novel approach to nonstationary seismic deconvolution using the Gabor transform. This nonstationary transform represents a signal as a superposition of sinusoids that are localized by time-shifted windows. The resulting time-frequency decomposition is a suite of local Fourier transforms that facilitates nonstationary spectral analysis or filtering. In a result that generalizes the seismic convolutional model, we show that the Gabor transform of a nonstationary seismic signal is the product of source signature, Q filter, and reflectivity effects. We use this spectral factorization theorem as a basis for a new deconvolution algorithm in the Gabor domain. We estimate the Gabor spectrum of the underlying reflectivity directly from the Gabor spectrum of an attenuated seismic signal. Tests on synthetic and real data show that our method works well and combines the effects of sourcesignature inversion and a data-driven inverse Q filter. In comparison with a stationary Wiener deconvolution, our Gabor deconvolution is similar within the Wiener design gate and superior elsewhere.

Phase‐shift time‐stepping for reverse‐time migration

As a result of the numerical performance of finite-difference operators, reverse-time migration (RTM) produces images which are typically low frequency or require large computational resources. We consider an alternative to wavefield propagation with finite differences, a two-way high-fidelity time-stepping equation based on the Fourier transform which is exact for homogeneous media if an aliasing condition is met. The technique is adapted to variable velocity using a localized Fourier transform (Gabor transform). The feasibility of using the time-stepping equation for RTM is demonstrated by studying its stability properties, and by migrating the Marmousi data set. We show that a high frequency wavefield can be time-stepped with no loss of frequency content and with a much larger time step than is commonly used.

Stieltjes representation of the 3D Bruggeman effective medium and Padé approximation

The paper deals with Bruggeman effective medium approximation (EMA) which is often used to model effective complex permittivity of a two-phase composite. We derive the Stieltjes integral representation of the 3D Bruggeman effective medium and use constrained Pade approximation method introduced in [39] to numerically reconstruct the spectral density function in this representation from the effective complex permittivity known in a range of frequencies. The problem of reconstruction of the Stieltjes integral representation arises in inverse homogenization problem where information about the spectral function recovered from the effective properties of the composite, is used to characterize its geometric structure. We present two different proofs of the Stieltjes analytical representation for the effective complex permittivity in the 3D Bruggeman effective medium model: one proof is based on direct calculation, the other one is the derivation of the representation using Stieltjes inversion formula. We show that the continuous spectral density in the integral representation for the Bruggeman EMA model can be efficiently approximated by a rational function. A rational approximation of the spectral density is obtained from the solution of a constrained minimization problem followed by the partial fractions decomposition. We show results of numerical rational approximation of Bruggeman continuous spectral density and use these results for estimation of fractions of components in a composite from simulated effective permittivity of the medium. The volume fractions of the constituents in the composite calculated from the recovered spectral function show good agreement between theoretical and predicted values.

Gabor deconvolution revisited

Gabor deconvolution has been updated and experience has been obtained on real data. The updates are (1) a new method of spectral smoothing called hyperbolic smoothing, (2) a Gabor transform using compactly supported windows that improves run times by one to two orders of magnitude (3) a post-deconvolution time-variant bandpass filter whose maximum frequency tracks along a hyperbola in the timefrequency plane. We discuss the technical details of these improvements and present a data example.

Inverting the Pascal Matrix Plus One

Immediately we notice this is basically the same matrix, except with every other subdiagonal multiplied by a factor of negative one. Obviously this result generalizes to any size Pascal matrix, and there is something so beautiful and natural about this result that it hardly needs a proof (see Call and Velleman [2] for a particularly elegant demonstration). In fact, it is interesting to note that there is no reason to stop at a finite matrix, for one can extend the Pascal matrix to an infinite lower triangular matrix and verify row-by-row that the inverse is given as

Improving explicit seismic depth migration with a stabilizing Wiener filter and spatial resampling

We present a new approach to the design and implementation of explicit wavefield extrapolation for seismic depth migration in the space-frequency domain. Instability of the wavefield extrapolation operator is addressed by splitting the operator into two parts, one to control phase accuracy and a second to improve stability. The first partial operator is simply a windowed version of the exact operator for a half step. The second partial operator is designed, using the Wiener filter method, as a band-limited, least-squares inverse of the first. The final wavefield extrapolation operator for a full step is formed as a convolution of the first partial operator with the complex conjugate of the second. This resulting wavefield extrapolation operator can be designed to have any desired length and is generally more stable and more accurate than a simple windowed operator of similar length. Additional stability is gained by reducing the amount of evanescent filtering and by spatially downsampling the lower tempor...

An FDTD scheme on a face-centered-cubic (FCC) grid for the solution of the wave equation

A method is proposed to improve the numerical dispersion characteristics for simulations of the scalar wave equation in 3D using the FDTD method. The improvements are realized by choosing a face-centered-cubic (FCC) grid instead of the typical Cartesian (Yee) grid, which exhibits non-physical distortions of the wavefront due to the FD stencil. FCC grids are the logical extension of hexagonal grids in 2D, and have been shown previously to provide optimal sampling of space based on close packing of spheres (highest density). The difference equations are developed for the wave equation on this alternative grid, and the dispersion relationship and stability for grids of equal and non-equal aspect ratios are derived. A comparison is made between FCC and Cartesian formulations, based upon having an equal volume density of gridpoints in each method (i.e. the computational storage requirements of each method would be the same for the same simulated space). The comparison shows that the FCC grid exhibits a much more isotropic dispersion relation than the Cartesian grid of equivalent density. Furthermore, for an equivalent density, the FCC method has a more relaxed stability criterion by a factor of approximately 1.35, resulting in a further reduction in computational resources.

The poorman's transform: approximating the Fourier transform without multiplication

A time-domain to frequency-domain transformation for sampled signals which is computed with only additions and trivial complex multiplications is described. This poormans transform is an approximation to the usual Fourier transform, obtained by quantizing the Fourier coefficients to the four values (+or-1, +or-j), and is especially useful when multiplication is expensive. For the general case of an N-point quantization, an analytic formula is given for the error in the approximation, which involves only contributions from aliased harmonics. Continuous-time signals are considered; in this case the approximation is exact for bandlimited signals. >

Numerical decomposition of a convex function

Given then×p orthogonal matrixA and the convex functionf:Rn→R, we find two orthogonal matricesP andQ such thatf is almost constant on the convex hull of ± the columns ofP, f is sufficiently nonconstant on the column space ofQ, and the column spaces ofP andQ provide an orthogonal direct sum decomposition of the column space ofA. This provides a numerically stable algorithm for calculating the cone of directions of constancy, at a pointx, of a convex function. Applications to convex programming are discussed.