Wavefield extrapolation in the domain provides a tool for depth migration with strong lateral variations in velocity. Implicit formulations of depth extrapolation have several advantages over explicit methods. However, the simple 3-D extension of conventional 2-D wavefield extrapolation by implicit finite-differencing requires the inversion of a 2-D convolution matrix which is computationally difficult. In this paper, we solve the 45 wave equation with helical boundary conditions on one of the spatial axes. These boundary conditions reduce the 2-D convolution into an equivalent 1-D filter operation. We then factor this 1-D filter into causal and anti-causal parts using an extension of Kolmogoroff's spectral factorization method, and invert the convolution operator efficiently by 1-D recursive filtering. We include lateral variations in velocity by factoring spatially variable filters, and non-stationary deconvolution. The helical boundary conditions allow the 2-D convolution matrix to be inverted directly without the need for splitting approximations, with a cost that scales linearly with the size of the model space. Using this methodology, a whole range of implicit depth migrations may now be feasible in 3-D.

Kirchhoff imaging beyond aliasing (ps 1755K) (src 2618K)

I present a method for anti-aliasing Kirchhoff imaging operators that improves the resolution of the image by properly imaging aliased components of the recorded data that would be suppressed by standard anti-aliasing methods for Kirchhoff operators. The proposed method succeeds in imaging ``beyond aliasing'' without generating aliasing noise because it exploits a priori knowledge on the dip-spectrum of the data. Therefore, the proposed method is not of general applicability but it successfully improves the image resolution when a priori assumptions on the dip-spectrum of the data are realistic. The imaging of a salt-dome flanks in the Gulf of Mexico has been enhanced by the application of the proposed method.

Imaging under the edges of salt bodies: Analysis of an Elf North Sea dataset (ps 710K) (src 1782K)

Depth imaging techniques still face difficulties under the edges of salt bodies. Elf encountered such a problem in a 3-D survey in the North Sea where they hoped to image a reflector that lays beneath a salt body. Based on the 3-D data, Elf created 2-D synthetic model of the salt body and surrounding subsurface. By analyzing this synthetic model using Kirchoff depth migration and raytracing techniques, it became clear that the problem is largely due to poor illumination. The portion of the reflector under the edge of the salt dome only received illumination from a limited range of offsets. In addition to further analysis of both the synthetic and real datasets, investigations of possible solutions such as weighting appropriate offsets, inversion and/or the use of common reflection angle gathers seem warranted.

Prestack Kirchhoff time migration for complex media (ps 2273K) (src 9K)

Constructing the seismic image in vertical time, as opposed to depth, eliminates the inherent ambiguity of resolving the vertical

Accurate linear interpolation in the extended split-step migration (ps 2412K) (src 11591K)

We present an implementation of the extended split-step migration method for strong lateral velocity gradients. Steep events are imaged by using the evanescent energy in the linear interpolation of different downward continued wavefields in every depth step and between reference velocities. The 3-D algorithm is a 3-D prestack operator based on common-azimuth continuation at every depth step. We show the extended split-step impulse response with a strong lateral velocity gradient of 0.5

Prestack depth migration in anisotropic heterogeneous media (ps 904K) (src 9875K)

This paper presents a 2-D and 3-D prestack depth migration in anisotropic media for P-waves. Assuming an acoustic VTI medium, the double square root (DSR) equation becomes dependent only on the migration velocity field and the parameter . In order to handle lateral velocity variation, I use the extended split-step approximation of the double square root. I tested this algorithm with two different 2-D synthetic seismic data sets in a VTI medium, and the results are encouraging. The first VTI synthetic model consisted of a set of dipping reflectors, from to with Thomsen's parameters and . The second model is the anisotropic Marmousi model characterized by strong lateral velocity and variation. In order to handle the lateral variation, I define a number of reference 's that in the same fashion that reference velocities are defined. The resulting anisotropic prestack Marmousi section correctly imaged dipping events. In addition, I show other possible implementations of this anisotropic migration, in order to handle variation as a function of depth and lateral coordinates. In the case where the anisotropic medium has non-zero , the extended split-step migration algorithm works in pseudo-depth to avoid the explicit dependency of the DSR operator on the vertical velocity.

