Reflection Coefficient Sequence Research Articles

Tridiagonal Approach to the Algebraic Environment of Toeplitz Matrices, Part II: Zero and Eigenvalue Problems

The central subject of this paper is the three-term recurrence formula satisfied by the symmetric (first-kind) and antisymmetric (second-kind) polynomials relative to a given sequence of reflection coefficients, the last element of which has unit modulus. The theory of these polynomials is shown to have interesting analogies with the classical theory of orthogonal polynomials on the real line. In the case of real data, the former is equivalent to a special case of the latter (by a change of variable). The main application considered here is the problem of computing the zeros of the highest-degree symmetric polynomial, which is identified as a predictor polynomial. This problem occurs not only in some modelling techniques for digital signal processing, but can also be interpreted as the eigenvalue problem for a unitary Hessenberg matrix. Attractive solution methods are derived from the “tridiagonal approach,” based on the three-term recurrence relation.

Synthesis of acoustic matching networks by discrete space Fourier transform method

This paper describes a method used to analyze and synthesize acoustsic plane-wave matching structures between a transducer and a medium of arbitrary characteristic impedance. The analysis is applicable to multilayer structures with periodic discontinuities in characteristic impedance. It is shown that the frequency response of such structures is approximately equal to the discrete Fourier transform of the individual layer reflection coefficients. The synthesis procedure uses the inverse discrete time Fourier transform to generate a sequence of layer reflection coefficients, and applies a boundary condition to calculate the characteristic impedances. Using this method, wide bandwidth matching structures can be obtained using finite impulse response (FIR) filter topologies as a model. A numerical example is presented to show how the procedure is implemented, and a linear analysis is performed which models the layers of the synthesized structure as ideal transmission lines. The results of the analysis are in good agreement with the response predicted by the Fourier method.

Tridiagonal Approach to the Algebraic Environment of Toeplitz Matrices, Part I: Basic Results

This paper contains a thorough investigation of a family of symmetric “predictor polynomials” associated with a nonnegative-definite Toeplitz matrix. These polynomials are constructed from the classical predictors and from the values assumed by some dual predictors in a fixed point of unit modulus; the appropriate duality is induced by changing the sequence of reflection coefficients into its conjugate mirror image, within a unit modulus factor. The central theme of the paper is a well-defined three-term recurrence relation satisfied by these symmetric polynomials; it motivates the “tridiagonal” terminology. The properties of the recurrence are studied in detail; special attention is paid to the important issue of computing the recurrence coefficients from the reflection coefficients. It is shown how this three-term recurrence formula produces an efficient solution method, called the split Levinson algorithm, for the linear prediction problem.

Effects of averaging and sampling on the statistics of reflection coefficients

The reflection coefficients derived from sonic and density logs are frequently used in seismic exploration. Even though they measure the in‐situ formation slowness and density, sonic and density tools do not measure the exact, continuous formation properties but locally averaged properties sampled at discrete depth points. Furthermore, the logs are frequently reinterpolated to form a Goupillaud medium for many applications such as synthetic seismogram computation. Both the logging tools and the Goupillaud interpolation introduce averaging and sampling effects into the reflection coefficients and significantly alter the autocorrelation of the reflection coefficient sequence. Analytical formulas are derived to show how the autocorrelation is altered and to calculate how the autocorrelation depends on the averaging and sampling intervals. Essentially, these effects impose sincsquared envelopes on the power spectrum of the reflection coefficient sequence and alias high‐frequency components to low‐frequency components in the spectral domain. These findings are verified using synthetic and real examples.

On the geometric sequences of reflection coefficients

It is shown that a decaying geometric sequence of reflection coefficients corresponds to a predictor polynomial that has all of its roots on a circle centered at the origin of the complex z plane. It is shown that the radius of this circle is equal to the magnitude of the ratio of any two consecutive terms in the geometric sequence.< <ETX xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink">&gt;</ETX>

Optimal seismic deconvolution

Deconvolution is one of the most important aspects of seismic signal processing e.g. for oil exploration, etc. The objective of the deconvolution procedure is to remove the obscuring effect of the wavelet's replicas making up the seismic trace and consequently obtain an estimate of the reflection coefficient sequence. In this paper generic formulas for obtaining optimal estimates of the reflection coefficient sequence are derived. These formulas can utilize any estimation/smoothing scheme and yield optimal fixed-point, fixed-lag and fixed-interval estimates of the reflection coefficient sequence. New deconvolution algorithms are derived by use of the new formulas and Lainiotis' partitioning filter and smoother. The new algorithms are compared to previous ones with respect to computational complexity, memory requirements, performance and parallelism efficiency. Additional results concerning the algorithms' asymptotic behavior are given and validated through extensive computer simulation.

