Gary F. Margrave

Theory of nonstationary linear filtering in the Fourier domain with application to time-variant filtering

A general linear theory describes the extension of the convolutional method to nonstationary processes. This theory can apply any linear, nonstationary filter, with arbitrary time and frequency variation, in the time, Fourier, or mixed domains. The filter application equations and the expressions to move the filter between domains are all ordinary Fourier transforms or generalized convolutional integrals. Nonstationary transforms such as the wavelet transform are not required. There are many possible applications of this theory including: the one-way propagation of waves through complex media, time migration, normal moveout removal, time-variant filtering, and forward and inverse Q filtering. Two complementary nonstationary filters are developed by generalizing the stationary convolution integral. The first, called nonstationary convolution, corresponds to the linear superposition of scaled impulse responses of a nonstationary filter. The second, called nonstationary combination, does not correspond to such a superposition but is shown to be a linear process capable of achieving arbitrarily abrupt temporal variations in the output frequency spectrum. Both extensions have stationary convolution as a limiting form and, in the discrete case, can be formulated as matrix operations. Fourier transformation shows that both filter types are nonstationary filter integrals in the Fourier domain as well. This result is a generalization of the convolution theorem for stationary signals because, as the filter becomes stationary in one domain, the integral in the other domain collapses to a scalar multiplication. For discrete signals, stationary filters are a matrix multiplication of the input signal spectrum by a diagonal spectral matrix, while nonstationary filters require off-diagonal terms. For quasi-stationary filters, a computational advantage is obtained by computing only the significant terms near the diagonal. Unlike stationary theory, a mixed domain of time and frequency is also possible. In this context, the nonstationary filter is applied simultaneously with the transform from time to frequency or the reverse. Nonstationary convolution becomes a generalized forward Fourier integral and nonstationary combination is a generalized inverse Fourier integral.

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.

The equivalent offset method of prestack time migration

A prestack time migration is presented that is simple, efficient, and provides detailed velocity information. It is based on Kirchhoff prestack time migration and can be applied to both 2-D and 3-D data. The method is divided into two steps: the first is a gathering process that forms common scatterpoint (CSP) gathers; the second is a focusing process that applies a simplified Kirchhoff migration on the CSP gathers, and consists of scaling, filtering, normal moveout (NMO) correction, and stacking. A key concept of the method is a reformulation of the double square‐root equation (of source‐scatterpoint‐receiver traveltimes) into a single square root. The single square root uses an equivalent offset that is the surface distance from the scatterpoint to a colocated source and receiver. Input samples are mapped into offset bins of a CSP gather, without time shifting, to an offset defined by the equivalent offset. The single square‐root reformulation gathers scattered energy to hyperbolic paths on the appropri...

Joint PP and PS seismic inversion

We present a case history of joint inversion of PP and PS reflection seismic data using a weighted stacking technique. Our example comes from Blackfoot Field, owned and operated by PanCanadian Petroleum, in southeastern Alberta, Canada. The exploration target at Blackfoot is a Lower Cretaceous channel system approximately 1.4 km deep (Figure 1). Figure 1. The Glauconitic channel system at Blackfoot Oil Field, Alberta, is a sequence of sand and shale filled valleys incised into Lower Cretaceous and Mississippian carbonates. The Blackfoot interpretation has an upper and lower channel that are prospective and separated by a nonporous lithic channel. These Glauconitic channels, with sand or shale fill, are found throughout the region, and, as there were many episodes of channel formation, can be stacked on top of one another. At Blackfoot, the channel interval is about 40 m thick and 100 m wide. There tends to be good porosity in an upper channel and a lower channel that are separated by a tight, lithic channel. The upper channel, where present, is usually gas-prone, while the lower channel is generally oil-prone. When the pore fluid in the channel sands is a compressible hydrocarbon instead of incompressible water, the bulk compressibility is reduced and this modifies the signature of seismic reflection data. Because pressure waves and shear waves sense different rock and pore-fluid properties, joint use of PP and PS data can provide superior lithologic discrimination. The conversion of one elastic wave, either P or S , into another upon reflection or transmission at an interface is described by the Zoeppritz equations. These equations are algebraically quite complex and it is not practical to reproduce them here. Instead, we will present useful concepts and approximate forms. (We invite the reader to visit our Web site, http://www.crewes.org, and interactively examine the equations using our …

Interpreting channel sands with 3C-3D seismic data

A 3C-3D seismic survey was acquired over the Blackfoot Field (near Strathmore, Alberta, Canada) in 1995. The survey, sponsored by a group of exploration companies, was planned and conducted by the CREWES Project (Department of Geology and Geophysics, The University of Calgary) and Boyd Exploration Consultants. Simultaneously with the surface data acquisition, a five‐level 3-C downhole tool (from Western Atlas International) was deployed in a well, and a 3C-3D VSP was recorded.

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.

Planned seismic imaging using explicit one-way operators

A method is presented to compare the accuracy and computational cost of explicit one-way extrapolation operators as used in seismic imaging. For a given model, accuracy is measured in terms of lateral positioning error, and cost is calculated relative to the cost of the spatial fast Fourier transform. The result is a planned imaging scheme that achieves the greatest accuracy with respect to the velocity model for a fixed cost. To demonstrate, we use a 2D section of the EAGE/SEG salt model and assemble a suite of the most common operators. The data are imaged using each operator individually, and the results are compared to the result from the plan-based algorithm. The planned image is shown to return improved accuracy for no additional cost.

Delineation of Steam Flood Using Seismic Attenuation

The possible use of attenuation measurements for time-lapse seismic monitoring of an EOR steam flood project in Saskatchewan, Canada is investigated. A VSP survey was used to calculate Q. These values were input to a synthetic seismogram attenuation modeling program that showed there should be an observable increase in attenuation after steam injection. Two seismic lines shot at the same location nine years apart were analyzed to see if attenuation anomalies were apparent. The results indicate a strong anomaly on the recent seismic line that is consistent with the location of steam injection. A weaker anomaly on the older line is consistent with the amount of steam injected at that time. Theoretical and laboratory analyses of compressional velocity as a function of changes in temperature, pressure, fluid type and fluid phase suggest there should be a measurable effect on compressional wave amplitude, isochron, and frequency response.

Analyzing the effectiveness of receiver arrays for multicomponent seismic exploration

This paper uses an experimental seismic line recorded with three‐component (3C) receivers to develop a case history demonstrating very little benefit from receiver arrays as compared to point receivers. Two common array designs are tested; they are detrimental to the P‐S wavefield and provide little additional benefit for P‐P data. The seismic data are a 3C 2‐D line recorded at closely spaced (2 m) point receivers over the Blackfoot oil field, Alberta. The 3C receiver arrays are constructed by summing five (one group interval) and ten (two group intervals) point receivers. The shorter array emphasizes signal preservation while the longer array places priority on noise rejection. The effectiveness of the arrays versus the single geophones is compared in both the t−x and f−k domains of common source gathers. The quality of poststack data is also compared by analyzing the f−x spectra for signal bandwidth on both the vertical receiver component (P‐P) and radial receiver component (P‐S) structure stacks produc...