A finite element code on the Connection Machine (ps 74K) (src 36K)
**Nichols D.**

I have written a finite element
code for the Connection Machine (CM) that models heat flow in a
circuit element. The code uses the conjugate gradient method to solve
the global Galerkin matrix associated with the finite element problem.
The design of the algorithm depends on the
connectivity of the elements. For regularly connected elements, the
application of the global matrix can be coded as a banded matrix
multiply. This is very efficient on the CM. If a regular mesh is
used the solution of the finite element system runs at 275MFlops on 4K
CM. For irregularly connected elements the communication
is more complicated and the global data transfer is much more difficult to
implement efficiently on the CM. The basic design considerations
that apply to this simple problem should be the same for more
complicated problems such as solution of the wave equation.

Porting a simple 2-D migration program to the Connection Machine (ps 65K) (src 14K)
**Cole S.**

As a learning exercise, I ported a simple 2-D implicit finite-difference
poststack migration program to the Connection Machine.
Getting the program up and running was not difficult,
but figuring out how to get it to perform well was something of a challenge.
Currently the program runs roughly 25 times faster than the
equivalent well-vectorized program running on a Convex C-1.
Since the algorithm itself is very familiar to most readers,
those with no prior Connection Machine experience
might find this review of my experiences a useful introduction.

Three-dimensional time-slice migration on the Connection Machine (ps 48K) (src 29K)
**Karrenbach M.**

Three-dimensional migration of zero-offset data using a depth-variable
velocity can be performed in one pass
using Fourier transforms of time slices. The migration process is
carried out entirely in the two-dimensional spatial Fourier domain.
Three-dimensional time-slice migration is inherently a parallel algorithm
and maps perfectly onto the natural layout of processors on the Connection
Machine.
This method is especially suitable for target-oriented migration of a
subset of a three-dimensional data volume.
...

3-D depth migration by a predictor-corrector method (ps 67K) (src 561K)
**Nichols D.**

3-D depth migration can be cast as an explicit or implicit finite
difference problem. The explicit problem is cheap to solve, but
unstable. The stable implicit problem can be solved by direct or
iterative methods. Direct methods are expensive and iterative methods
are slow to converge at low frequencies. I experiment with the use of
a predictor-corrector method, using iterative refinement of an
approximate solution to the implicit problem. The approximate solution
is obtained by solving the explicit problem. The iterative, conjugate
gradient method is still slow to converge at low frequencies. A
preconditioned conjugate gradient method would converge faster.

Singular value decomposition on the Connection Machine (ps 69K) (src 10K)
**Rub L.**

The singular value decomposition of a matrix provides us with valuable information about the structure of the matrix. The matrices arising in geophysical data analysis are very large and consequently the resulting singular value decompositions are computationally expensive to obtain. The following implementation on a Connection Machine shows that with appropriate computer hardware and programming, matrix size will not be as much of an issue.

Experience with CM Fortran (ps 62K) (src 639K)
**Muir F.**

I had originally intended to call this Short Note ``Experience with the
Connection Machine'', but on reflection it is the Fortran interface that I
know, not the machine itself, and this is how it should be. I am in no
sense a system programmer, but I have many years experience developing
algorithms and translating them into Fortran subroutines, so it is with this
background that I discuss my chosen topic.
...

Wave-equation algorithms on massively parallel computers (ps 126K) (src 158K)
**Biondi B.**

The wave-equation algorithms used in reflection seismology
are based on the combination of the following three methods:
Fourier domain methods, finite differencing, and Kirchhoff
integrals.
These basic methods can be efficiently implemented on
a massively parallel computer using simple algorithms.
The proposed algorithms are easily expressed
in the new Fortran 90 language.
I implemented and ran the Fortran 90 codes on the SEP CM
and measured good performances.

Residual depth migration (ps 266K) (src 463K)
**Zhang L.**

Residual depth migration is a process of converting an image migrated
with one slowness model to an image migrated with another slowness model.
The kinematic operators of this process are defined implicitly
in a pair of nonlinear equations. This nonlinear system can be efficiently
solved by first transforming it into the initial value problem of a pair
of partial differential equations, and then solving the partial differential
equations with a finite-difference scheme. Once the operators are calculated,
the residual depth migration can be accomplished with Kirchhoff-type integrals.
The proposed algorithm handles general variable-slowness models, and it is
applicable to both pre- and post-stack images. A synthetic example of residual
depth migration of a post-stack image shows that the algorithm
gives kinematically accurate results.

