SEP -- TABLE OF CONTENTS |

When solving a missing data problem, geophysicists and geostatisticians have very similar strategies. Each use the known data to characterize the model's covariance. At SEP we often characterize the covariance through Prediction Error Filters (PEFs) (Claerbout, 1998). Geostatisticians build variograms from the known data to represent the model's covariance (Issaks and Srivastava, 1989). ...

Iterative least-square inversion for amplitude balancing [pdf 296K]

Variations in source strength and receiver amplitude can introduce a bias in the final AVO analysis of prestack seismic reflection data. In this paper we tackle the problem of the amplitude balancing of the seismic traces from a marine survey. We start with a 2-D energy map from which the global trend has been removed. In order to balance this amplitude map, we first invert for the correction coefficients using an iterative least-square algorithm. The coefficients are calculated for each shot position along the survey line, each receiver position in the recording cable, and each offset. Using these coefficients, we then correct the original amplitude map for amplitude variations in the shot, receiver, and offset directions.

Least-square inversion with inexact adjoints. Method of conjugate directions: A tutorial [pdf 288K]

This tutorial describes the classic method of conjugate directions: the generalization of the conjugate-gradient method in iterative least-square inversion. I derive the algebraic equations of the conjugate-direction method from general optimization principles. The derivation explains the ``magic'' properties of conjugate gradients. It also justifies the use of conjugate directions in cases when these properties are distorted either by computational errors or by inexact adjoint operators. The extra cost comes from storing a larger number of previous search directions in the computer memory. A simple program and two examples illustrate the method.

Asymptotic pseudounitary stacking operators [pdf 344K]

Stacking operators are widely used in seismic imaging and seismic data processing. Examples include Kirchhoff datuming, migration, offset continuation, DMO, and velocity transform. Two primary approaches exist for inverting such operators. The first approach is iterative least-squares optimization, which involves the construction of the adjoint operator. The second approach is asymptotic inversion, where an approximate inverse operator is constructed in the high-frequency asymptotics. Adjoint and asymptotic inverse operators share the same kinematic properties, but their amplitudes (weighting functions) are defined differently. This paper describes a theory for reconciling the two approaches. I introduce a pair of the

Forward interpolation [pdf 780K]

As I will illustrate in later chapters, the crucial part of data regularization problems is in the choice and implementation of the regularization operator or the corresponding preconditioning operator . The choice of the forward modeling operator is less critical. In this chapter, I discuss the nature of forward interpolation, which has been one of the traditional subjects in computational mathematics. Wolberg (1990) presents a detailed review of different conventional approaches. I discuss a simple mathematical theory of interpolation from a regular grid and derive the main formulas from a very general idea of function bases. Forward interpolation plays only a supplementary role in this dissertation, but it has many primary applications, such as trace resampling, NMO, Kirchhoff and Stolt migrations, log-stretch, and radial transform, in seismic data processing and imaging. Two simple examples appear at the end of this chapter.

Spitz makes a better assumption for the signal PEF [pdf 128K]

In real-world extraction of signal from data we are not given the needed signal prediction-error filter (PEF). Claerbout has taken , the PEF of the signal, to be that of the data, . Spitz takes it to be . Where noises are highly predictable in time or space, Spitz gets significantly better results. Theoretically, a reason is that the essential character of a PEF is contained

Seismic reflection data interpolation with differential offset and shot continuation [pdf 440K]

I propose a finite-difference offset continuation filter for interpolating seismic reflection data. The filter is constructed from the offset continuation differential equation and is applied on frequency slices in the log-stretch frequency domain. Synthetic and real data tests demonstrate that the proposed method succeeds in structurally complex situations where more simplistic approaches fail.

Applications of plane-wave destruction filters [pdf 776K]

Plane-wave destruction filters originate from a local plane-wave model for characterizing seismic data. These filters can be thought of as a - analog of - prediction-error filters and as an alternative to - prediction-error filters. The filters are constructed with the help of an implicit finite-difference scheme for the local plane-wave equation. On several synthetic and real-data examples, I demonstrate that finite-difference plane-wave destruction filters perform well in applications such as fault detection, data interpolation, and noise attenuation.

Inverse B-spline interpolation [pdf 1.3M]

B-splines provide an accurate and efficient method for interpolating regularly spaced data. In this paper, I study the applicability of B-spline interpolation in the context of the inverse interpolation method for regularizing irregular data. Numerical tests show that, in comparison with lower-order linear interpolation, B-splines lead to a faster iterative conversion in under-determined problems and a more accurate result in over-determined problems. In addition, they provide a constructive method for creating discrete regularization operators from continuous differential equations.

