SEP-97 (1998)

Migration

Implicit 3-D depth migration by wavefield extrapolation with helical boundary conditions

(ps 587K) (src 453K)
Rickett J., Claerbout J., and Fomel S.
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)
Biondi B. L.
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)
Prucha M. L., Clapp R. G., and Biondi B. L.
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)
Alkhalifah T.
Constructing the seismic image in vertical time, as opposed to depth, eliminates the inherent ambiguity of resolving the vertical P-wave velocity from surface seismic data in transversely isotropic media with a vertical axis of symmetry (VTI media). By ray tracing in the space-time (x-)-domain, a traveltime map is built by interpolating the traveltime information along the rays onto a regular grid in space and time. This traveltime map is used by the prestack Kirchhoff time migration to obtain the migration summation trajectories. Since the traveltime map is extracted using ray tracing, the migration can practically handle any lateral velocity variations. Specifically, the prestack time migration yields good images of the isotropic and anisotropic Marmousi models.

Accurate linear interpolation in the extended split-step migration

(ps 2412K) (src 11591K)
Malcotti H. and Biondi B.
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 s^-1, and we apply the migration algorithm over two synthetic data sets: a 2-D salt dome model and the 3-D SEG-EAEG data set. In addition, we show the resulting prestack migrated image of the real data set that used to built the 2-D salt dome model. Our results show that the evanescent energy is important in the linear interpolation to preserve steep events in a strong lateral velocity gradient.

Prestack depth migration in anisotropic heterogeneous media

(ps 904K) (src 9875K)
Malcotti H.
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.

Prestack Imaging

Azimuth moveout vs. dip moveout in inhomogeneous media

(ps 1495K) (src 2358K)
Biondi B.
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 V(z) DMOs to be less sensitive to the given velocity function.

The azimuth moveout operator for vertically inhomogeneous media

(ps 784K) (src 595K)
Alkhalifah T. and Biondi B. L.
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 v(z) one is extremely small, which suggests that the impact of such v(z) variations on the AMO operator is generally small. Complex vertical velocity variations, on the other hand, result in more complicated AMO operators that include, among other things, triplications at moderate angles. Regardless of the complexity of the model, the v(z) operator has the same first order behavior as its homogeneous counterpart. As a result, for small dip angles the homogeneous AMO, as a tool for partial stacking, often enhances the image. Moderate to steep dips in complex v(z) media requires the application of an algorithm that honors such velocity variations.

Discrete Kirchhoff theory and irregular geometry

(ps 563K) (src 2564K)
Chemingui N. and Biondi B.
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)
Rickett J. and Bachrach R.
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 ...

Velocity Estimation

Regularizing time tomography with steering filters

(ps 632K) (src 659K)
Clapp R. G. and Biondi B. L.
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 (z,x) coordinates, we can uncouple velocity estimation from reflector position estimation. In this coordinate system we avoid many of the problems associated with the depth tomography problem, and converge quickly to a reasonable solution.

Interval velocity estimation with a null-space

(ps 930K) (src 26621K)
Clapp R. G., Sava P., and Claerbout J. F.
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)
Fomel S.
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.

Applied Inversion

Shot interpolation for radon multiple suppression

(ps 2566K) (src 7399K)
Crawley S.
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)
Crawley S., Clapp R., and Claerbout J.
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 x/t in (t,x) space). Break a common-midpoint gather into pie shaped regions bounded by various values of x/t. Such a pie-shaped region tends to have constant dip spectrum throughout the region so it is a natural region for smoothing estimates of dip spectra or of gathering statistics (via 2-D PEFs). In this paper we use smoothly variable PEFs to interpolate missing traces, though obviously they have many other uses.

Horizon refinement by synthesis of seismic and well log data

(ps 603K) (src 234K)
Brown M.
Subsurface rock strata are often conceptualized in three dimensions as depth surfaces, or horizons. By picking the corresponding reflection event from 3-D seismic data, an approximate ``seismic horizon'' can be obtained, but in general the inconsistencies in the approximation are nontrivial. Picks from well log data give the depth to the horizon to high accuracy, but only at a few points. To achieve the final goal, which is a representation of the actual horizon more accurate than the seismic horizon, I present a least squares optimization scheme to optimally ``warp'' the seismic data to match the well log measurements. I then test the method on a small group of data from offshore West Africa, consisting of two picked seismic horizons data from about 15 well logs. The method shows promise in more far-reaching problems such as anisotropic parameter estimation and global velocity update.