A comparison of stochastic gradient and minimum entropy deconvolution algorithms

This paper examines the connection between a stochastic gradient method for deconvolution, originally formulated in the electrical engineering literature , and minimum entropy type (ME-type) deconvolution methods, now well known in geophysics. It is shown that Gray's favoured ME-type deconvolution algorithm has a direct stochastic gradient equivalent which differs from the method examined here only in the presence of the inverse autocovariance matrix in Gray's method. The two procedures can be made more alike if the stochastic gradient algorithm is used merely for phase estimation, and is preceded by standard whitening deconvolution. In any event, the stochastic gradient approach is deficient in two respects, namely the assumption of an independent and identically distributed (iid) input sequence, and knowledge of its first and second absolute moments, assumptions which might be reasonable in the setting of data transmission problems, but not for seismic reflection coefficient sequences.

SOLUTION OF THE VSP ONE‐DIMENSIONAL INVERSE PROBLEM*

ABSTRACTThe acoustical impedance distribution of the substratum, or equivalently, the reflection coefficient sequence, is determined from VSP data. This nonlinear inverse problem is solved by a least‐squares method.As the wavelet is unknown, the impedance distribution and the Neumann boundary condition (which characterizes the excitation of the medium) are simultaneously identified.The inversion method is applied to synthetic and field VSP's; the result is satisfactory, even when strong noise corrupts the data, provided that a suitable constraint on the impedance distribution is introduced in order to ensure the stability of the inverse problem.The reliability of the inversion result in the case of field VSP, is confirmed. Some ways in which this result may be used are illustrated (calibration of the seismic surface data, multiple identification, prediction ahead of the bit).

Inverse Scattering with Noisy Data

In this paper we analyze error propagation in layer-peeling inversion methods. A bound for the error in recovering the reflection coefficient at a certain depth is given in terms of the estimated reflection coefficients. The error propagation results are then used to discuss some practical inversion algorithms that exploit available prior information on the reflection coefficient sequence.

Seismic time‐invariant convolutional model

A layered‐earth seismic model is subdivided into two subsystems. The upper subsystem can have any sequence of reflection coefficients but the lower subsystem has a sequence of reflection coefficients which are small in magnitude and have the characteristics of random white noise. It is shown that if an arbitrary wavelet is the input to the lower lithologic section, the same wavelet convolved with the white sequence of reflection coefficients will be the reflected output. That is, a white sedimentary system passes a wavelet in reflection as a linear time‐invariant filter with impulse response given by the reflection coefficients. Thus, the small white lithologic section acts as an ideal reflecting window, producing perfect primary reflections with no multiple reflections and no transmission losses. The upper subsystem produces a minimum‐delay multiple‐reflection waveform. The seismic wavelet is the convolution of the source wavelet, the absorption effect, this multiple‐reflection waveform, and the instrument effect. Therefore, the seismic trace within the time gate corresponding to the lower subsystem is given by the convolution of the seismic wavelet with the white reflection coefficients of the lower subsystem. The linear time‐invariant seismic model used in predictive deconvolution has been derived. Furthermore, it is shown that any layered subsystem which has small reflection coefficients acts as a linear time‐invariant filter. This explains why time‐invariant deconvolution filters can be used within various time gates on a seismic trace which at first appearance might look like a continually time‐varying phenomenon.

Non‐Gaussian reflectivity, entropy, and deconvolution

Standard deconvolution techniques assume that the wavelet is minimum phase but generally make no assumptions about the amplitude distribution of the primary reflection coefficient sequence. For a white reflection sequence the assumption of a Gaussian distribution means that recovery of the true phase of the wavelet is impossible; however, a non‐Gaussian distribution in theory allows recovery of the phase. It is generally recognized that primary reflection coefficients typically have a non‐Gaussian amplitude distribution. Deconvolution techniques that assume whiteness but seek to exploit the non‐Gaussianity include Wiggins’ minimum entropy deconvolution (MED), Claerbout’s parsimonious deconvolution, and Gray’s variable norm deconvolution. These methods do not assume minimum phase. The deconvolution filter is defined by the maximization of a function called the objective. I examine these and other MED‐type deconvolution techniques. Maximizing the objective by setting derivatives to zero results in most cases in a deconvolution filter which is the solution of a highly nonlinear Toeplitz matrix equation. Wiggins’ original iterative approach to the solution is suitable for some methods, while for other methods straightforward iterative perturbation approaches may be used instead. The likely effects on noise of the nonlinearities involved are demonstrated as extremely varied. When the form of an objective remains constant with iteration, the most general description of the method is likelihood ratio maximization; when the form changes, a method seeks to maximize relative entropy at each iteration. I emphasize simple and useful link between three methods and the use of M-estimators in robust statistics. In attempting to assess the accuracy of the techniques, the choice between different families of distributions for modeling the distribution of reflection coefficients is important. The results provide important insights into methods of constructing and understanding the statistical implications and behavior of a chosen nonlinearity. A new objective is introduced to illustrate this, and a few particular preferences expressed. The methods are compared with the zero‐memory nonlinear deconvolution approach of Godfrey and Rocca (1981); for their approach, two distinctly different yet statistically comparable models for reflection coefficients are seen to give surprisingly similarly shaped nonlinearities. Finally, it is shown that each MED‐type method can be viewed as the minimization of a particular configurational entropy expression, where some suitable ratio plays the role of a probability.