Migration to zero-offset in variable velocity media (ps 281K) (src 256K)
**Popovici A. M.**

Migration to Zero-Offset (MZO) or Prestack Partial Migration (PSPM)
transforms prestack data into zero-offset data and is
equivalent in a constant velocity medium
to normal moveout correction (NMO)
followed by the dip moveout correction (DMO) applied
as a single step process.
I present an algorithm to construct the MZO operator in
variable velocity media as a
potentially efficient alternative to full prestack migration.
The MZO impulse response is computed using finite-difference
travel-time maps, by considering the MZO process
as the combination of
full prestack migration followed by
zero-offset modeling.
The generalized impulse response can
be applied to the data as an integral operator.
I apply the operator to different synthetic models with
depth variable velocity and show the improvements
in stacking alignment over conventional DMO and NMO processing.

Depth migration by an unconditionally stable explicit finite-difference method (ps 175K) (src 1187K)
**Ji J. and Biondi B.**

We develop an unconditionally stable, explicit depth migration method.
The downward continuation operator derived by a finite-difference
approximation of the one-way equation is given by the exponential
of a banded matrix.
We approximate this exponential by decomposing the banded matrix
into block diagonal matrices, of which the
exponential can be computed analytically.
The derived downward continuation operator
is explicit and unconditionally stable,
and thus it may be efficiently implemented on either
vector or parallel computers.
First, we apply this algorithm to a 15-degree migration.
To increase the accuracy at steep dip, we add
more terms to the Taylor-series expansion of the square-root operator.
However, the improvement in accuracy gained adding
more terms to the Taylor series,
it is at the expense of higher computational cost.

Uniform-amplitude DMO (ps 249K) (src 550K)
**Zhang L.**

In SEP-59 (Zhang, 1988), I describe a new derivation of the
dip-moveout (DMO) operator
in the Fourier domain. The new operator is kinematically
equivalent to Hale's DMO operator but has a different amplitude spectrum.
From the results of synthetic experiments, I noticed that the
impulse responses of the operator have properties
that agree qualitatively with wave theory: the operator tends to give
a uniform amplitude response for reflectors of varying dips.
However, in that study I did not attempt to make a quantitative explanation of
such an amplitude phenomenon.
...

Anisotropic scalar imaging using the double elliptic approximation (ps 322K) (src 975K)
**Karrenbach M.**

The double elliptic approximation to the transverse isotropic dispersion
relation can be used to image a subsurface scalar eigenfield.
The rationale behind the
double elliptic approximation is the use of only four parameters for 2D
data. These parameters are found by conventional velocity analysis
after a vector wavefield is converted to its scalar eigenfields.
The dispersion relation is parameterized by
horizontal and vertical direct-wave velocities and
horizontal and vertical normal-moveout velocities.
Estimation of those four parameters
typically requires a combination of
surface seismic data and cross-well or VSP data.
If the algorithm is implemented in the phase domain,
existing scalar imaging techniques are easily extended to anisotropy
As an example a transverse
isotropic medium (Greenhorn shale) is modeled and migrated using both
the exact anisotropic dispersion relation and its double elliptic
approximation. The approximation works very well paraxially, but it fails
in areas of triplication.

Prestack reverse-time migration in anisotropic media (ps 1082K) (src 6375K)
**Karrenbach M.**

The finite-difference method used with the full two-way wave equation
can model all wave phenomena exhaustively and provides correct amplitude
information at reflectors.
Reverse time migration,
or conjugate-gradient inversion, uses conjugate modeling to arrive at an
approximate inverse to the forward modeling process.
Reverse time migration is expensive compared to one-way wave
equations or algorithms with simplified reflection geometry, but it
handles all amplitudes properly and includes all wave propagation phenomena.
Reverse time migration is suggested as the preferred tool for detailed
medium property identification.I show a two-dimensional example, where
elastic siffnesses in an anisotropic
media can be imaged and determined by crosscorrelating strain components of
a forward-modeled shot wave field with strain components of the
backward-propagated data wavefield.