Dip moveout (DMO) is often applied to prestack data to better preserve dipping events when performing partial stacks over ranges of offset. The tests presented in this paper, conducted on the SEG-EAGE salt data set, indicate that the application of azimuth moveout (AMO) in place of constant-velocity DMO yields better partial stacks in two important cases: first, when the velocity increases with depth. Second, when a salt body causes NMO-velocity conflicts between deeper flat reflectors and shallower dipping reflectors. AMO is less sensitive to velocity variations than DMO because it is a residual operator, and thus it has less tendency to overcorrect reflections that have been moved out with too high NMO velocity. AMO has also the potential advantage over

The azimuth moveout operator for vertically inhomogeneous media (ps 784K) (src 595K)

The azimuth moveout (AMO) operator, unlike the DMO operator, has a 3-D structure in homogeneous isotropic media, with an out-of-plane (crossline) component. In general, this component is concaved downward giving the operator an overall skewed-saddle shape. The AMO operator is typically smaller in size than conventional DMO operators. When velocity varies vertically, the operator shape changes depending on how the velocity varies. The general shape of the operator, however, remains overall saddle. In fact, for smooth velocity increases with depth, similar to those found in the Gulf of Mexico, the AMO operator does not vary much from its homogeneous counterpart. The residual AMO operator, constructed by cascading a forward homogeneous AMO operator with an inverse

Discrete Kirchhoff theory and irregular geometry (ps 563K) (src 2564K)

Discrete implementations of integral operators are essentially matrix-vector multiplications. In this paper, we study the structure of Kirchhoff-type matrices where each row corresponds to a summation surface and each column corresponds to an impulse response. Due to irregularities in seismic coverage, the columns and rows are generally badly scaled. To balance the coefficients of the matrix, we propose two formulations for normalization: a data-space formulation based on row scaling of pull (sum) operators, and a model-space normalization based on column scaling of push (spray) operators. In both approaches, the final image is normalized by a reference model that is the operator's response to an input vector with all components equal to one. We apply this normalization technique to approximate a data covariance matrix based on the definition of a data-space pseudo inverse. This data covariance is an AMO matrix. It represents an equalization filter that corrects the imaging operator for the interdependencies among data parameters. We investigate the use of the normalization operators as preconditioners for the iterative solution of multichannel inversion. The diagonal transformation ensures common magnitude to all the variables and accelerates the convergence of the linear solver. Beyond the goal of fold normalization, the advantage of the iterative solution is to interpolate for missing data and reduce the artifacts related to data aliasing.

Imaging the very near surface with a high-resolution 3-D seismic survey (ps 784K) (src 707K)

We designed and completed a very high-resolution 3-D seismic survey at Moss Landings, California. The survey was part of an ongoing series of experiments being conducted by one of the authors (Ran Bachrach of the Stanford Rock Physics and Borehole group, SRB) exploring the resolution limits of shallow/environmental reflection seismology particularly with regards to detection and monitoring of viscous fluid contaminants. This project is also part of an on-going collaboration between SRB and SEP. Figure 1 shows a 2-D shot gather in the study ...