Identification and filtering

A theory of the discrete Riccati equation asymptotic behavior for a degenerate filtering problem corresponding to prediction. As an application, ARMA identification is shown to be determined from the asymptotic behavior of the reflection coefficients. Explicitly the time constants of the reflection coefficient sequence determine the moving average portion of the process. This work completes earlier work of the author on “invariant directions” of the Riccati equation.

Deconvolution applied to a near-bottom seismic profiler

Recent work attempting to relate the characteristics of the acoustic return from a near-bottom seismic profiler and the physical properties of sediment core samples motivated this study. The energy envelope of the seismogram does not reproduce the detailed structure of the underlying stratigraphy and in fact, can be misleading. Accordingly, we have modeled the physics of the process in terms of a convolution system and undertaken a comparison of the tradeoffs between three deconvolution techniques: (1) deterministic, (2) linear predition, and (3) homomorphic. The resolution and stability of each technique as a function of signal-to-noise ratio are illustrated for recovery of an underlying reflection coefficient sequence. Further, since the signals returned by significantly different stratigraphic columns may be strikingly similar in appearance, we have assessed the reliability of each processing scheme by pursuing two reflection coefficient sequences in parallel. The first, derived from physical properties of the core samples, contained no prominent features. The second was identical to the first with the addition of four large reflectors placed such that the energy envelopes of the two resulting seismograms were identical.

ZERO MEMORY NON‐LINEAR DECONVOLUTION*

AbstractA type of iterative deconvolution that extracts the source waveform and reflectivity from a seismogram through the use of zero memory, non‐linear estimators of reflection coefficient amplitnde is developed. Here, we present a theory for iterative deconvolution that is based upon the specification of a stochastic model describing reflectivity. The resulting parametric algorithm deconvolves the seismogram by forcing a filtered version of the seismogram to resemble an estimated reflection coefficient sequence. This latter time series is itself obtained from the filtered seismogram, and so a degree of iteration is required. Algorithms utilizing zero memory non‐linearities (ZNLs) converge to a family of processes, which we call Bussgang, of which any colored Gaussian process and any independent process are members. The direction of convergence is controlled by the choice of ZNL used in the algorithm. Synthetic and real data show that, generally, five to ten iterations are required for acceptable deconvolutions.

A sum autoregressive formula for the reflection response

The normal incidence unit impulse reflection response of a perfectly stratified medium is expressible as an autoregressive‐moving average (ARMA) model. In this representation, the autoregressive (AR) component describes the multiple patterns generated within the medium. The moving average (MA) component, on the other hand, bears a simple relation to the sequence of reflection coefficients (i.e., primaries only) of the layered structure. An alternate representation of the reflection response can be formulated in terms of a superposition of purely AR time‐varying minimum‐delay wavelets. Each successive addition of a deeper interface to the layered system gives rise to an AR wavelet whose leading term is equal to the magnitude of the primary reflection originating at this interface. We accordingly call these wavelets “generalized primaries.” The AR component of every generalized primary contains only those multiple reflections that arise from the addition of its particular interface to the layered medium.

Simultaneous Spherical Divergence Correction and Optima Deconvolution

The effects of spherical divergence on reflection seismograms destroys the stationarity assumption implicit in Weiner smoothing by distorting the amplitudes of the reflection coefficient sequence. The traditional correction for spherical divergence consists of an immediate scaling of the observed signal by an expected attenuation factor. It is this correction that violates the stationarity assumption, for it distorts the amplitudes of the observation noise process. Attempts to adapt the Weiner smoother to time-varying signals by time gating are cumbersome and suboptimal. In this paper, we show how to implement a Kalman deconvolution smoother which compensates for the effects of spherical divergence. The implementation is simple and straightforward, and estimates of the reflection coefficient sequence are optimal and corrected in amplitude. A computational example which compares our approach to an ad hoc Weiner smoother is provided for a simulated seismic signal. By means of this example, we demonstrate better performance for our approach over the Weiner smoother approach.