Madagascar revisited: A missing data problem

(ps 1019K) (src 10504K)
Lomask J.
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.

Traveltimes

The offset-midpoint traveltime pyramid in transversely isotropic media

(ps 1044K) (src 45K)
Alkhalifah T.
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)
Sun Y. and Fomel S.
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), tetragonal (trigonal) fast-marching eikonal solver (TFMES), which can significantly reduce the first-order approximation error without greatly increasing the computational complexity. In the trigonal coordinates, there are six equally-spaced points surrounding one specific point and the number is twelve in the tetragonal coordinates, whereas the numbers of points are four and six respectively in the Cartesian coordinates. This means that the local traveltime updating space is more densely sampled in the tetragonal ( or trigonal) coordinates, which is the main reason that TFMES is more accurate than its counterpart in the Cartesian coordinates. Compared with the fast-marching eikonal solver in the polar coordinates, TFMES is more convenient since it needs only to transform the velocity model from the Cartesian to the tetragonal coordinates for one time. Potentially, TFMES can handle the complex velocity model better than the polar fast-marching solver. We also show that TFMES can be completely derived from Fermat's principle. This variational formulation implies that the fast-marching method can be extended for traveltime computation on other nonorthogonal or unstructured grids.

The fast marching method in Spherical coordinates: SEG/EAGE salt-dome model

(ps 1345K) (src 10K)
Alkhalifah T.
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.

Reservoir Characterization

A cross-equalization processing flow for off-the-shelf 4-D seismic data

(ps 1920K) (src 4521K)
Rickett J. and Lumley D.
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 L₂ matched-filtering, warping to compensate for different stacking/migration velocities, and amplitude balancing. A balance is struck between global and local operators to ensure a clean cross-equalization result without inadvertently removing changes due to fluid production.

AVO inversion in V(x,z) media

(ps 3540K) (src 5355K)
Sun Y. and Dong W.
We implement a new Kirchhoff-typed AVO inversion scheme in V(x,z) media. The WKBJ Green's function is calculated using a finite-difference scheme. We propose a pair of Kirchhoff inversion operators which have more obvious physical meaning. By analyzing the Kirchhoff inversion operator, we find out an unique relationship between the weighting function and the kinematic equation, which is very important to recover the true amplitude of the reflection coefficient. Our scheme is a two-step AVO inversion approach. Common-image gathers (CIG) are generated in the first step. These common-image gathers can be used to update the velocity model and reduce the influence of velocity error in the final AVO inversion results. AVO intercepts and slopes are estimated in the second step using a least-squares procedure. Finally, a fluid-line section is generated to highlight the existence of V_p/V_s anomaly. One dipping-layered synthetic example demonstrates the accuracy of our scheme and the influence of NMO stretch on the estimated AVO coefficients. The result from a field data example, the Mobil AVO dataset, shows a strong V_p/V_s anomaly in the middle of the section that may be a potential hydrocarbon indicator.

Rocks as poroelastic composites

(ps 536K) (src 43K)
Berryman J. G.
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.

Helix Theory

Multidimensional recursive filters via a helix

(ps 2249K) (src 1780K)
Claerbout J.
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)
Claerbout J.
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)
Sava P., Rickett J., Fomel S., and Claerbout J.
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)
Rickett J. and Claerbout J.
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.

Software Modules

"SEP" module: A Fortran-90 interface to SEPlib

(ps 349K) (src 10K)
Fomel S.
A simplified interface to SEPlib is implemented with a Fortran-90 module.

Genkir3D: A toolkit for Kirchhoff imaging

(ps 336K) (src 631K)
Biondi B.
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.

References

• Claerbout, J. F., 1997, Geophysical exploration mapping: Environmental soundings image enhancement: Stanford Exploration Project.
• Clapp, R. G., Fomel, S., and Claerbout, J., 1997, Solution steering with space-variant filters: SEP-95, 27-42.
• Fomel, S., Clapp, R., and Claerbout, J., 1997, Missing data interpolation by recursive filter preconditioning: SEP-95, 15-25.
• Nash, S. G., and Sofer, A., 1996, Linear and nonlinear programing: McGraw-Hill.