Simplifying 3-D migration operators using Givens rotations (ps 56K) (src 141K)
**Cole S.**

The pentadiagonal form of the 2-D Laplacian makes 3-D migration
schemes inaccurate (with splitting) or expensive. Using Givens rotations,
and assuming that velocity is constant or varies
only with depth, it is easy and inexpensive to
reduce the pentadiagonal system to a tridiagonal form,
since the algorithm lends itself well to a parallel implementation.
Unfortunately, for the case of lateral velocity variation, the matrix
rapidly becomes less sparse as rotations are performed, making
the method too expensive. But, even in the presence of lateral
velocity variation, this scheme may be a useful preconditioner for
other solution methods.

Some aspirations in seismic reservoir characterization (ps 42K) (src 5K)
**Lumley D. E.**

Having recently arrived as a Ph.D. student of the SEP as of a few days ago,
I will now take the liberty of
introducing myself by discussing some
problems which interest me and some thoughts I have regarding their solution.
My interests reside generally in the realm of migration/inversion (m/i)
of seismic reflection data, with a particular emphasis on subsurface target
characterization.
By m/i, I mean hybridizing separate but related
concepts from migration and inversion camps, in order to produce
...

Elastic parameter estimation by Kirchhoff prestack depth migration/inversion (ps 540K) (src 481K)
**Lumley D. E. and Beydoun W. B.**

Obtaining multiparameter elastic depth images is a key step toward
seismic reservoir characterization. We develop a prestack depth
migration/inversion (m/i) method which produces such images in two steps.
The first step involves a ``true amplitude''
Kirchhoff elastic prestack depth migration, which provides migrated
depth images of the elastic reflectivity *R*_{pp}, and
the associated specular reflection angles .
The second step involves fitting three elastic
parameters to the migrated elastic specular reflectivity by use of the
Bortfeld linearized approximation to the nonlinear Zoeppritz
relation, which results in three migrated elastic parameter depth images
of the target zone. We investigate the m/i in terms of
stability and choice of parameterization, and calculate quantitative
``confidence'' images to appraise our results. We find that relative
changes in P and S impedances are the most robust parameters for standard
surface seismic reflection acquisition geometries, and that little or no
unambiguous information regarding density variation is contained in the
standard specular illumination range of 0-30. Our method
is tested with synthetic and field data examples, the latter to be
presented at a later date pending proprietary release.

Model decomposition with Walsh functions in nonlinear traveltime inversion (ps 124K) (src 573K)
**Filho C. A. C.**

In Cunha (1990) I used a set of *sine*-like and *cosine*-like square
functions to describe the velocity model in a
traveltime inversion scheme for elliptically anisotropic media.
This specific decomposition of the model allowed the implementation
of a fast nonlinear algorithm that retrieves a different
frequency-component of the model at each step. A drawback of this
decomposition is that the basic set is neither orthogonal nor
complete. Although the non-orthogonality can be somewhat useful
(in this nonlinear scheme) for correcting errors in the
components estimated by previous steps, the lack of completeness
...

Anisotropic tomography (ps 250K) (src 254K)
**Michelena R. J. and Muir F.**

Tomographic estimation of velocities is usually performed by fitting
picked travel-time data to a set of time/distance equations using a slowness
function which is not dependent on angle. If the medium is assumed to be
transversely isotropic with vertical symmetry axis, an elliptic form can be
used to approximate the traveltime-distance relationship. When this
relationship is used to fit the data, the result is two images representing
the vertical and horizontal components of slowness. The inversion is a simple
extension of the well known isotropic schemes and whether it is successful
depends on the range of ray angles available. This is illustrated with
synthetic and field-data examples.

Interval velocity estimation with event picking (ps 761K) (src 6628K)
**Filho C. A. C.**