Maximum entropy and the relationship of the partial autocorrelation to the reflection coefficients of a layered system

A horizontally stratified half space bounded by a perfect reflector at the top gives rise to a seismogram which, when completed by the direct downgoing pulse at zero time and by symmetry about the origin for negative time, produces an autocorrelation function. If this autocorrelation is convolved with the corresponding prediction error operators of increasing length, we obtain a "gapped function," which deviates more and more from the perfect symmetry exhibited by the autocorrelation. This gapped function consists of the downgoing and upgoing waveforms at the top of each layer. The gap separates the two waveforms, and the gapwidth increases as deeper and deeper layers are reached. In particular, the width of the gap is a measure of the entropy of the seismogram at a given depth level-the deeper we go into the sub-surface, the higher the entropy of the corresponding gapped function. We explore the nature of the gapped function as it relates to the Toeplitz recursion generating the prediction error operators, and we re-derive the synthetic seismogram in terms of wave motion measured in units proportional to the square root of energy. We obtain an explicit relationship between the partial autocorrelation function on the one hand, and the reflection coefficient sequence on the other. This formulation allows us to generalize earlier results, so that we can treat the case for which the surface reflection coefficient is less than unity in magnitude. We investigate both the physical as well as the mathematical foundations of the stratified model, and the relation that this model bears to maximum entropy spectral analysis [1], [6]. In particular, we discuss Burg's fundamental result on the relationship of the maximum entropy spectrum to the reflection coefficient sequence characterizing the given subsurface model.

Iterative identification of non-invertible autoregressive moving-average systems with seismic applications

Let A be an invertible (i.e. minimum-delay) polynomial of degree m and let E be any polynomial of degree n. Then S = E A represents the transfer function of an autoregressive moving-average (ARMA) system. The numerator E may also be invertible, in which case the ARMA system S is invertible. In seismic exploration the sequence with z-transform given by S represents a model of an observable seismogram, where E is the z-transform of the reflection coefficient sequence of the unknown (and desired) deep reflecting interfaces and where A represents the (invertible) z-transform of the unknown inverse reverberation filter due to the surface layers. Because there is no physical reason for the reflection coefficient sequence to be invertible, the seismic model S = E A will in general be a non-invertible ARMA system. Given such a non-invertible sequence with z-transform S, the problem is to estimate E and A under the restriction that A is invertible. We give an iterative identification scheme based on the following algorithm. Initially we let the initial estimate E ( o) be unity. We then apply an iterative algorithm which for the k th step has inputs S and E ( k−1) and outputs A ( k) and E ( k) , where A ( k) is required to be invertible. Iteration is terminated when S and E (k) A (k) are approximately the same, and then E ( k) and A ( k) are the required estimates of E and A, respectively.

SINGLE‐CHANNEL WHITE‐NOISE ESTIMATORS FOR DECONVOLUTION

The Wiener filtering approach to deconvolution is limited by certain modeling assumptions, which may not always be valid. We develop a Kalman filtering approach to deconvolution which permits more flexible modeling assumptions than the Wiener filtering approach. Our approach is applicable to time‐varying or time‐invariant wavelets as well as to nonstationary or stationary noise processes. We develop equations herein for minimum‐variance estimates of the reflection coefficient sequence, as well as error variances associated with these estimates. Our estimators are compared with an ad hoc “prediction error filter,” which has recently been reported on in the geophysics literature. We show that our estimators perform better than the prediction error filter. Simulation results are included, for both time‐invariant and time‐varying situations, which support our theoretical developments.

White-noise estimators for seismic data processing in oil exploration

This paper is motivated by a problem from seismic data processing in oil exploration. We develop a Kalman filtering approach to obtaining optimal smoothed estimates of the so-called reflection coefficient sequence. This sequence contains important information about subsurface geometry. Our problem is shown to be equivalent to that of estimating white-plant noise for a linear dynamic system. By means of the equations which are derived herein, it is possible to compute fixed-interval, fixed-point, or fixed-lag optimal smoothed estimates of the reflection coefficient sequence, as well as respective error covariance-matrix information. Our optimal estimators are compared with an ad hoc error filter, (PEF) which has recently been reported on in the geophysics literature. We show that one of our estimators performs at least as well as, and in most cases, better than the prediction error filter. Simulation results are given which support and illustrate the theoretical developments.