SEP-70 (1991)

Parallel Computing

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.

Imaging and Signal Processing

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.

Inversion and Velocity Estimation

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<sub>pp</sub>, 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.

Modeling

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

Picking and Interpolation

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.

Anisotropy

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)
Dellinger J. and Muir F.
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)
Muir F.
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 TI media that lie behind the dispersion relations that I approximate.

Novel Seismic Sources

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

The Electronic Book and SEP Software

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^ATEX, 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.