Most interval-velocity estimation methods are formulated
as an iterative optimization scheme which requires an initial model to
start the inversion. In addition, ray-tracing or approximate
equations are used to relate the model to the objective function.
For the horizontally-layered earth, it is possible to directly
relate the interval velocity between two interfaces to the
difference in shapes of the two reflected wavefronts corresponding
to those interfaces. I have devised a velocity estimation
method based on this relation,
using beam-stacks to define a probability
density function for the horizontal slowness, and an automatic
picking algorithm to define the reflection events. The results
show that good resolution is achieved when the method is applied
to synthetic and real data.

Finite-difference traveltime maps (ps 344K) (src 298K)
**Popovici A. M.**

The finite-difference traveltime algorithm
presented by Van Trier and Symes (1989; 1990)
uses a first-order upwind finite-difference
scheme described by Engquist and Osher (1980).
The Engquist-Osher scheme has an
underlying physical minimization condition important
for maintaining the stability of the algorithm;
it also introduces approximations in
the estimation of the traveltime field in
order to obtain greater computational speed at
the expense of accuracy.
The error can be reduced in
a slower algorithm.
I extend the algorithm in 3-D and show
results for a constant velocity medium and slowly
varying velocity medium.
I find there are greater stability problems in 3-D than in 2-D
for rapidly varying velocity
models.

Applying finite-difference method to ray tracing (ps 276K) (src 485K)
**Zhang L.**

For a given slowness model, the traveltime field can be calculated, on a regular
grid, with finite-difference methods. The contour lines of this field define
the trajectories of wavefronts. From the orthogonal relations between
wavefronts and rays, one can calculate another field whose contour lines
define the trajectories of rays.
Two-point ray tracing can then be accomplished by finding the contour line of a
given value. The method is limited to rays traveling in one direction and
associated with first arrivals. The algorithm can efficiently trace rays to
destinations on a regular grid, and thus can be applied to tomographic
inversion and time-to-depth conversion of migrated images

Absorbing boundary conditions for wave equation calculation (ps 247K) (src 1070K)
**Mo L. and Shan B.**

Because of the limitations of our seismic survey profile and computer core
size, we can only solve the infinite domain wave equation on a finite domain.
The computational edges we introduce generate artificial reflections. These
reflections, when propagated inward, mask the true solution of the problem.
So we must establish absorbing boundary conditions through which the waves
propagate with as little reflection as possible.
The subject of absorbing boundary conditions has been extensively studied
by many authors, e. g. , Clayton et al. (1977), Reynolds (1978) for modeling,
and Clayton et al. (1980) for migration. In this paper we present a method
...

Automatic picking and its applications (ps 249K) (src 507K)
**Zhang L.**

I have developed a new automatic picking algorithm that can be applied to many
picking problems in seismic data processing.
This algorithm does the initial picking
by solving a constrained non-linear optimization problem with a fast
search algorithm. Then, by linearizing the model of data, the objective
function defined in the optimization process is approximately reduced to a
quadratic form. The residual corrections are obtained by solving a linear
equation.
I applied this algorithm to several practical problems, including
dip-picking, event-picking, well-log interpolation and velocity picking.
Examples with field data show that the results are reasonably
accurate and reliable.

Well-to-well log correlation: choosing the matching attribute (ps 615K) (src 1476K)
**Filho C. A. C.**

The geological factor that one wishes to resolve must be the guiding factor
for choosing the correlation-attribute in well-to-well log matching.
Opting for a smoothed version of the sonic logs leads to
the matching of equivalent lithologic units, whereas using an attribute that
measures the local relative variability of the logs results in the
correlation of iso-chronological events.

A simple modeling scheme for layered anisotropic media (ps 377K) (src 1023K)
**Nichols D.**

Wave propagation in an anisotropic medium can be modeled by a two-way
phase shift scheme. The purpose of the method I present is to
calculate the particle velocities at the surface of a stack of
anisotropic plane layers generated by a traction at the surface. This
method is simple but flexible. It can be used to generate primary-only
sections or to model any order of multiples. The user can also choose
to ignore propagation of waves of a particular wavetype in any layer.
This permits, for example, the creation of output sections containing
only P-waves even though the full elastic reflection coefficients have
been modeled.

A simplified canonical form for monoclinic systems (ps 37K) (src 3K)
**Muir F.**