Standard depth tomography methods often do not converge, converge slowly, or converge to a model that is geologically unrealistic. By using dip based steering filters rather than standard isotropic regularization operators, convergence speed and final model quality are improved. By formulating the tomography operator in two-way vertical traveltime (,) rather than in depth (

Interval velocity estimation with a null-space (ps 930K) (src 26621K)

Many powerful velocity estimation programs have been written based on intuition and exhaustive experience. It is not easy to improve such programs. Recently we presented a more formal approach showing how directional interpolation in interval velocity can exploit null-spaces in order to build more ``geologically pleasing'' velocity models Clapp et al. (1997). ...

Velocity continuation by spectral methods (ps 1873K) (src 4323K)

I apply Fourier and Chebyshev spectral methods to derive accurate and efficient algorithms for velocity continuation. As expected, the accuracy of the spectral methods is noticeably superior to that of the finite-difference approach. Both methods apply a transformation of the time axis to squared time. The Chebyshev method is slightly less efficient than the Fourier method, but has less problems with the time transformation and also handles accurately the non-periodic boundary conditions.

Decreased CMP fold, such as that found in multi-source, multi-streamer acquisition geometries, can hinder processing steps which benefit from well sampled CMP gathers, such as radon transforms. In two steps of linear least squares, multiscale prediction-error filters (PEF)s can estimate local dips from the recorded data and then use the dip information to fill in unrecorded shot or receiver gathers. In this paper I use multiscale, volumetric PEFs to interpolate sparse, multiple-contaminated data. The increase in fold improves multiple-suppression results.

Decon and interpolation with nonstationary filters (ps 542K) (src 366K)

Building on the notions of time-variable filtering and the helix coordinate system, we develop software for filters that are smoothly variable in multiple dimensions. Multiscale prediction-error filters (PEFs) can estimate dips from recorded data and use the dip information to fill in unrecorded shot or receiver gathers. The data are typically divided into patches with approximately constant dips. Instead, we estimate a set of smoothly varying filters, up to one PEF for every data sample. They are more memory-intensive to estimate, but the smoothly varying filters do give more accurate interpolation results than discrete patches. Finally, we offer an improved method of controlling the smoothness of the filters. We design filters like directional derivatives that we call ``flag filters''. They destroy dips in easily adjusted directions. We use them in residual space to encourage dips in the specified directions. We develop the notion of ``radial-flag filters'', i.e., flag filters oriented in the radial direction (lines of constant

Horizon refinement by synthesis of seismic and well log data (ps 603K) (src 234K)

Subsurface rock strata are often conceptualized in three dimensions as depth surfaces, or

Madagascar revisited: A missing data problem (ps 1019K) (src 10504K)

The Madagascar satellite data set provides images of a spreading ridge off the coast of Madagascar. This data set has two regions: the southern half is densely sampled and the northern half is sparsely sampled. The sparsely sampled region presents a missing data problem. I am using prediction-error filters estimated on the dense southern half to fill data on the sparse half. The prediction-error filters effectively spread the texture of the spreading ridge to the sparse tracks.

Prestack Kirchhoff time migration, for transversely isotropic media with a vertical symmetry axis (VTI media), is implemented using an offset-midpoint traveltime equation; Cheop's pyramid equation for VTI media. The derivation of such an equation required approximations that pertain to high frequency and weak anisotropy. Yet, the resultant offset-midpoint traveltime equation for VTI media is highly accurate for even strong anisotropy. It is also strictly dependent on two parameters: the normal-moveout (NMO) velocity and the anisotropy parameter, . It reduces to the exact offset-midpoint traveltime equation for isotropic media when .In vertically inhomogeneous media, the NMO velocity and parameters in offset-midpoint traveltime equation are replaced by their effective values; the velocity is replaced by the root-mean-squared velocity, whereas is given by a more complicated equation that includes summation of the fourth power of velocity.

Fast-marching eikonal solver in the tetragonal coordinates (ps 566K) (src 1555K)

Accurate and efficient traveltime calculation is an important topic in seismic imaging. We present a fast-marching eikonal solver in the tetragonal coordinates (3-D) and trigonal coordinates (2-D),

The fast marching method in Spherical coordinates: SEG/EAGE salt-dome model (ps 1345K) (src 10K)

Applying the fast marching method to solve the eikonal equation on the 3-D SEG/EAGE salt-dome model demonstrates two key features of the method, stability and efficiency. Such an application, also reveals some of the accuracy deficiencies of the Cartesian-coordinate implementation of the fast marching method. The accuracy is improved by applying the fast marching method in spherical coordinates. Obviously, this domain better represents waves emanating from a point source than the Cartesian coordinates. However, the finite-difference solution of the eikonal equation, in any domain, provides traveltimes corresponding only to the fastest arrivals. These arrivals, in inhomogeneous media, include typically head-waves and other low-energy waves. The eikonal solution of the salt-dome model includes a lot of low energy waves, such as head-waves emanating from the top of the salt structure. These low-energy waves have replaced the more important direct waves in many regions of the solution. Using such a traveltime solution for imaging will result in a less than ideal image.

Seismic reservoir monitoring is a new technology that involves interpreting differences between different generations of 3-D seismic data in terms of changes in reservoir properties over time. A vast number of 3-D seismic surveys are repeated, in whole or in part, for reasons other than reservoir monitoring. In this paper, we cross-equalize two off-the-shelf migrated datacubes over a Gulf of Mexico field with favorable reservoir properties, to see whether their time-lapse nature can be exploited for reservoir monitoring purposes. Key elements in the processing flow include spatial realignment to a common grid, wavelet equalization by

AVO inversion in

We implement a new Kirchhoff-typed AVO inversion scheme in

Rocks as poroelastic composites (ps 536K) (src 43K)

In Biot's theory of poroelasticity, elastic materials contain connected voids or pores and these pores may be filled with fluids under pressure. The fluid pressure then couples to the mechanical effects of stress or strain applied externally to the solid matrix. Eshelby's formula (for the response of a single ellipsoidal elastic inclusion in an elastic whole space to a strain imposed at infinity) is a very well-known and important result in elasticity. The hardest technical part of Eshelby's work was in computing the elliptic integrals needed to evaluate the fourth-rank tensors for inclusions shaped like spheres, oblate and prolate spheroids, needles and disks. Having a rigorous generalization of Eshelby's results valid for poroelasticity means that the hard part of Eshelby's work can be carried over from elasticity to poroelasticity - and also thermoelasticity - with only trivial modifications. Effective medium theories for poroelastic composites such as rocks can then be formulated easily by analogy to well-established methods used for elastic composites. An identity analogous to Eshelby's classic result has been previously derived by the author for use in these more complex and more realistic problems in rock mechanics analysis. Using these results as the starting point for new methods of estimation, I apply these techniques to the Biot-Willis parameter, which is the technical name for the effective-stress coefficient for total volume strain. The results show that poroelastic parameters can now be estimated as easily as elastic parameters for arbitrary ellipsoidal inclusions using any of the standard effective medium theories.

Wind a wire onto a cylinder to create a helix. I show that a filter on the 1-D space of the wire mimics a 2-D filter on the cylindrical surface. Thus 2-D convolution can be done with a 1-D convolution program. I show some examples of 2-D recursive filtering (also called 2-D deconvolution or 2-D polynomial division). In 2-D as in 1-D, the computational advantage of recursive filters is the speed with which they propagate information over long distances. We can estimate 2-D prediction-error filters (PEFs), that are assured of being stable for 2-D recursion. Such 2-D and 3-D recursions are general-purpose preconditioners that vastly speed the solution of a wide class of geophysical estimation problems. The helix transformation also enables us the partial-differential equation of wave extrapolation as though it was an ordinary-differential equation.

Factorization of cross spectra (ps 466K) (src 5K)

A solved problem is the factorization of a positive real autospectrum into a minimum-phase wavelet and its adjoint. The most practical method is that of Kolmogoroff. Here I extend the Kolmogoroff method to cross-spectra. This problem arises in the extrapolation of 3-D wavefields where we need to factor an operator like .We get a band matrix to solve. In principle, we factor it into lower and upper triangular band matrices ...

Wilson-Burg spectral factorization with application to helix filtering (ps 465K) (src 19K)

Spectral factorization methods are used for the estimation of minimum - phase time series from a given power spectrum. We present an efficient technique for spectral factorization, based on Newton's method. We show how to apply the method to the factorization of both auto and cross-spectra, and present a simple example of 2-D deconvolution in the helical coordinate system.

Helical factorization of the Helmholtz equation (ps 673K) (src 3441K)

The accuracy of conventional explicit wavefield extrapolation algorithms at high dips is directly related to the length of the convolution filters: increasing the dip range leads to increased cost. Recursive filters have the advantage over convolutional filters in that short filters can move energy long distances. We discard both Crank-Nicolson and McClellan transforms, and extrapolate waves by factoring the 3-D Helmholtz equation in a helical coordinate system. We show that one of the minimum-phase factors provides a 90extrapolator, that can be applied recursively in the (domain. By developing a purely recursive wavefield extrapolator, we hope to achieve accuracy at high dips with shorter filters than is possible with explicit methods.

A simplified interface to SEPlib is implemented with a Fortran-90 module.

Genkir3D: A toolkit for Kirchhoff imaging (ps 336K) (src 631K)

Genkir3D is a software package that facilitates the implementation of new 3-D Kirchhoff algorithms. It enables the user to specify a new Kirchhoff operator by only writing a function that computes the kinematics and the amplitudes of the summation surfaces, in addition to the maximum data frequency for anti-aliasing. All the practical issues related to data interpolation, data filtering, geometry handling, and limited memory usage, are efficiently and flexibly taken care by the package. The input and the output are data sets in Sep3D data format, and therefore can have arbitrary trace geometry. The package is developed in Fortran 90.

7/5/1998