The two-stage missing data interpolation approach of Claerbout (1998) (henceforth, the

Multiple suppression using prediction-error filter [pdf 240K]

I present an approach to multiple suppression, that is based on the moveout between primary and multiple events in the CMP gather. After normal moveout correction, primary events will be horizontal, whereas multiple events will not be. For each NMOed CMP gather, I reorder the offset in random order. Ideally, this process has little influence on the primaries, but it destroys the shape of the multiples. In other words, after randomization of the offset order, the multiples appear as random noise. This ``man-made'' random noise can be removed using prediction-error filter (PEF). The randomization of the offset order can be regarded as a random process, so we can apply it to the CMP gather many times and get many different samples. All the samples can be arranged into a 3-D cube, which is further divided into many small subcubes. A 3-D PEF can then be estimated from each subcube and re-applied to it to remove the multiple energy. After that, all the samples are averaged back into one CMP gather, which is supposed to be free of multiple events. In order to improve the efficiency of the algorithm of estimating the PEF for each subcube, except for the first subcube which starts with a zero-valued initial guess, all the subsequent subcubes take the last estimated PEF as an initial guess. Therefore, the iteration count can be reduced to one step for all the subsequent subcubes with little loss of accuracy. Three examples demonstrate the performance of this new approach, especially in removing the near-offset multiples.

Multidimensional recursive filter preconditioning in geophysical estimation problems [pdf 1.3M]

Constraining ill-posed inverse problems often requires regularized optimization. We consider two alternative approaches to regularization. The first approach involves a column operator and an extension of the data space. It requires a regularization operator which enhances the undesirable features of the model. The second approach constructs a row operator and expands the model space. It employs a preconditioning operator, which enforces a desirable behavior, such as smoothness, of the model. In large-scale problems, when iterative optimization is incomplete, the second method is preferable, because it often leads to faster convergence. We propose a method for constructing preconditioning operators by multidimensional recursive filtering. The recursive filters are constructed by imposing helical boundary conditions. Several examples with synthetic and real data demonstrate an order of magnitude efficiency gain achieved by applying the proposed technique to data interpolation problems.

Exploring three-dimensional implicit wavefield extrapolation with the helix transform [pdf 1.1M]

Implicit extrapolation is an efficient and unconditionally stable method of wavefield continuation. Unfortunately, implicit wave extrapolation in three dimensions requires an expensive solution of a large system of linear equations. However, by mapping the computational domain into one dimension via the helix transform, we show that the matrix inversion problem can be recast in terms of an efficient recursive filtering. Apart from the boundary conditions, the solution is exact in the case of constant coefficients (that is, a laterally homogeneous velocity.) We illustrate this fact with an example of three-dimensional velocity continuation and discuss possible ways of attacking the problem of lateral variations.

The Wilson-Burg method of spectral factorization with application to helical filtering [pdf 840K]

Spectral factorization is a computational procedure for constructing minimum-phase (stable inverse) filters required for recursive inverse filtering. We present a novel method of spectral factorization. The method iteratively constructs an approximation of the minimum-phase filter with the given autocorrelation by repeated forward and inverse filtering and rearranging the terms. This procedure is especially efficient in the multidimensional case, where the inverse recursive filtering is enabled by the helix transform. To exemplify a practical application of the proposed method, we consider the problem of smooth two-dimensional data regularization. Splines in tension are smooth interpolation surfaces whose behavior in unconstrained regions is controlled by the tension parameter. We show that such surfaces can be efficiently constructed with recursive filter preconditioning and introduce a family of corresponding two-dimensional minimum-phase filters. The filters are created by spectral factorization on a helix.

Spectral factorization revisited [pdf 168K]

In this paper, we review some of the iterative methods for the square root, showing that all these methods belong to the same family, for which we find a general formula. We then explain how those iterative methods for real numbers can be extended to spectral factorization of auto-correlations. The iteration based on the Newton-Raphson method is optimal from the convergence stand point, though it is not optimal as far as stability is concerned. Finally, we show that other members of the iteration family are more stable, though slightly more expensive and slower to converge.

Plane wave prediction in 3-D [pdf 772K]

The theory of plane-wave prediction in three dimensions is described by Claerbout (1999,1993). Predicting a local plane wave with T-X filters amounts to finding a pair of two-dimensional filters for two orthogonal planes in the 3-D space. Each of the filters predicts locally straight lines in the corresponding plane. The system of two 2-D filters is sufficient for predicting all but purely ...

Solution steering with space-variant filters [pdf 388K]

Most geophysical problem require some type of regularization. Unfortunately most regularization schemes produce ``smeared'' results that are often undesirable when applying other criteria (such as geologic feasibility). By forming regularization operators in terms of recursive steering filters, built from a priori information sources, we can efficiently guide the solution towards a more appealing form. The steering methodology proves effective in interpolating low frequency functions, such as velocity, but performs poorly when encountering multiple dips and high frequency data. Preliminary results using steering filters for regularization in tomography problems are encouraging.

Random lines in a plane [pdf 132K]

Locally, seismic data is a superposition of plane waves. The statistical properties of such superpositions are relevant to geophysical estimation and they are not entirely obvious. Clearly, a planar wave can be constructed from a planar distribution of point sources. ...

Texture synthesis and prediction error filtering [pdf 476K]

The spectrum of a prediction-error filter (PEF) tends toward the inverse spectrum of the data from which it is estimated. I compute 2-D PEF's from known ``training images'' and use them to synthesize similar-looking textures from random numbers via helix deconvolution. Compared to a similar technique employing Fourier transforms, the PEF-based method is generally more flexible, due to its ability to handle missing data, a fact which I illustrate with an example. Applying PEF-based texture synthesis to a stacked 2-D seismic section, I note that the residual error in the PEF estimation forms the basis for ``coherency'' analysis by highlighting discontinuities in the data, and may also serve as a measure of the quality of a given migration velocity model. Last, I relate the notion of texture synthesis to missing data interpolation and show an example.

Multi-dimensional Fourier transforms in the helical coordinate system [pdf 632K]

For every two-dimensional system with helical boundary conditions, there is an isomorphic one-dimensional system. Therefore, the one-dimensional FFT of a 2-D function wrapped on a helix is equivalent to a 2-D FFT. We show that the Fourier dual of helical boundary conditions is helical boundary conditions but with axes transposed, and we explicitly link the wavenumber vector, , in a multi-dimensional system with the wavenumber of a helical 1-D FFT, . We illustrated the concepts with an example of multi-dimensional multiple prediction.

It can be shown that for a 1-D Earth model illuminated by random plane waves from below, the cross-correlation of noise traces recorded at two points on the surface is the same as what would be recorded if one location contained a shot and the other a receiver. If this is true for real data, it could provide a way of building `pseudo-reflection seismograms' from background noise, which could then be processed and used for imaging. This conjecture is tested on synthetic data from simple 1-D and point diffractor models, and in all cases, the kinematics of observed events appear to be correct. The signal to noise ratio was found to increase as , where is the length of the time series. The number of incident plane waves does not directly affect the signal to noise ratio; however, each plane wave contributes only its own slowness to the common shot domain, so that if complete hyperbolas are to be imaged then upcoming waves must be incident from all angles.

When is anti-aliasing needed in Kirchhoff migration? [pdf 700K]

We present criteria to determine when numerical integration of seismic data will incur operator aliasing. Although there are many ways to handle operator aliasing, they add expense to the computational task. This is especially true in three dimensions. A two-dimensional Kirchhoff migration example illustrates that the image zone of interest may not always require anti-aliasing and that considerable cost may be spared by not incorporating it.

Imaging complex structures with first-arrival traveltimes [pdf 680K]

I present a layer-stripping Kirchhoff migration algorithm which is capable of obtaining accurate images of complex structures by downward continuing the data and imaging from a lower datum. I use eikonal traveltimes in a Kirchhoff datuming algorithm for the downward continuation. After downward continuation, I perform Kirchhoff migration. The method alternates steps of datuming and imaging. Because traveltimes are computed for each step, the adverse effects of caustics, headwaves, and multiple arrivals do not develop. In principal, this method only requires the same number of traveltime calculations as a standard migration. Tests on the Marmousi data set produce excellent results.

Evaluating the Stolt-stretch parameter [pdf 972K]

The Stolt migration extension to a variable velocity case describes the velocity heterogeneity with a constant parameter, which is related to the stretch transformation of the time axis. We exploit a connection between modified dispersion relations and nonhyperbolic traveltime approximations to derive an explicit expression for the stretch parameter. This analytical expression allows one two achieve the highest possible accuracy within the Stolt stretch approximation. Using a real data example, we demonstrate an application of the explicit Stolt stretch formula for an optimal partitioning of the migration velocity in the method of cascaded migrations.

Antialiasing of Kirchhoff operators by reciprocal parameterization [pdf 496K]

I propose a method for antialiasing Kirchhoff operators, which switches between interpolation in time and interpolation in space depending on the operator dips. The method is a generalization of Hale's technique for dip moveout antialiasing. It is applicable to a wide variety of integral operators and compares favorably with the popular temporal filtering technique. Simple synthetic examples demonstrate the performance and applicability of the proposed method.

Angle-gather time migration [pdf 440K]

Angle-gather migration creates seismic images for different reflection angles at the reflector. We formulate an angle-gather time migration algorithm and study its properties. The algorithm serves as an educational introduction to the angle gather concept. It also looks attractive as a practical alternative to conventional common-offset time migration both for velocity analysis and for AVO/AVA analysis.

A method for computing first arrival traveltimes and amplitudes in a general two-dimensional (2-D) velocity model is presented. The method is the result of merging two recently published ray tracing methods. The product is a very robust algorithm that is able to produce broadband wave phenomena, such as dispersion and wavelength dependent scattering. Its ability to produce broadband wave phenomena, is achieved by performing a wavelength-dependent smoothing of the velocity model across wavefronts. In the limit of high frequency, the method reduces to geometrical ray theory. The method is able to illuminate areas of large geometrical spreading where conventional ray tracing methods may give no arrivals. The method is tested on synthetic complex velocity models.

Traveltime computation with the linearized eikonal equation [pdf 192K]

Traveltime computation is an important part of seismic imaging algorithms. Conventional implementations of Kirchhoff migration require precomputing traveltime tables or include traveltime calculation in the innermost computational loop . The cost of traveltime computations is especially noticeable in the case of 3-D prestack imaging where the input data size increases the level of nesting in computational loops. ...

Huygens wavefront tracing: A robust alternative to conventional ray tracing [pdf 840K]

We present a method of ray tracing that is based on a system of differential equations equivalent to the eikonal equation, but formulated in the ray coordinate system. We use a first-order discretization scheme that is interpreted very simply in terms of the Huygens' principle. The method has proved to be a robust alternative to conventional ray tracing, while being faster and having a better ability to penetrate the shadow zones.

A variational formulation of the fast marching eikonal solver [pdf 540K]

I exploit the theoretical link between the eikonal equation and Fermat's principle to derive a variational interpretation of the recently developed method for fast traveltime computations. This method, known as fast marching, possesses remarkable computational properties. Based originally on the eikonal equation, it can be derived equally well from Fermat's principle. The new variational formulation has two important applications: First, the method can be extended naturally for traveltime computation on unstructured (triangulated) grids. Second, it can be generalized to handle other Hamilton-type equations through their correspondence with variational principles.

A second-order fast marching eikonal solver [pdf 444K]

The fast marching method (Sethian, 1996) is widely used for solving the eikonal equation in Cartesian coordinates. The method's principal advantages are: stability, computational efficiency, and algorithmic simplicity. Within geophysics, fast marching traveltime calculations (Popovici and Sethian, 1997) may be used for 3-D depth migration or velocity analysis. ...

Azimuth moveout (AMO) transforms 3-D prestack seismic data from one common azimuth and offset to different azimuths and offsets. AMO in the time-space domain is represented by a three-dimensional integral operator. The operator components are the summation path, the weighting function, and the aperture. To determine the summation path and the weighting function, we derive the AMO operator by cascading dip moveout (DMO) and inverse DMO for different azimuths in the time-space domain. To evaluate the aperture, we apply a geometric approach, defining AMO as the result of cascading prestack migration (inversion) and modeling. The aperture limitations provide a consistent description of AMO for small azimuth rotations (including zero) and justify the economic efficiency of the method.

Theory of differential offset continuation [pdf 364K]

I introduce a partial differential equation to describe the process of prestack reflection data transformation in the offset, midpoint, and time coordinates. The equation is proved theoretically to provide correct kinematics and amplitudes on the transformed constant-offset sections. Solving an initial-value problem with the proposed equation leads to integral and frequency-domain offset continuation operators, which reduce to the known forms of dip moveout operators in the case of continuation to zero offset.

Amplitude preservation for offset continuation: Confirmation for Kirchhoff data [pdf 132K]

Offset continuation (OC) is the operator that transforms common-offset seismic reflection data from one offset to another. Earlier papers by the first author presented a partial differential equation in midpoint and offset to achieve this transformation. The equation was derived from the kinematics of the continuation process with no reference to amplitudes. We present here a proof that the solution of the OC partial differential equation does propagate amplitude properly at all offsets, at least to the same order of accuracy as the Kirchhoff approximation. That is, the OC equation provides a solution with the correct traveltime and correct leading-order amplitude. ``Correct amplitude'' in this case means that the transformed amplitude exhibits the right geometrical spreading and reflection-surface-curvature effects for the new offset. The reflection coefficient of the original offset is preserved in this transformation. This result is more general than the earlier results in that it does not rely on the two-and-one-half dimensional assumption.

Effective AMO implementation in the log-stretch, frequency-wavenumber domain [pdf 204K]

Azimuth moveout (AMO), introduced by Biondi et al. (1998), is used as part of the styling goal (in conjunction with a derivative as a roughener) in Biondi and Vlad (2001). This paper describes the implementation of AMO for the above-stated purpose, with a historical background, proof, and discussion of pitfalls and practical steps.

...

I show Shearer's earthquake stacks over all source-receiver locations at constant offset and compare them to exploration seismic data. This electronic document simply reads the stacks and plots them.

A prospect for super resolution [pdf 120K]

Wouldn't it be great if I could take signals of 10-30 Hz bandwidth from 100 different offsets and construct a zero-offset trace with 5-100 Hz bandwidth? This would not violate Shannon's sampling theorem which theoretically allows us to have a transform from 100 signals of 20 Hz bandwidth to one signal at 2000 Hz bandwidth. ...

The double-elliptic approximation in the group and phase domains [pdf 184K]

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

Velocity continuation and the anatomy of residual prestack time migration [pdf 340K]

Velocity continuation is an imaginary continuous process of seismic image transformation in the post-migration domain. It generalizes the concepts of residual and cascaded migrations. Understanding the laws of velocity continuation is crucially important for a successful application of time migration velocity analysis. These laws predict the changes in the geometry and intensity of reflection events on migrated images with the change of the migration velocity. In this paper, I derive kinematic and dynamic laws for the case of prestack residual migration from simple geometric principles. The main theoretical result is a decomposition of prestack velocity continuation into three different components corresponding to residual normal moveout, residual dip moveout, and residual zero-offset migration. I analyze the contribution and properties of each of the three components separately. This theory forms the basis for constructing efficient finite-difference and spectral algorithms for time migration velocity analysis.

Velocity continuation by spectral methods [pdf 520K]

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.

Time migration velocity analysis by velocity continuation [pdf 992K]

Time migration velocity analysis can be performed by velocity continuation, an incremental process that transforms migrated seismic sections according to changes in the migration velocity. Velocity continuation enhances residual normal moveout correction by properly taking into account both vertical and lateral movements of events on seismic images. Finite-difference and spectral algorithms provide efficient practical implementations for velocity continuation. Synthetic and field data examples demonstrate the performance of the method and confirm theoretical expectations.

Traveltime sensitivity kernels: Banana-doughnuts or just plain bananas? [pdf 132K]

Estimating an accurate velocity function is one of the most critical steps in building an accurate seismic depth image of the subsurface. In areas with significant structural complexity, one-dimensional updating schemes become unstable, and more robust algorithms are needed. Reflection tomography both in the premigrated (Bishop et al., 1985) and postmigrated domains (Stork, 1992; Kosloff et al., 1996) bring the powerful ...

Modeling 3-D anisotropic fractal media [pdf 592K]

This paper presents stochastic descriptions of anisotropic fractal media. Second order statistics are used to represent the continuous random field as a stationary zero-mean process completely specified by its two-point covariance function. In analogy to the two-dimensional Goff and Jordan model for seafloor morphology, I present the von Karman functions as a generalization to media with exponential correlation functions. I also compute a two-state model by mapping the random field from continuous realizations to a binary field. The method can find application in modeling impedances from fractal media and in fluid flow problems.

Marine seismic data from the Blake Outer Ridge offshore Florida show strong ``bottom simulating reflections'' (BSR) associated with methane hydrate occurence in deep marine sediments. We use a detailed amplitude versus offset (AVO) analysis of these data to explore the validity of models which might explain the origin of the bottom simulating reflector. After careful preprocessing steps, we determine a BSR model which can successfully reproduce the observed AVO responses. The P- and S-velocity behavior predicted by the forward modeling is further investigated by estimating the P- and S-impedance contrasts at all subsurface positions. Our results indicate that the Blake Outer Ridge BSR is compatible with a model of methane hydrate in sediment, overlaying a layer of free methane gas-saturated sediment. The hydrate-bearing sediments seem to be characterized by a high P-wave velocity of approximately 2.5 km/s, an anomalously low S-wave velocity of approximately 0.5 km/s, and a thickness of around 190 meters. The underlaying gas-saturated sediments have a P-wave velocity of 1.6 km/s, an S-wave velocity of 1.1 km/s, and a thickness of approximately 250 meters.

SEP -- TABLE OF CONTENTS |

2016-03-17