Anisotropy can be introduced into earth modeling in a useful and natural way
by recognizing a correspondence between anisotropic complexity and the
complexity of the fabric of an otherwise relatively undisturbed sedimentary
basin. The correspondence I have in mind looks something like this:

**isotropic**- Uniform depositional environment. Sediments are not much altered by overburden.
**transverse isotropic**- Layering effects. Variable deposition. Alteration ...

The double-elliptic approximation in the group and phase domains (ps 117K) (src 203K)

Elliptical anisotropy has found wide use as a simple approximation to transverse isotropy because of a unique symmetry property (an elliptical dispersion relation corresponds to an elliptical impulse response) and a simple relationship to standard geophysical techniques (hyperbolic moveout corresponds to elliptical wavefronts; NMO measures horizontal velocity, and time-to-depth conversion depends on vertical velocity). However, elliptical anisotropy is only useful as an approximation in certain restricted cases, such as when the underlying true anisotropy does not depart too far from ellipticity ...

Various equations for TI media (ps 44K) (src 3K)

Waves and rays are necessary components of any suite of seismic processes. In an isotropic medium the transition is seamless-not so in the simplest anisotropic material. Here, not only do velocities depend on direction, but worse, rays and waves have different velocities, and worse still, there is no local transformation from one kind of velocity to the other. In an earlier paper (Muir, 1990) I described a partial solution where approximate plane wave equations and ray equations have simple, interchangeable properties. In this note I relate the parameters of these approximants to the elastic moduli of the

...

A drill-bit source experiment using a 2-D array (ps 341K) (src 11642K)
**Cole S.**

I conducted a drill-bit source experiment using a
2-D array of 240 geophones. The directional resolving
power of the array offers the possibility to detect
drill-bit signals in the presence of strong surface noise sources.
A stacking-based scheme for locating sources in 3-D
does an excellent job of finding direct arrivals from
the drill bit.

Analysis of some aspects of locating sources with passive seismic data: synthetic examples (ps 46K) (src 95K)
**Vanyan L.**

Locating sources of seismic energy by computing the semblance
along moveout trajectories requires knowledge of the velocity
structure.
Since our knowledge of the medium velocities is always imperfect,
there will be errors in the location of sources.
Using an incorrect average velocity results in a vertical shift of
the semblance maximum as compared with the actual source location.
The presence of random travel time fluctuations
causes a decrease in the maximum semblance value without
changing its location.

Sources can be located even if the data are aliased.
Synthetic tests show that while some artifacts are introduced by
the aliasing, meaningful results can still be obtained, especially
when the source signal is broadband.
For any proposed experiment, such synthetic tests should be
run to see how well one can expect to locate sources given the
proposed array geometry.

Subtracting sources in the frequency domain (ps 42K) (src 55K)
**Vanyan L.**

When using passive seismic data to detect and locate
sources of seismic energy related to various endogenous
processes the typical situation is that these sources
have to be determined against a strong background, usually
of surface origin. Vanyan and Cole (1990) tried
to cope with this problem using a least squares
method for subtracting unwanted signals in frequency domain.
Here I give some formal basis for applying this
method to passive data, describe a more general multichannel
approach, and show the results of applying both methods to
...

Electronic book update (ps 47K) (src 8K)
**Claerbout J. F.**

I am writing a new textbook
(and revising an old one)
that will appear in electronic form as well as paper form.
On the workstation screen appear pages with pasted-in figures
and push buttons in the figure captions.
A pushbutton takes you to a directory where the figure gets made.
From there, some figures spring up into interactive programs
and others into batch programs.
Reading is simplest when remote sites use networks and X
to read my Stanford copy of the book.
I am also attempting to distribute the book on tape
but the main drawback there
is the volume of complicated software required to be installed,
including cake, ratfor, seplib, vplot, TEX and L^{A}TEX, XView, and
the X window system.

Introduction to seplib and SEP utility software (ps 52K) (src 6K)
**Claerbout J.**

Most of the seismic utility software at the Stanford Exploration Project
handles seismic data as a rectangular lattice or ``cube'' of numbers.
Each cube processing program appends to the history file for the cube.
Preprocessors extend fortran
to enable it to allocate memory at run time,
to facilitate input and output of data cubes,
and to facilitate self-documenting programs.

12/18/1997