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)

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.

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
...

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.

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.

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.

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 _{2}* 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

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.

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.

"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.

7/5/1998