THE RECOVERY OF SUBSURFACE REFLECTIVITY AND IMPEDANCE STRUCTURE FROM REFLECTION SEISMOGRAMS BY TIM ELLIS SCHEUER B.Sc, The University Of Utah, 1979 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF GEOPHYSICS AND ASTRONOMY We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1981 (c) Tim Ellis Scheuer, 1981 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Geoplysi'cs and A&Wonomy The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date CU. I<j . ?/ DF-6 (2/79) i i Abstract This thesis is concerned with the problem of estimating broadband acoustic impedance from normal incidence reflection se i sinograms. This topic is covered by following the linear inverse formalisms described by Parker (1977) and Oldenburg (1980). The measured seismogram is modelled as a convolution of subsurface reflectivity with a source wavelet. Then an appraisal of the seismogram is performed to obtain unique bandlimited reflectivity information. This bandlimited reflecitivity information is then utilized in two different construction algorithms which provide a broadband estimate of reflectivity; from which a broadband impedance function may be computed. The first construction method is a maximum entropy method which uses an autoregressive representation of a small portion of the reflectivity spectrum to predict spectral values outside that small portion. The second and most versatile construction method is the linear programming approach of Levy and Fullagar (1981) which utilizes the unique bandlimited spectral information obtained from an appraisal and provides a broadband reflectivity function which has a minimum 1( norm. Both methods have been tested on synthetic and real seismic data and have shown good success at recovering interpretable broadband impedance models. Errors in the data and the uniqueness of constructed reflectivity models play important roles in estimating the impedance function and in assessing its uniqueness. The Karhunen-Loeve transformation is discussed and applied on real data to stabilize the construction results in the presence of noise. The generally accepted idea that low frequency impedance information must be supplied from well log or velocity analyses because of the bandlimited nature of seismic data has been challenged. When accurate, bandlimited reflectivity information can be recovered from the seismic trace, then an interpretable, broadband impedance model may be recovered using the two construction algorithms presented in this thesis. iv TABLE OF CONTENTS Page Abstract ii List Of Illustrations , vAcknowledgements ix Background And Introduction 1 CHAPTER 1 Reflectivity And Impedance 9 CHAPTER 2 The Convolutional Model 15 CHAPTER 3 Inverse Theory And Appraisal Deconvolution .. 21 Time Domain Appraisal Deconvolution 21 Frequency Domain Appraisal Deconvolution . . . .- 27 Impedance Computation And Appraisal 29 An Example Of Appraisal Deconvolution 32 CHAPTER 4 Reflectivity Construction 37 AR Extension Of The Reflectivity Spectrum 40 Broadband Reflectivity Construction Using The L Norm And A Linear Programming Algorithm 50 General LP Forumulation 52 The Constrained LP Solution 55 Data Errors And The LP Solution 60 Inequality Constraints 6Equality Constraints 1 Stability And Efficiency Considerations 67 CHAPTER 5 Multitrace Reflectivity Construction Using Real Seismic Data 78 V The Real Data 78 The Karhunen-Loeve Transformation 87 Conclusion 99 References 104 Appendix A: Derivation Of The Common Trace 106 vi List of Illustrations Page Fig. 1 Seismic reflection geometry 3 Fig. 2 Common midpoint source-receiver pairs resulting from multifold coverage 4 Fig. 3 Processing sequence on a CDP gather 5 Fig. 4 Normal incidence reflection from a compositional boundary 9 Fig. 5 Subsurface plane layered model 11 Fig. 6 Error in the linear approximation 13 Fig. 7 The convolutional model of reflection seismograms 18 Fig. 8 a) the frequency domain representation of the convolutional model, and b) deconvolution by spectral division 20 Fig. 9 Synthetic models for an example, of appraisal deconvolution 34 Fig. 10 Appraisal deconvolution 5 Table 1 Relative resolution and variance measures for the appraisals shown in fig. 10 36 Fig. 11 The nonuniqueness inherent in reflectivity measurements 39 Fig. 12 Autoregressive reconstruction of a reflectivity function containing few reflectors 48 Fig. 13 Autoregressive reconstruction of a reflectivity function containing many reflectors 51 Fig. 14 Flow diagram for the constrained LP algorithm 8 Fig. 15 Deconvolved seismic section; represents subsurface structure in parts of western Alberta .63 vii Fig. 16 Zero phase wavelet (averaging function) estimated from data inside the box in figure 15 63 Fig. 17 Linear programming reconstruction of a reflectivity function containing a few reflectors 64 Fig. 18 Linear programming reconstruction of a reflectivity function containing many reflectors 66 Fig. 19 Synthetic test data. Wavelet has band range of 10-50HZ 68 Fig. 20 Performance of the AR and LP methods of reflectivity construction in the presence of additive noise 70 Fig. 21 Performance of the construction methods using an innaccurate wavelet 72 Fig. 22 Stability of reconstructions with respect to the frequency band used 74 Fig. 23 Stability of the AR method with respect to order and the stability of order with respect to the number of reflection coefficients; using accurate and innaccurate data 75 Fig. 24 The accuracy and efficiency of the LP solution with respect to the weighting exponent q, using accurate data 77 Fig. 25 Data section A; • deconvolved vibrator data along with a geological, interpretation 79 Fig. 26 Data section B; deconvolved explosion data 7Fig. 27 AR reflectivity and impedance reconstructions from section A; b) using unsmoothed, and c) using smoothed spectral estimates 82 Fig. 28 LP reflectivity and impedance reconstructions from section A; b) using unsmoothed, and c) using smoothed spectral estimates 83 VI 1 1 Fig. 29 AR reflectivity and impedance reconstructions from section B using smoothed spectral estimates 85 Fig. 30 LP reflectivity and impedance reconstructions from section B using smoothed spectral estimates 86 Fig. 31 Application of the K-L transformation to stabilize the choice of AR order. Also showing the principle of the mixing algorithm 92 Fig. 32 Application of the K-L transformation and mixing algorithm to the noisy section B 93 Fig. 33 Application of the K-L algorithm to LP reflectivity and impedance reconstructions 95 Fig. 34 The final reflectivity and impedance reconstructions from section A along with the original geologic interpretation 96 Fig. 35 Comparison of final reflectivity reconstructions with the original deconvolved section (A) 97 Fig. 36 Comparison of final impedance reconstructions with impedance averages obtained from the deconvolved data . 98 Fig. 37 The final reflectivity and impedance reconstructions from section B. The given velocity log has been inserted for comparison 100 Fig. 38 Comparison of final reflectivity reconstructions with the original deconvolved section (B) 101 Fig. 39 Comparison of final impedance reconstructions with impedance averages obtained from the deconvolved data 102 i ix Acknowledgements I acknowledge Ms. Kate Aasen for providing interesting diversions to thesis work; Ellis and Dorothy; Mr. Patience, my office mate, Kerry Stinson; the members of Geophysics House - Rob Conraads, Jim Horn, and Don Plenderleith; Prof. Tad Ulrych for being an inspiration in time series analysis and in particular I appreciate his help with the AR parts of this thesis; Shlomo Levy for leaving several of his then (1979) unpublished papers lying around the department which inspired parts of this thesis - Shlomo suggested use of the K-L transformation; Mr. Computer, Colin Walker; Ken 'Wizard' Whittall; and Peter Fullagar. Additional thanks go to Kate and Kerry along with Don King for improving the language of this thesis. I wish to acknowledge all members of this department for providing a dynamic environment in which to study. I have been fortunate in receiving financial assistance from Imperial Oil of Canada and I deeply appreciate the time and effort they have spent in providing a suitable data set. I am also grateful to MRO Associates Inc. for providing a second data set. Finally, I would like to thank and acknowledge my advisor Prof. Doug Oldenburg. Doug has provided me with endless motivation and support during the past two years. Dougs' enthusiastic and preceptive style of didactic guidance as well as his strong desire to seek further knowledge have contributed much energy to all aspects of this work. 1 Background and Introduction Reflection seismologists and stratigraphic geologists are generally interested in the problem of estimating subsurface acoustic impedance from normal incidence seismograms. The estimation is usually carried out in two steps. First the measured seismogram is processed to resemble the subsurface reflectivity function; and then this result is converted into an estimate of acoustic impedance. The second step normally requires the introduction of low frequency impedance information obtained either from nearby well logs (Galbraith and Millington, 1979), or from a particular velocity analysis carried out on a suite of common midpoint seismograms (Lavergne and Willm, 1977; Lindseth, 1979). This step has been considered necessary because of the inherent bandlimited nature of seismic recordings. Because of experimental limitations, seismic reflection information is typically recoverable only within the frequency band, range of 10-60 Hz. The low frequency information (0-lOHz) is vital to an accurate interpretation of a recovered impedance log. In this thesis, it is proposed that the estimation of acoustic impedance be carried out with a somewhat different philosophy; in a way which essentially follows the Backus-Gilbert framework of linear inverse theory (Oldenburg, 2 1981). First, the measured seismogram is modelled and then appraised to obtain unique information about the subsurface reflectivity function. This unique information is called 'reflectivity averages' and the method used here to obtain it is called 'appraisal deconvolution'. A broadband reflectivity function is then constructed using this unique information along with some physical constraints and a suitable algorithm. From this result, an interpretable broadband estimate of subsurface impedance may be obtained. The importance of this work lies primarily in the insights gained by approaching impedance recovery using a linear inverse formalism; and in the efficacy of the various procedures which have been refined to construct broadband reflectivity and impedance models. A brief review of some important aspects of the seismic reflection method and associated terminology is now presented. The basic seismic reflection survey entails recording subsurface boundary reflections arising from an artificial disturbance (shot or source) created near the ground surface. These reflections are recorded at a line of equally spaced receivers (geophones) which extend a short distance from the source point (see figure 1). A disturbance arriving at the geophone from a single reflection will be called the source wavelet (shown in the 3 seismic traces as a coherent wiggle - figure 1). According to the simple laws of reflection (for horizontal boundaries) the reflection point (RP) is located vertically beneath the source-receiver midpoint. Each disturbance created and recorded will have a single source geometry similar to that shown in figure 1. By overlapping the single source geometry at successive shots, it is possible to arrange such that different source-receiver offsets have a common midpoint (see figure 2). In the case of horizontal boundaries, these common midpoint source- receiver depth Boundary Source Wavelet .Surface Resulting Seismic Traces Fig. 1. Diagram showing a typical source-receiver geometry used in seismic reflection prospecting and the recorded seismic trace from each receiver. 4 pairs will be associated with a common depth point (CDP) as shown in figure 2. A collection of seismic traces for which the source-receiver offsets have a common midpoint is called a CDP gather (see figure 3). The number of traces in the gather defines the magnitude of multifold coverage obtained by the overlapping shots. A CDP gather with N traces provides N-fold CDP coverage and is said to cover the subsurface by NxlOO percent. The stacking procedure which will be described in the following paragraph attempts to make use of this redundant information to increase the Midpoint Resulting Traces Normal Incidence Ray Path COP Fig. 2. Common midpoint source-receiver pairs resulting from multifold coverage and the recorded trace from each receiver. 5 signal to noise ratio of the data. Observing the CDP gather in figure 3 it is evident that there is differential travel time (or normal moveout-NMO) to each reflective event across the gather for different source-receiver offsets. After correcting this differential travel time in each trace to match the travel time of the normal incidence ray, all traces may be summed together algebraically (i.e., stacked) enhancing the signal to noise ratio of the result (assuming random noise). The resulting stacked trace is then plotted at the source-receiver A**^"* Midpoint Stacked Result CDP Gather NMO Corrected Fig. 3- Processing sequence on a common depth point gather, resulting in a stacked normal incidence seismic trace. 6 midpoint, thereby producing a normal incidence seismogram in which the reflective event is located vertically beneath the source-receiver midpoint. This seismogram is in suitable form for treatment by the methods described in this thesis. Obtaining unique reflectivity information from- a normal incidence seismogram involves many assumptions and sequential processing steps. A proper treatment of the observed seismogram must include the effects of: (1) energy losses through geometrical spreading, anelastic absorption, and transmission losses across boundaries in a layered media (2) dispersion (3) source type and coupling (4) field geometries and array responses (5) multiple reflections and other coherent noise (6) the mathematical model imposed on. the seismogram (7) other The solutions of these complications will not be adressed in this thesis. It is assumed that the seismic data have been processed so that the seismogram, denoted by s(t), may be 7 considered to result from the convolution of the subsurface reflectivity function r(t) with the source wavelet w(t) plus random noise n(t). That is, the normal incidence seismogram is described by the following convolutional model 5CO-ra)* wet) +• nco where the symbol * denotes the convolution operation. Acoustic impedance estimation from a normal incidence seismogram will be developed from basic concepts. A complete treatment of the problem begins in chapter 1 with the developement of linear relationships between reflectivity and acoustic impedance according to normal 'incidence, plane wave theory. A discussion of the convolutional model is then provided in chapter 2 to establish simple grounds for understanding the complete experimental problem. In chapter 3, the traditional method of deconvolving the seismogram when properties of the source wavelet are known is viewed under the umbrella of inverse theory. Two distinct methods of reflectivity construction will be presented in chapter 4 and their efficacy shown using synthetic seimic data. The first construction method is a maximum entropy method which uses an autoregressive representation of a small portion of the reflectivity spectrum to predict spectral values outside that small 8 portion. The second and most versatile construction method is a linear programming approach which utilizes the unique spectral information from an appraisal deconvolution and provides a reflectivity function which has a minimum 1| norm, leading naturally to an impedance model with blocky character. The final chapter is devoted to examples of the above construction methods using real seismic data. Also, the Karhunen-Loeve transformation, a weighted averaging procedure which plays a strong role in stablizing the impedance reconstructions, is described and then utilized on real data. Errors in the data and the uniqueness of constructed reflectivity models will play important roles in estimating the acoustic impedance function and in assessing its uniqueness. These features will be developed in some detail throughout this thesis. The generally accepted idea that low frequency impedance information must be supplied from well logs or velocity analyses because of the bandlimited nature of the seismogram, is challenged. If accurate, bandlimited reflectivity information can be recovered from the seismic trace, then an interpretable, broadband impedance model may be constructed using the two algorithms presented in this thesis. 9 Reflectivity and Impedance A review of the basic concepts in normal incidence, plane wave theory will now be presented to establish the relationship between reflectivity and acoustic impedance to be used throughout this thesis. At a compositional boundary, incident and reflected wave amplitudes are related by a reflection coefficient r which depends on the density p and velocity v of the adjacent media where A; is the incident and Ar is the reflected wave amplitude (see figure 4). I.I Aj Ar Transmitted Energy Fig. 4. Normal incidence reflection from a compositional boundary. 10 The product of density and velocity is the acoustic impedance of the medium and will be denoted z. Extending this concept to a layered subsurface with many boundaries, an expression for the reflectivity series is obtained in which there is a reflection coefficient at each boundary described by (see figure 5) Rearranging this expression, a simple formula for the computation of acoustic impedance results where z, is the acoustic impedance in the initial layer. Expressions 1.2 and 1.3 represent the main objectives of the reflection seismologist. First the seismologist must obtain an accurate reflectivity measurement. This measurement is then converted into an estimate of subsurface acoustic impedance so that rock types and perhaps other information may be inferred. Defining as the logarithm of normalized acoustic impedance, where the normalization is with respect to the initial acoustic impedance Zj , expression 1.3 becomes 11 Surface Layer k»l zk»l Reflectivity Series iiiiiiiiiri Layer 1 r II iiii i /1 > Layer 2 z2 •1 r Layer 3 " '2 t-• • « • - 13 i i i i i t i /1 Layer k II / / ; i / i *k ' rK-l Acoustic Impedance Fig. 5. Subsurface plane layered model showing the reflectivity series and acoustic impedance funtion defined for this type of model. In order to simplify further discussion of it will be from here on referred to as impedance, distinguishing it from 2 which is the acoustic impedance. For |r| <1 , 00 FIN-1 12 This approximation is valid for |r| <.25 as shown in figure 6. Using this simplification, expression 1.4 reduces to 1= 1 which gives a linear relationship between impedance and reflectivity. This expression may also be written or In the time domain it can be shown that relation 1.6 has the continuous analogue (eg. Peterson et al, 1955) 4^ = 2 ret) ,7 or ^(O = 1 $trcu) du i.g 0 where the reflectivity function r(t) and the reflection coefficients rj are related by rv rc-O-Z r^a-r,) 1-1 13 where N is the number of boundaries and S(t- ) is the Dirac delta function located at the time occurence of each reflection coefficient. Equation 1.8 may be also written as (eg., Peacock, 1979) ^(O = 2 ret) * uit) 1.1 where u(t) is the unit step function. Convolving the reflectivity function with the unit step performs the role of integration. Fig. 6. The error between f( and the linear approximation fx for values of |r| less than .9. 14 In this thesis, the linear relationships between reflectivity and impedance (equations 1.5-1.9) will be used in place of the nonlinear relationship (equation 1.4), to simplify further mathematical transformations. It must be remembered however, that all linear relations have the restriction that reflection coefficients |r| must be less than .25. If a reflection coefficient exceeds .25 then the impedance should be calculated via equation 1.4. 15 2 The Convolutional Model The purpose of this chapter is to briefly explain the convolutional model in seismic reflection theory and to make clear the limitations of conventional deconvolution (inverse filtering / spectral division) techniques at recovering broadband reflectivity information from a seismic trace. A seismic trace is generally looked upon as a sequence of overlapping source wavelets that result from a variety of reflecting boundaries (Ricker, 1940 and 1953). Information about the subsurface geology is hidden in the traveltimes and amplitudes of these wavelets, that is, they conceal the subsurface reflectivity function. The ultimate goal of seismic deconvolution is to replace each overlapping wavelet with a corresponding reflection coefficient representative of a subsurface boundary. It will be shown that this goal can be only partially fulfilled using a conventional deconvolution technique. Several good examples of conventional deconvolution techniques can be found in Deconvolution, Geophysics Reprint Series 1, volumes I and II . The simple convolutional model without additive noise states that the measured seismogram s(t) (whether it be the result before or after the stacking procedure), is equal to the reflectivity function r(t) convolved with a source wavelet w(t), that is, 16 The source wavelet is defined to take on a variety of characteristics depending on the experimental procedures used and the .physical effects of wave propogation. The source wavelet is sometimes defined as the particle velocity of a single reflected disturbance as recorded in time by a geophone. However, it is necessary to include processing effects in the wavelet when considering a deconvolution of stacked seismic data. Probable effects contributing to the source wavelet are summarized below: (a) source type and coupling (b) transmission effects (c) field geometries and geophone arrays (d) ghosting effect (e) near surface weathering (f) instrument response (g) processing effects (eg., stacking) (h) etc. By following a convolutional model where all experimental effects listed above are attached to the source wavelet, a strong foundation is set for recovering an estimate of subsurface reflectivity from the measured seismogram. Reflection seismologists realize the importance of the source wavelet in seismic data processing and are actively engaged in devising improved methods which provide source 17 wavelet estimates. These methods may be classified into three groups: (1) direct knowledge of the source signature (eg., vibrating source) (2) direct measurement of the source signature (Mayne and Quay, 1971; Kramer et al, 1968) (3) estimation of a wavelet from the measured seismogram based on assumed statistical or physical properties of the reflectivity function and wavelet (Robinson, 1967; Ulrych, 1971; Otis and Smith, 1977; Lines and Ulrych, 1977, and Oldenburg, et al., 1981). For all discussions relating to the convolutional model, it will be assumed that the source wavelet is known and that determination of a broadband reflectivity function is desired. This section will conclude with a pictorial representation of the convolutional model and an example of spectral deconvolution which indicates what reflectivity information may be uniquely recovered from the seismogram. Figure 7a shows the result of a seismic experiment (where the data are measured in time) in which the source wavelet 18 is a delta function. The convolution of a delta function with subsurface reflection coefficients results in a perfect recovery of the true reflectivity function. Convolving the subsurface structure with a more realistic source wavelet leads to a poorly resolved, uninterpretable measurement of reflectivity as shown in figure 7b. This effect may be observed in the frequency domain by applying the convolution theorem to figure 7b. The convolution in the time domain becomes a multiplication of spectral components in the frequency domain as shown in figure 8a (where only the Subsurface Measured Reflectivity Source Reflectivity Function Wavelet Function Fig. 7. a) The convolution model in the time domain using an ideal wavelet and, b) using a more realistic wavelet. 19 amplitude spectra and positive frequencies are shown: phase effects have been left out for simplicity). Note that the bandlimited character of the measured reflectivity is due to the bandlimited nature of the source wavelet. In general, the reflectivity function is not bandlimited and thus it is seen that a complete recovery of the original reflectivity spectrum is inhibited by the spectral nulling effect of the wavelet. Unique reflectivity information may be obtained from the measured data by removing the spectral smoothing effects of the wavelet within its nonzero bandrange. In other words, a rough appraisal of uniqueness may be made by performing a bandlimited spectral division as shown in figure 8b. Since the wavelet contains insignificant energy at both low and high frequencies, no reflectivity information is recoverable at these frequencies. Only those frequencies which are contained within the bandwidth of the wavelet can be recovered (see figure 8b). This band of frequencies represents the only unique reflectivity information recoverable through spectral division or inverse filtering techniques. The nonuniqueness inherent in recovering broadband reflectivity is now apparent from the above analysis. There are infinitely many reflectivity functions which will satisfy the resultants of figures 8a and 8b. The only 20 information that all such functions must have in common is the band of frequencies shown on the right hand side of figure 8b. However, recovering broadband reflectivity information from seismic data can be a tractable problem if the missing frequencies are accurately replaced. Standard methods for replacing low frequencies involve well log or velocity analysis. Alternatively, a reflectivity construction technique, which utilizes the information provided by a single deconvolved trace, may be used to fill in the unknown frequencies. Two such reflectivity construction methods are presented in chapter 4. Reflectivity. Spectrum Wavelet Spectrum Measured Reflectivity Spectrum Measured Reflectivity Spectrum Wavelet Recovered Reflectivity Spectrum Fig. 8. a) The convolution model (fig. 7b) in the frequency domain and, b) deconvolution by spectral division. 21 3 Inverse Theory and Appra i sal Deconvolut ion The purpose of this section is to provide the reader with the basic concepts of inverse theory and its relation to deconvolution. This heuristic approach requires that the notation be somewhat simplified; a more rigourous treatment of this topic can be found in Oldenburg, (1981). This section is in most part a precis of that paper. The convolutional model of a seismogram neglecting the additive noise term-is written When the wavelet is known, equation 3.1 may be appraised in either the. time or frequency domain. Each of these approaches will be discussed seperately although the results under similar conditions are identical. A) Time domain appraisal deconvolution In the time domain, an inverse filter v(t) is desired that contracts the wavelet into a delta function; i.e., 3.1 3.2 But in the seismic case, w(t) is bandlimited, so there is no 22 filter which will contract the wavelet into a delta function. Finding an inverse filter which contracts the wavelet into a zero phase approximation to a delta function is the best alternative. By zero phase, it is meant that the approximated delta function has its energy peak located at the true delta function position and is symmetrical about its energy peak. This can be done by finding a filter which minimizes the integrated squared difference between the two sides of 3.2, that is, the filter minimizes where the optimum position t0 of the desired spike approximation may depend on the length or phase characteristics of the wavelet (Treitel and Robinson, 1966: Oldenburg, 1981). Obtaining the best v(t) by minimizing 0 , both sides of 3.1 are then convolved with the inverse filter to perform the appraisal deconvolution soo*vet) - rco*MO* VCO = rco* a CO. = <ra)> 3,4 23 where <r(t)> denotes reflectivity averages. Pictorially, the appraisal might look like: SCO VCO <rco> The convolution of the wavelet with the inverse filter produces an averaging function which is in some sense more localized than the original wavelet. For example: WCO VCO QCt) The averaging function plays a fundamental role in inverse theory because its character immediately indicates the resolving power of the wavelet. This type of averaging 24 function is often called a zero phase wavelet by seismic data processors. From 3.4, the convolution of the averaging function with the true reflectivity model will produce reflectivity averages. For example: Ail A rA I rco aco <r(t)> Knowledge about the true reflectivity, recoverable from the seismogram, is completely summarized by the average values <r(t)> and the averaging function a(t) unless extra information is provided about the form of r(t). Deconvolution utilizing the most resolved averaging function is seldom useful on noisy seismic data because the inverse filter will enhance the noise in the data to unacceptably large values. The noisy seismogram is written set) = ret) * wet) * rut) 25 Convolving this with an inverse filter v(t) produces SCO * VCO = rCO* wet) * vft) + nco * vet) = rot)* act) + net)* VCt) = <ra)7 -h £<rft)> 3,5" where £<r(t)> is an error term arising from the additive noise. In this case it will be desirable to find an inverse filter which produces the most resolved averaging function and also keeps the error term minimum. Data errors can be accomodated in the inverse filter solution by minimizing a statistical norm of £<r(t)>. The statistical form of the error term which is mathematically complimentary to the functional form of the resolution measure 3.3, is its variance which is denoted By simultaneously minimizing the error variance Y and the resolution measure (f) , an optimum inverse filter v(t) may be derived. Let where O is called the tradeoff parameter and varies between 0 and IT/2. A tradeoff between resolution and statistical 26 error produced by the resulting inverse filter may be obtained by minimizing 3.6 using different values of 0 . When © = 0, the v(t) which minimizes 3.6 yields the most resolved averaging function, but the reflectivity averages derived from s(t)*v(t) will have the greatest uncertainty. As © increases, resolution will be lost but greater statistical accuracy will be achieved. By using this philosophy, an inverse filter may be derived at any value of © and that one which provides the best compromise between resolution and accuracy must be chosen. Since (f) is not calculable, it is convenient to define a resolution measure 'L' as being the inverse height of the averaging function at its energy peak. Then, for example, as Q increases from 0 to we may plot a tradeoff curve to monitor resolution L and variance Y. The result would look something like: increasing * eso Tradeoff Curve variance \ * \ \e >o y xe=*/2 decreasing resolution 27 In practical appraisal deconvolution, the inverse filter is derived at a nonzero value of the tradeoff parameter so that better statistical accuracy will be attained. Very often a dramatic decrease in the variance can be achieved with only a small sacrifice in resolution as indicated by a steep initial decline of the tradeoff curve. B) Frequency domain appraisal deconvolution Deconvolution is carried out in the frequency domain by spectral division. If S(f), R(f), and W(f) are respectively the Fourier transforms of s(t), r(t), and w(t), where the Fourier transform pair is defined as, -oo then the convolution theorem applied to 3.1 yields A spectral division gives 28 where W*(f) is the complex conjugate of W(f). It is well known that this equation will produce erroneous results at frequencies where W(f) is small if the data are inaccurate, so 3.7 is necessarily altered by adding a stablizing term to the denominator. For inaccurate data, the stablizing term has been derived in the frequency domain by Oldenburg, (1981). In a manner quantitatively consistent with the time domain formulation, the frequency domain appraisal equation may be written where o( is a constant dependent on the magnitude of random noise in the data . and is the tradeoff parameter. The instabilities in the spectral division will now be overcome if © is sufficiently large. Equation 3.8 can be written where V(f) is defined as the quantity in curly brackets. Taking the inverse Fourier transform of 3.9 yields where v(t) is the inverse filter in the time domain. Thus 29 the results of appraisal in the frequency domain are completely analagous to those in the time domain. Each appraisal carried out at a tradeoff value Q yields unique reflectivity averages, a statistical variance of those averages, and an averaging function. All appraisal deconvolutions presented in this work were carried out in the frequency domain because of the computational and conceptual ease of working in that domain. C) Impedance computation and appraisal An appraisal shows that the reflectivity averages, their statistical error, ' and the associated averaging function completely summarize our knowledge of the true reflectivity function when the wavelet is known. However, in the event that the output from an appraisal is used to compute an impedance function, some care must be taken to differentiate between an interpretable, broadband reflectivity model, and averages of the reflectivity model. An interpretable, broadband reflectivity model is one which produces an adequate fit to the data and provides a meaningful impedance estimate. In general, reflectivity averages cannot be successfully used alone to provide a useful impedance function because they lack important low frequencies. For example, integrating the basic model equation 3.1 by convolving it with twice the unit step 30 function u(t), gives Denoting the left hand side s^t), as the integrated seismogram and making use of equation 1.9 gives Next convolving both sides by the best inverse filter of the wavelet produces $&) * vct) = 7a) * woo *va) - '^tt)*atO = <7<0? This result indicates that by appraising the integrated seismogram, impedance averages, <^(t)>, are obtained. The convolution of the integrated data with the inverse filter produces an output equal to the convolution of the true 31 impedance model with the same averaging function obtained from the reflectivity appraisal. For example: act) The resolving capability of the averaging function, in this case, is not as important as its low frequency content. Since the frequency content of the averaging function is dependent on that of the wavelet, it is missing both low and high frequencies; and it is the lack of low frequencies in the averaging function which leads to the uninterpretable impedance estimate <^(t)> shown above. Clearly there are troughs and peaks in the impedance averages which do not exist in the true impedance model and the major structural trends are absent. Reflection seismologists' have generally known that bandlimited impedance averages, when considered as a final model, will usually not produce correct structural interpretations. Thus, techniques of recovering the low 32 frequency impedance information from well log or velocity analyses have been developed (Galbraith and Millington, 1979; Lindseth,1979). In addition (or alternatively), the operating philosophy of the seismologist should include the idea of constructing broadband models from the appraised seismic data, to help clarify the interpretations. In the next chapter, reflectivity averages obtained from an appraisal are used in two algorithms which construct broadband reflectivity (and thus impedance) estimates. A comment about uniqueness may be made by reflecting on the basic model equation s(t)=r(t)*w(t). There are infinitely many functions r(t) which satisfy this equation (given the bandlimited nature of the wavelet) and this number is further enlarged if we accept that the data are inaccurate and,permit functions which generate statistically allowable misfits to those data. This group of permissible functions may well give rise to a wide range of computed impedance functions, unless some way of constraining the nature of the constructed reflectivity functions is.devised. D) An example of appraisal deconvolution derived from an interpretation of well log data and represents a reasonable picture of general subsurface impedance in parts of Alberta. This impedance model will be The impedance shown in figure 9a was 33 used as a starting point for synthetic examples of appraisal deconvolution. A reflectivity function (panel b) was computed from the impedance model via equation 1.7 and convolved with a wavelet (panel c), to produce a synthetic seismogram s(t). Random noise having a STD of 20% of max|s(t)| was added to this result, providing the noisy seismogram shown in panel d. The time sampling interval for this example is 4msec. The corresponding frequency domain representation of this procedure is shown in panels e-g; where only the amplitude spectrum of positive frequencies are displayed. The appraisals, carried out in the frequency domain at seven different values of tradeoff parameter Q, are shown in figure 10. The relative resolution and variance measures 1 (where LQ simply refers to the first value computed). The normalized averaging functions for each value of © are plotted in the first column of figure 10. The normalized reflectivity averages along with their amplitude spectra are plotted in column 2; and the normalized impedance averages are shown in column 3. Note the change in character of the averages and in their spectral content as the tradeoff parameter is increased. Stability is achieved around where frequencies outside the bandrange of the wavelet have been significantly reduced in importance. defined by L/L~ and are listed in table 3^ 1-1 ( a D TIME (sec) 2i56 1OO ret") b (D 0.0 O.B 1 1 1 1 i i i i i 1RCOI Y— wa) c o oo _ o o f iwcoi a CO 0.0 4.0 - 1 i i i 3 FREQUENCY 125Hz iscol Fig. 9. Synthetic models for an example of appraisal deconvolution. 35 1 2 3 k A . A A A— JV— ) FHEOUCMCT 125 Ha norm a(l-©) norm <ttt;e)> lRtf)l norm <^(t;e)> Fig. 10. Appraisal de convolutions shown using seven values of tradeoff parameter 36 Although the resulting averages are undesireable for impedance interpretation, they are unique and may be utilized in a construction procedure which leads to broadband forms of reflectivity and impedance. The topic of reflectivity construction will be covered next. e LA0 .001 1,0 1 .0 ,01 1 .ST .1 7.0 .014 .sr 2.5" .002) i.o 3.5- .QOOiX 1.4 5.IS .0002S 1 S2 10.0 .00006Z Table 1. Relative resolution and variance measures for the appraisals shown in fig. 10. 37 4 Reflectivity Construct ion Seismic reflection analysis, thus far, has led to the formation of unique information about the reflectivity function based upon the method of appraisal deconvolution. It has been shown that the only unique information obtainable through use of the convolutional model is in the form of reflectivity averages <r(t)> = r(t)*a (t) , where a(t) is an averaging function. It was also shown that the only unique impedance information recoverable by appraisal methods is in the form of impedance averages < r*Jit)>= v£(t)*a(t) , where a(t) is the same averaging function as above. It was noted that impedance averages are unacceptable as a final model of subsurface structure because of the bandlimited nature of a(t). Since appraisal methods (i.e., inverse filtering methods) must always fail to generate an acceptable impedance function if the source wavelet is missing low frequencies, it is worthwhile investigating whether a particular construction method might successfully fill in the missing frequency information about r(t). In fact, since the unique information about r(t) recoverable by appraisal methods represents only a bandlimited portion of its spectrum (eg., figures 8&10), the recovery of missing portions of R(f) will require some sort of construction technique. 38 The problem of broadband reflectivity recovery as described herein may be thought of as follows. First an appraisal of the seismogram is performed which replaces each overlapping wavelet present in the data with a zero phase averaging function. Then by use of a construction method, each zero phase averaging function is replaced by a corresponding reflection coefficient representative of a subsurface boundary. Two construction procedures will be described which have shown success in accomplishing this goal, but first an example is presented which emphasizes the nonuniqueness inherent in any reflectivity construction procedure. Consider the dual example shown in figure 11. In panel (a) is a blocky (or minimum structure) impedance function. A more realistic impedance function (a noisy version of the of the blocky impedance model) is shown in panel (b). Directly below each impedance model is the reflectivity function computed from each via equation 1.6. Adjacent to each reflectivity function is its associated amplitude spectrum. The bandlimited (10-50HZ) averages of rj (t) and r^(t) shown in panels (h) and (i), are very similar. Therefore, it is reasonable to expect that any construction algorithm carried out identically on both <rj(t)> and <r^(t)>, utilizing only the bandlimited information supplied by each, would lead to similar broadband reflectivity Fig. 11. The nonuniqueness inherent in reflectivity measurements. 40 models. Further, these constructed models may not resemble either (t) or ^(t) ^ no constrai-nt is placed on the form of the constructed output because there are infinitely many reflectivity models which have averages similar to those shown in panels (h) and (i). The size of model space may, however, be restricted by considering only a particular class of models. Both construction methods discussed in this thesis tend to favour, in theory, the construction of models similar to rj(t), that is, models with as few significant reflection coefficients as possible. The first of these methods, an autoregressive (AR) technique, makes the assumption that the reflectivity spectrum R(f) behaves in a predictable manner. It will be shown that this method can accurately reconstruct the missing low frequency portions of a reflectivity spectrum. A) Autoregressive extension of the reflectivity spectrum It is well known that the maximum entropy method (MEM) of power spectrum analysis, when applied to a short portion of a stationary data series, will in some cases, provide a more resolved spectral estimate than that from any conventional Fourier transform method. This is because conventional spectral estimators assume that the data are zero outside the known interval, whereas the MEM estimate assumes that the data are 'continued in a sensible way' 41 outside that known interval. By 'continued in a sensible way', it is meant that the process which produced the data series is modelled in a manner which allows prediction of unknown data. For example, only a fraction of one period of a sine wave need be known in order to infer its complete unending form. We must only specify that the process indicated by the known fraction is monochromatic. A general process which requires only a small portion of the data series to derive an adequate model for the whole is an autoregressive (AR) process. It can be shown that the MEM spectral estimator assumes an AR model for the stationary data series to be examined (eg., Ulrych and Bishop, 1975). The AR representation of a stationary data series x with zero mean, is P where e is a white noise series with zero mean and variance OQ . And the are prediction filter coefficients of order p. This expression indicates that a future value x^ can be predicted to within an error e^, from XK-1 ' xk~2 x k-p' th-at is, from 'p' past values of x. For this reason, equation 4.1 is called a forward prediction 42 representation. A reverse prediction expression may be written past value x^ can be predicted to within an error e£ , from xX+i 'x K*l' * * " 'x K*p' t^iat *s ' £rom 'P' future values of x. It follows, then, that for an AR process x, past and future values can be predicted from any set of N>p known values of x. If the prediction filter coefficients are determined by minimizing the variance of the error series, the reverse and forward solutions are identical. Thus, a single set of filter coefficients may be found by minimizing the sum of the forward and reverse error variances. This sum is estimated by the following forumula P H.2 where the are reverse prediction filter coefficients and er also has variance (J^. This representation shows that a N 43 Differentiating 4.3 with respect to the filter coefficients leads to the following matrix equation C^rq 7.7 where 4 - .-.p and r Co(Vjo(2 }. . . . o^p ) The o('s may be obtained by inverting equation 4.4. Alternatively, the prediction coefficients can be determined efficiently by minimizing 4.3 using the Burg-Levinson recursion technique described by Ulrych and Bishop (1975). Both solutions effect a least squares fit of an AR model of order 'p' to the data series x. The first method of constructing a broadband reflectivity function uses an AR model of the reflectivity spectrum to predict values of R(f) outside the frequency 44 band range of the source wavelet. Consider the basic model equation s(t)=r(t)*w(t) and its Fourier transform S(f)=R(f)W(f). If W(f) contains significant energy in the band of frequencies between fj and f2> then the reflectivity spectrum R(f)=S(f)/w(f) can be adequately recovered only in that frequency band. This band of recovered frequencies represents a small portion of the complex data series R(f). Thus, if R(f) can be modelled as an AR process of some order, then the recovered portion of R(f) may be used to estimate that AR model which will allow a sensible extension of R(f) outside frequencies fj and f^. Importantly, a proper choice of order is necessary for good prediction results. For the reflectivity construction problem, it will be shown later how the choice of p is related to the number of major subsurface boundaries or reflection coefficients. Attempting to show that the reflectivity spectrum R(f)' may be represented as an AR process, the reflectivity r(t) is first written in terms of a sampling function n-o where At is the sampling interval, N is the number of data samples, rn is the reflection coefficient at each sample interval, and ^(t-nat) is the Dirac sampling function. The 45 discrete Fourier transform of 4.5 is Rj=|,rnc^,r^;M.J....W-. 16 where the subscript j refers to the jth frequency fy=j/N4t and where values of j larger than N/2 refer to negative frequencies (N/2 being the Nyquist index). The inverse discrete Fourier transform is defined as Expression 4.6 may also be written KM ruo rV-l These relations will be referred to as reflectivity constraint equations. From equation 4.6, the reflectivity spectrum is a sum of complex sinusoids. Ulrych (1973), and Ulrych and Clayton (1976) have shown that harmonic functions may be adequately modelled as AR processes. The order of a single complex sinusoid is p=1, and in general, the order of a complex 46 harmonic process is equal to the number of complex sinusoids of different period included in the data. From 4.6, the number of complex sinusoids of different period (T^=N&t/n) included in the reflectivity spectrum is equal to the number of nonzero reflection coefficients r^. Thus, the order 'p' of the discrete reflectivity spectrum and the number of reflection coefficients are directly related. For computational considerations, the order of R must be less than the number of spectral values recovered between fj and f2« This puts a limit on the number of significant reflection coefficients which may exist within the data window being considered. If 20 complex frequency values of R are known, then there should be less than 20 significant reflection coefficients contributing energy to those known frequencies to insure a proper determination of the AR model. The AR representation of the reflectivity spectrum may be written P where all terms are complex. This expression is used to extend values of R to frequencies outside the known range f^ and f^. Part of the success of this AR extension method depends on accurately determining the order 'p' of the 47 reflectivity spectrum, i.e., determining the optimum number of significant reflection coefficients expected in the constructed result given M known complex frequency values. Optimum criteria for determining the order of a data series have been discussed by Ulrych and Bishop (1975), and Ulrych and Clayton (1976). However, choosing an order of 2/3 to 3/4 the number of recovered frequency samples provides good results if that number exceeds the number of significant reflection coefficients expected in the solution (an averaging procedure is described in the next chapter which stablizes the choice of AR order). Having chosen a good order, the prediction filter coefficients are determined and then equation 4.8 is implemented in a forward manner to extend the high frequency end of R , and in a reverse manner to fill in the missing low frequencies. The synthetic example shown in figure 12 demonstrates the ability of the AR technique to recover low frequency impedance structure. An ideal starting point for the AR method is an accurate portion of a reflectivity spectrum. The 20 estimates of R between 10 and 50 Hz shown in panel (a), represent an accurate portion of the true spectrum shown in panel (j). Panels (b) and (c) show the corresponding bandlimited reflectivity and impedance functions. The bandlimited and true impedance functions are compared in panel (c). 48 CD O ~ o o 10 1 Jo 1 1 1 . AA A (L In .. 1 \l \ \ vvu b) <ra)>. CD o ~ • • w — rl o d) 1 1 1 1 1 CD o — 1 CJ o I 1 1 1 1 1 1 i; CD O ~ CD „—.i,y -36 TIME IO r<° O i i i i i i fREQUENCY 125 |RCe)| 0 I) —•< / .36 TIME Fig. 12. Autoregressive reconstruction of a reflectivity function containing a few reflectors. 49 Panel (d) shows the result of extending the spectrum into the low frequency range only, using AR prediction and an order of 15 (which is greater than the number of spikes in the true reflectivity shown in panel (k)). Panel (e) illustrates that the resulting constructed reflectivity function has not improved in appearance, however, the resulting impedance function (panel (f)) contains the proper low frequency information needed to give it geologic meaning. Panels (g)-(i) show the result of predicting both the low and high missing frequencies. The reflectivity function (panel (h)) is plotted as a spike series simply because it has a complete spectrum within the Nyquist frequency. The resulting impedance structure (panel (i)) indicates that the predicted high frequency values are somewhat erroneous. This is a common occurrence among the various reflectivity models that we have tested using the AR method. It appears less profitable to predict the very high frequencies because of small uncertainties in 'determining the proper prediction filter coefficients. These uncertainties may propagate into large errors in the predicted spectral value as the prediction filter slides off the edge of the known data and begins to predict spectral values exclusively from previously predicted values. Fortunately, the small band of missing low frequencies can, in many cases, be accurately 50 recovered. Figure 13 shows the results of applying the AR method to a more complicated example. Panels (j)-(l) show the results of adding 10% white noise to the simple synthetic impedance model presented in the previous example. The 18 spectral estimates of R (shown in panel (a)), between 14 and 50Hz, are accurate. The results of predicting only the missing low frequencies are very encouraging as has been shown in panels (d)-(f). The results of also predicting the missing high frequencies (panels (g)-(i)), show a simpler and perhaps more interpretable reflectivity function but the resulting impedance estimate has not been improved. An order of 12 was chosen for this example. More examples of this method are presented in section C) of this chapter and in the final chapter where it is applied to real data. B) Broadband reflectivity construction using the ljnorm and a linear programming algorithm Broadband reflectivity recovery in the time domain by minimizing the lj norm leads naturally to a linear programming (LP) formulation (Claerbout and Muir, 1 973; Taylor et al, 1979), but lj norm minimization has also been used for reflectivity construction in the frequency domain. Levy and Fullagar (1981) have formulated a procedure to 51 o 114 I J 1 1 1 b) <rC0> o~l I I I I I rcco 9) i 1 1 1 1 rcCt) 0 ; ?c(0 36 .36 i 1 r—i 6 f REQUENCt 125 j) IRCO| ret) 0 Fig. 13. Autoregressive reconstruction of a reflectivity function containing many reflectors. 52 construct a reflectivity function containing a minimum number of reflection coefficients. Their method, carried out in the frequency domain, assumes only that processing of the seismogram has resulted in unique reflectivity information (reflectivity averages), i.e., a bandlimited portion of the true reflectivity spectrum has been recovered. This method has an advantage in that the solution is constructed from only the most reliably known spectral values of R(f). However, as formulated, it does not take full advantage of the unique reflectivity averages. An outline of Levy and Fullagar's general LP approach will be followed by an explanation of how their approach can be made more efficient utilizing information from an appraisal of the seismogram and information from velocity or well log analyses. (a) General LP formulation An outline of this method begins with the sampled reflectivity function given previously as The objective is then to find that set of reflection coefficients rn which satisfy the constraint equations 4.7, N-l 53 and which have a minimum lf norm, that is, llrcOU, -2 lrn| is minimum Minimization of the 1| norm will produce a set of reflection coefficients with a minimum number of nonzero values. As an aside, utilizing expression 1.7, equation 4.9 may be rewritten Minimizing the 1( norm of the reflectivity function is thus equivalent to minimizing the 1| norm of impedance gradients. This minimization will produce an impedance model with a minimum of structural variations. One way to formulate the LP solution is by first obtaining estimates of the reflectivity spectrum R by spectral division of the seismogram where j refers to the jth frequency as before. But only those spectral components, where W is sufficiently large, are used. Estimates of R may also be obtained through many other methods of zero phase deconvolution. For example, they may be obtained from an appraisal deconvolution. vo Rj=Sj/Wj III 54 Reliably determined spectral values are then substituted into the constraint equations 4.7 and an LP algorithm is selected to find those reflection coefficients r^ which satisfy equations 4.7 and which have a minimum lj norm. Note that if both 4.11 and 4.7 are satisfied in the above procedure, the constructed reflectivity function, denoted by rc(t), should adequately reproduce the observed seismogram i.e., s(t)2Vc(t)*w(t). Most LP algorithms are designed to locate only positive unknowns, whereas both positive and negative reflection coefficients are expected. The solution to this problem is to express each r^ as the difference between two positive quant i ties, i.e., This relation is then substituted into equations 4.9 and 4.7, and an LP algorithm will search for the positive unknowns a^ and bn. The number of unknowns is doubled by this procedure, thereby increasing the computing costs. However, if no extra information can be incorporated into the LP algorithm, the formulation outlined above must be followed. 55 No attempt is made in the above formulation to utilize that unique information which could be supplied by an appraisal of the seismogram and by well log or velocity analyses. This information, when available, can make a valuable contribution to the LP solution. Assimilating reflectivity averages and impedance estimates into the algorithm leads to a more efficient and reliable solution of the LP reflectivity construction problem. (b) The constrained LP solution In the constrained LP approach, the objective is to find that set of reflection coefficients r^ which satisfy the constraint equations, and which have a minimum weighted 1^ norm, that is, N-l where the are weights meant to emphasize those . rn which are prominent in the average values <r(t)>. A natural weighting function for this approach is the arithmetic inverse of the reflectivity averages, <r(t)>~'. This is because the LP algorithm using lj norm minimization, will tend to reduce an r^ if it has a large relative weight. Thus, if the average values <r(t)> indicate that a 56 particular reflection coefficient should be near zero, that interpretation is assisted in the LP algorithm by giving it where <r^> refers to the digitized averages, and q is a weighting exponent ranging from 0 (no weighting) to 2 or more (0=1 or 2 is usually sufficient to provide good results). The reflectivity averages may also be used to provide polarity information to the LP algorithm. Defining a polarity function as a large weight relative to other possible rn's. A natural form of 4.13 is then -l i( I if *rn> 2 O 57 the constraint equations may be rewritten Real [Rj J -1 r1 coiC2TTi n/w) r3nC<v) OM r\=0 where r^ =sgn (<rn> ) r^ and j ranges over only the most reliable frequencies as before. Then the LP algorithm may search for the strictly positive unknowns r^ , eliminating need for the difference equation 4.12. The true polarity and magnitude of the constructed reflection coefficients may be recovered by applying rA = sgn (<rA> ) r^ . Additional constraint equations may be added to those in 4.16 if impedance estimates are available from velocity or well log analyses. From equation 1.5, impedance constraint equations may be written -^1 S9h«rn>)rn' 4,17 n-o where k+1 is the index of the impedance that is known and equation 4.17 is in the proper form for input to the constrained LP algorithm. If only the basement impedance is known (i.e., ), this relation becomes essentially a 58 constraint on the dc frequency R^. For example, the following two equations provide identical constraints: rV-t i\=o The constrained LP algorithm is summarized in figure 14. REFLECTIVITY Weighting Information Polarity Information AVERAGES Reliable Spectral Estimates R Impedance Estimates Objective Function Constraint Equations LP Algorithm Constructed Reflectivity Fig. 1^-. Flow diagram for the constrained LP algorithm. 59 Before presenting examples of the LP method, it is noted that by using the reflectivity averages to constrain the construction problem, certain restrictions are automatically placed on the resolution of dipole reflection coefficients. For example, a thin shale wedged in a porous sand formation might result in the following reflectivity funct ion Experimentally, the following reflectivity averages would be recovered (from equation 3.4) rCti <fCO> Therefore, this dipole can not be resolved by either appraisal or construction methods. However, if such dipoles are sufficiently seperated to be resolved by the appraisal, 60 then these features, complete with the low frequency information, may be recovered through construction procedures. (c) Data errors and the LP solution In practical cases, the observed seismogram contains a certain amount of noise which will affect subsequent processing. This data noise translates into an error for each recovered spectral value R^' of the reflectivity spectrum. In the presence of noise, it is undesireable to solve the constraint equations exactly. Two ways of handling errors in estimates of R have been discussed for this LP approach by Levy and Fullagar (1981) and will be briefly presented here. 1) Inequality constraints Most LP algorithms allow inequality constraints, where each constraint equation is required to be satisfied within a specified error bound. The reflectivity constraint equations may be easily converted into inequality constraints as follows ±Real[R.i] + Ej > fL^COLbxtin/H) 61 where E.j and *C>j are specified uncertainties at the jth frequency. In this inequality formulation, the LP algorithm is required to minimize the 1| norm of the r^, while satisfying the inequality relations 4.20. Error terms E^ and c-j can be most simply chosen empirically by assuming a data noise level and then requiring the uncertainties to be a fraction of this. Note that the number of constraint equations is doubled in this procedure. 2) Equality constraints In the process of determining those reflection coefficients which have a minimum ljnorm, the LP algorithm may also directly determine an estimate of the errors E-j and £^ at each frequency. This may be done by retaining the equality constraint relations, but introducing the uncertainties as additional unknowns J J r\=o ImqcjLR^ - £j -I r^sinUTfjn/AO 4,21 Also, an additional equation which constrains a statistical r 62 norm of the uncertainties is necessary to control the levels of noise recovered in the LP solution. The details of a constraint on the ljnorm of the error terms is covered by Levy and Fullagar (1981). Then the LP algorithm will search for those r^ which have a minimum lj norm and in so doing also determine error terms at each frequency which favour the minimization. Of course, the errors and can be positive or negative, so they must be cast in the form of 4.12 for input to the LP algorithm. Note that the number of unknowns is increased by 4 at each frequency, but the number of constraint equations needed in this approach is nearly half the number used in the inequality formulation. For the following examples of LP construction, an inequality approach was adopted. Synthetic data was generated using a zero phase wavelet (an averaging function) estimated from the correlated and stacked seimic data shown in figure 15 (both the data and wavelet estimate were provided by Imperial Oil of Canada). The wavelet and its amplitude spectrum are shown in figure 16. The first example of constrained LP construction is shown in figure 17. The reflectivity function shown in panel (b) was computed from the impedance model in panel (a). The reflectivity was then convolved with the zero phase wavelet to produce the synthetic seismic trace in panel (c). The appraisal deconvolution shown in panel (d) W 79E39 62921 63 Fig. 15. Deconvolved seismic section; represents subsurface structure in parts of western Alberta. \f 256sec w(t) IWCOI I25h; Fig. 16. Zero phase wavelet (averaging function) estimated from the data inside the box in fig. 15, and its amplitude spectrum. 64 <0 b) <r(t;©=0> 0 T .72 72 Fig. 1?. Linear programming reconstruction of a reflectivitv function containing a few reflectors. 65 indicates little improvement over the original data because the zero phase wavelet is already an ideal averaging function. The constructed reflectivity function in panel (e) was recovered by supplying the LP algorithm with accurate spectral values between 12 and 35 Hz. The averages <r(t)>'were used to weight the objective function and polarity constraints were also provided. The resulting impedance function, shown in panel (f), compares very well with the true impedance. No impedance constraints were included in this example. Figure 18 shows the results of applying the LP method using a noisy version of the impedance log from the previous example. The noisy log is shown in panel (a). The complexity of the resulting reflectivity and deconvolution (panels (b) and (d)) present a difficult problem for the LP algorithm since the method is best suited to find simple reflectivity functions. Accurate values of R between 10 and 35 Hz, supplied to the LP algorithm along with polarity constraints and weighting (from <r(t)>""') lead to the recovered reflectivity function in panel (e). The resulting impedance function (panel (f)) shows that the construction was reasonably successful at recovering the major trends of the original impedance. By including two impedance constraints, one midway and one on the basement impedance, a good recovery of 66 Fig. 18. Linear programming reconstruction of a reflectivity function containing many reflectors. 67 the original impedance structure was obtained as shown in panel (h). Further examples of this method are presented next. C) Stability and efficiency considerations Before applying the AR and LP construction methods to real seismic data, it is worthwhile investigating their performance on synthetic examples under the following complications. 1) the presence of additive noise in the data (both methods) 2) inaccurate knowledge of the source wavelet (both methods) 3) variance of the known frequency range fj to f*2 (both methods) 4) sensitivity of AR order to the number of reflection coefficients (AR method only) 5) effects of different weighting functions (LP method only) Figure 19 shows the synthetic data which will be the starting point for examples of the above considerations. The sampling interval is 4msec. For the first and second complications above., constructions will be performed using two slightly different spectral estimates. Since the wavelet is known, estimates 68 Fig. 19. Synthetic test data. Wavelet has a band range of 10-50Hz. of R between frequencies fj and fj may be obtained through a bandlimited spectral division; S(f)/W(f). Or estimates of R may be taken from an appraisal of the seismogram; S(f)V(f). More reliable estimates of R may be obtained from the appraisal because insignificant or noise contaminated frequencies have been winnowed from these results. Recalling equation 3.4, the appraised seismic trace may be written In the frequency domain this becomes DC0=RC0A(0 = 5(0V«) If A(f) is zero phase and approximately unity within the 69 bandrange to then the discrete spectral estimates D^', may be equated to R.J between fj and for use in the construction algorithms. For the first example, various amounts of random noise were added to the original synthetic trace s(t) in figure 19. The noisy seismograms are shown in panel (a) of figure 20 along with appraisal deconvolutions carried out on each trace at © =1. The amount of noise added to the data is indicated on the center line. The AR constructions are shown in panel (b). Fourteen recovered spectral values between 14 and 40Hz were used to predict the missing low frequencies and the high frequencies to 125Hz using an order of 11. The results on the left were obtained using spectral estimates obtained from .S(f)/W(f) and the results on the right were obtained using appraisal estimates S(f)V(f). Straight spectral division by the known wavelet has provided sharper reflectivity constructions, however, the resulting impedance models are not better than those obtained by using appraisal estimates. The LP constructions are shown in panel (c). Spectral values between 14 and 40Hz were derived from S(f)/W(f) to obtain the results on the left; and spectral values were derived from S(f)V(f) to obtain the constructions on the right. Clearly, using the appraisal estimates provides better LP constructions. 70 S -£o * $ < >5 •p -p o co rH CH (1) J-i CH o • CO CD "O CO o -p CD E o C CD > H! -P • H cd cc CC CH < O OJ 0) -p «H O o c CD to CD U -Ol CD O C CD rt .C e -p j-i o c CH -H CD C PH O • H -P • O O 3 CM t-i -P • CO W C •H O O 71 The second example, shown in figure 21, demonstrates the effect of inaccurate knowledge of the wavelet. Trace 1 in panel (a) shows the true mixed phase wavelet used to generate the original data and trace 5 shows a zero phase version of this wavelet which is used as an inaccurate wavelet. An appraisal of the original data using the inaccurate wavelet is shown in trace 6 (compare this with trace 3; an appraisal using the true wavelet). The AR constructions are shown panel (b). Again, 14 recovered spectral values between 14 and 40Hz were used to predict the missing low frequencies and the high frequencies to 125Hz using 5 different orders as indicated. The results using spectral division estimates S(f)/W(f) shown on the left are somewhat better than the results using appraisal estimates S(f)V(f) shown on the right. The LP constructions are shown in panel (c). The same known frequency range was used and different weighting exponents 'q' were tested. Also for the bottom 2 reflectivity constructions, a constraint was placed on the basement impedance. Again, using the appraisal output S(f)V(f), provides much better LP constructions. The third example tests the effects of varying the known frequency range f| to f2« The AR and LP results shown in figure 22 were obtained using accurate spectral estimates of R(f) in the ranges indicated. The LP constructions shown 1 V wet) 2 ——-v— ret) 3 ^^^^ sCt) u -JWW^ <rct;d2i)> ^ Y wet) 6 --A^vA^vV1- <rCt;© = 0? Q) AR —*—13 —»Y* A~vy~ orcier 0 .4 LP 1 c) Fig. 21. Performance of the construction methods using an innaccurate wavelet. 73 on the right remain stable, but the AR results shown on the left fail when the number of known frequency values becomes close to, and hence, the order becomes less than the number of spikes (9) in the original reflectivity r(t). The fourth example shown in figure 23 demonstrates the effect of varying the AR order using accurate and inaccurate data. Twenty spectral values between 10 and 50Hz were obtained from a spectral division of: 1) accurate data and 2) data with 10% added random noise. The AR constructions obtained using accurate spectral values are shown on the left and the constructions using inaccurate spectral values are shown on the right. The normalized squared error (NSE), plotted vs order, between the true and constructed models and defined by are shown at the bottom of figure 23. The NSE indicate that a long prediction filter (i.e., a high order) leads to construction may be obtained if the order is greater than the number of significant spikes expected in the solution. Because a greater range of reliable frequency values are spanned by a longer filter, a more reliable prediction of stable reconstructions. In other words, a good AR 74 —>^X-^j-/- roo -15 —K^Ji^A^^ \0-SOH -12—ibv^vX^^/u |0-«fO -7—V^-^N^^^-VSA- 10-30 -15—AAVl—A^yu^ IS-^O -10—Aj\rwwvw-~p^ /r-HcO -5—-A^VwV^ly-N^ |S"-30 - 1 LP I T T J-»-1v-WW o Fig. 22. Stability of reconstructions with respect to the frequency band used. 75 —h/t^A—A— —U^A-^Y^ ./L-J^-JU ret) H 7 1 10 I I M 11 IM I S" 17 IS order fNyVwA r^^yAvA~ Ax f/OlSE" FREE" \0% NOISE order Fig. 23. Stability of the AR method with respect to order or the stability of order with respect to the number of reflection coefficients; using accurate and innaccurate data. 76 unknown values may result. However, if the filter is too long (close to the length of the known data), noise on the data will also be predicted, perhaps giving udesireable results. This example suggests choosing the largest possible order consistent with assumed noise levels in the spectral data. If the signal to noise ratio of the data is very high, choosing a large order should lead to good results. Also, the NSE for the impedance constructions indicate that, out of 9 boundaries in the true ^(t), only 6 boundaries are significant with respect to the order. The last example in this section compares the use of different weighting exponents on the objective function in the LP algorithm. Figure 24 shows LP constructions using 14 spectral values between 14 and 40Hz and 3 different weighting exponents. The constructions are all very good but the ones' using objective function weighting (q=l,2) provide a good solution in 1/3 of the computing time required with no weighting (q=0). 7? JU- /I CPU TIME r<CO o I.Z Sec WO 1 .H 2 .H i 1 0 .4 Fig. 24. The accuracy and efficiency of the LP solution with aspect to the weighting exponent q, using accural 78 5 Multitrace Reflectivity Construction using Real Seismic Data A) The Real Data Data set A, shown in figure 25 along with a geological interpretation, consists of 24 adjacent traces taken from the deconvolved seismic section in figure 15 (vibrator data), where the data location is indicated by the box. Sonic logs are available near the seismic line, but they are at least 1 mile removed. Data set B, shown in figure 26, consists of a .5sec window of deconvolved explosion data (38 traces in all). ^ A close by velocity log is available for part of that time window near a trace location. The raw trace amplitudes for each data set represent measured voltages. Thus, if absolute reflection and impedance magnitudes are desired, then it is necessary to match the energy in each deconvolved trace with energy consistent with subsurface reflectivity. The standard method of energy scaling utilizes available well logs. A synthetic seismogram is constructed from a nearby velocity log (and a density log if available) and then the total energy in the deconvolved trace is matched to the energy in the synthetic. Energy scaling is necessary in the constrained LP approach if impedance or spectral constraints, obtained from well log or velocity analyses, Fig. 25. Data section A; deconvolved vibrator data along with a geological interpretation. Fig. 26. Data section B; deconvolved explosion data. 80 are included in the algorithm. Since each seismic trace in sections A and B has undergone an appraisal deconvolution, it represents reflectivity averages, <r(t)>=r(t)*a(t) , plus noise error £ <r(t)>=n(t)*v(t). The specific appraisal techniques used on each data set are unknown, but it is assumed that the original source wavelet has been localized into a zero phase averaging function, a(t)=w(t)*v(t) . Each deconvolved trace may be written Normally, A(f) is a smooth function and generally deviates from its ideal value of unity between f} and f^ (eg., fig. 16). Therefore, it may be necessary to first remove any spectral smoothing of R(f) produced by A(f) before supplying spectral estimates to a construction algorithm. Since an average A(f) is available for data set A (fig. 16), the sensitivity of the construction methods with respect to the smoothing effect of the averaging function will be tested using the following two slightly different estimates of cl(OrrCO*ac.o + ^rc^> In the frequency domain this is 5". I 81 R(f ) : 1) R(f) ^D(f)/A(f) (smoothing removed) 2) R(f) ~ D(f) for fj <. f >. Figure 27 shows construction results for data set A using the AR algorithm. The original data and the impedance averages obtained from them are shown in panel (a). Panel (b) shows the reflectivity and impedance constructions using spectral estimates D(f)/A(f); and panel (c) shows the results using estimates D(f). In both cases, 33 recovered frequencies between 13 and 46 Hz were used to predict the missing low frequencies and high frequencies to 62 Hz using an order of 24. No major differences exist in the two constructed sections, however, the results in panel (b), obtained by using estimates D(f)/A(f) are somewhat noisier. The LP constructions for data set A are shown in figure 28. The results using spectral estimates D(f)/A(f) are shown in panel (b); and the results using estimates D(f) are shown in panel (c). In both cases, reliable frequencies from 12 to 35 Hz were used and two impedance constraints were included; one midway and one on the basement impedance. Clearly, the constructed section in panel (c), obtained by using spectral estimates, D(f), provides the most consistent results. The above AR and LP construction results are generally Pig. 27. AR reflectivity and impedance reconstructions from section A; b) using unsmoothed, and c) using smoothed spectral estimates. Fig. 28. LP reflectivity and impedance reconstructions from section A; b) using unsmoothed, and c) using smoothed spectral estimates. 84 trace by trace comparable, establishing confidence in the reliability of both methods. Using smoothed spectral estimates, D(f), provided more consistent LP constructions and less noisy AR constructions. Thus, further constructions will be carried out using smoothed spectral estimates, D(f), obtained from the appraised data. Autoregressive constructions for data set B are shown in figure 29. The original data and the impedance averages obtained from them are shown in panel (a). Frequencies between 12 and 40 Hz were deemed reliable-by examining the amplitude spectra of several deconvolved traces. An order of 13 was used for the AR construction shown in panel (b), where the missing low frequencies and high frequencies to 65 Hz were predicted for each trace. The LP constructions for data set B are shown in figure 30. Frequencies between 12 and 40 Hz were used along with one impedance constraint (midway) to obtain the results shown in panel (b). The AR and LP constructed sections obtained from data set B provide similar and more localized reflectivity information. However, both constructed sections retain the noisy character of the original data. A somewhat noisy character is also observed in the constructed sections obtained from data set A. Although good interpretations could be made from these sections, the consistency of the Fig. 29. AR reflectivity and impedance reconstructions from section B using smoothed spectral estimates. CO Fig. 30. LP reflectivity and impedance reconstructions from section B using smoothed spectral estimates. CO ON 87 results must be improved before the AR and LP algorithms can be routinely applied to seismic data. The recovered broadband sections above were laterally inconsistent in some places, not because of the presence of structural changes but more likely due to the presence of noise in the data. A linear transformation known as the Karhunen-Loeve (K-L) transformation (eg. Ready and Wintz, 1973; Kramer and Mathews, 1956) can be used to extract the major similarities in a set of adjacent seismic traces for the purpose of imroving lateral consistency. Consider N traces of seismic data, where each trace SJ (t) is composed of a common signal x(t) plus a difference signal n;(t); i.e., If the nj(t) represent noise signals, then it is useful to find an x(t) which summarizes the most correlated components in the N traces. In general, the common signal may be written as a linear combination of the s,'(t) as follows B) The Karhunen-Loeve Transformation 5,2 88 The sj(t), i=1,...N, span at most N dimensions in Hilbert space and if the traces are highly correlated, they will have a preferred orientation in that space. This suggests rewriting the common trace as a linear combination of M orthonormal basis functions Y,-(t) where the directions of these basis functions are defined by the eigenvectors of the inner product or covariance matrix M 0-i The %*(t) may be defined as follows (see appendix A) r.6 where s ^ c^co, stcOj ....sNcty) \\is the jth eigenvalue arranged in descending order i.e., 89 X\ > >^-X|s/ an<3 b-j is the jth eigenvector of P . The in 5.4 may be found by minimizing the sum of integrated squared differences between x(t) and the sj(t), that is by minimizing i-.i 4 r-1 The result derived in appendix A is I—- N N ;fr 1 where Jjfj is the eigenvector matrix defined in appendix A. Substituting 5.6 and 5.8 into 5.4 gives M 1 where I-1 Each term in expansion 5.9 is called a principal component of the common trace. For example, 01 b, * £ is the first principal component and so on. If the N traces are highly correlated, then \\ and Yv will be relatively large and most of the information about • the common trace will be 90 contained in the first principal component. In fact if all N traces are identical, then the first principal component contains all of the information about the common trace and the other components are trivial. At best, the components of the first principal eigenvector provide optimum trace weights to construct the desired common trace. In any case, if M=N, i.e., if all principal components are summed, then 5.9 gives the mean trace (as shown in appendix A) A criterion could be established tp help determine how many principal components should be included in the common trace. However, if only the major similarities between 3 or 5 traces are desired, then only one or two principal components need be used to provide a common trace In this thesis, expression 5.9 is referred to as the K-L transformation and the common trace formed by this transformation is called a principal component or K-L trace. The K-L transformation can be applied in an overlapping manner to deconvolved seismic data, to reflectivity and impedance constructions, and to the choice of AR order, to help improve the reliability and interpretability of constructed reflectivity and impedance sections. For example, the left hand side of figure 31 shows the 91 results of applying AR prediction to the bandlimited reflectivity function, <r(t)> (12-46HZ), using several orders. The seventeen known frequencies between 12 and 46Hz were used to predict the missing low frequencies and the high frequencies to 125Hz. The results are slightly different for each order. Considering N=5 traces and using only the first principal component gives the K-L traces shown on the right hand side. These traces are very consistent over a wide range of orders indicating that the choice of AR order is less critical using the K-L approach. Figure 31 also shows the principle of the mixing algorithm. Traces 1-5 are used to obtain the top K-L trace, traces 2-6 are used to obtain the next K-L trace and so on. Figure 32 shows the results of applying the algorithm to the deconvolved data set B using N=3 and only the first principal component. The resulting principal component section shown below, indicates that significant information in the original data has been preserved and the noisy character is absent. Application of the K-L algorithm to an LP reflectivity construction (from figure 28) is shown in figure 33. The original constructed reflectivity and impedance sections are shown on the left and the principal component sections, obtained using N=3 traces and only the first principal component, are shown on the right. The interpretability of 92 11 10 8 ^A^~Ar^-^Y^^ 9 9 —Kj^j^f\—y^J^ 10 trace 8 7 order [f0; ^w] —IULJLJL^^L r(t) —^J^-y-^-1 A^vJ—A^^^Lv 16—j > L-K^A—^—Y^7^ 2 — ^-A-vy^Avv*rA^--^V^' 15 K-L 0 .4 Fig. 31* Application of the K-L transformation to stabilize the choice of AR order. Also showing the principle of the mixing algorithm. 93 K-L Fig. 32. Application of the K-L transformation and mixing algorithm to the noisy section B. 94 both sections has been markedly improved. Using this K-L algorithm, noisy traces can be discriminated against while preserving the most correlated information recovered by the construct ions. The final constructed reflectivity and impedance sections for data set A are shown in figure 34 along with the original geologic interpretation. Both the AR results, shown on the left, and the LP results, shown on the right, were obtained by applying the K-L transformation to the raw constructions shown previously using N=3 traces and the first principal component in the mixing algorithm. The AR and LP sections provide nearly identical interpretations. The positive impedance contrast at the PreCretaceous unconformity (which is a dominant feature in nearby well logs) clearly correlates in all traces and the oil bearing Redwater Reef shows up as an impedance low. Overall, it is easier to make structural interpretations from the recovered impedance sections. Comparisons between the original data and the reconstructions are shown in figures 35 and 36. The final constructed reflectivity and impedance sections for data set B are shown in figure 37. Again, these results were obtained by applying the K-L mixing algorithm to AR and LP constructions using N=3 traces and the first principal component. The available velocity log has been inserted into the impedance reconstructions for 95 Fig. 33- Application of the K-L algorithm to LP reflectivity and impedance reconstructions. AR LP If A PARK CoLOKArDO • WA&NM SUICROP BrAveR HILL LAKE" ELK POINT LEA PARK COLORADO WABAMU* SUBCROP ZZZZZZ RfD^ATR REEF Bwvep. HUL LAKE" ELK fO'NT Fig. 3>4-. The final reflectivity and impedance reconstructions from section A along with the original geologic interpretation. ON Fig' 35' Comparison of final reflectivity-reconstructions with the original deconvolved section (A). 98 Fig. 36. Comparison of final impedance reconstructions with impedance averages obtained from the deconvolved data. 99 comparison. The well log confirms the first .3sec of adjacent reconstructions and consistency between the AR and LP results provides additional confidence in the final models. Again, it is evident that subsurface structural interpretations would be aided by using the impedance section. Comparisons between the original data and the reconstructions are shown in figures 38 and 39. Conclusion The problem of estimating broadband acoustic impedance from reflection seismograms has been investigated using a linear inverse formalism. It was shown that a conventional deconvolution / inverse filtering technique (which recovers unique bandlimited reflectivity information from the seismogram) should precede the use of a particular construction algorithm (which attempts to recover the unknown frequency bands). Two construction methods presented in this thesis, the AR and the LP methods, have shown success at recovering missing spectral information. The K-L transformation was applied in real data examples to stabilize both methods in the presence of noise and impedance constraints provided additional stability to the LP solution. The AR method represents an inexpensive (computationally) way to obtain reconnaissance broadband AR LP Fig. 37. The final reflectivity and impedance reconstructions from section B, The given velocity log has been inserted for comparison. o o 101 Pig. 38. Comparison of final reflectivity reconstructions with the original deconvolved section (B). 102 Fig. 39. Comparison of final impedance reconstructions with impedance averages•obtained from the deconvolved data. 103 impedance information. Although the method works very well on synthetic and real data examples, at this time there is no way to constrain the output using information available from other analyses. Absolute impedance information, from available well logs should be utilized to improve the stability of AR constructions in the presence of noise. Other improvements to this method might include ways of obtaining the prediction filter coefficients by constraining the predicted dc frequency R0 or solving for the missing frequencies directly using maximum entropy considerations. The LP algorithm provides the most versatile approach to recovering broadband impedance. Robust in the presence of noise, this method allows much room for improvement; restricted only by the creativity of the designer and the efficiency of the linear programming algorithm. It is concluded that impedance sections obtained from both construction methods could greatly aid the interpretation of subsurface structure. The problem of reflection polarity may, in many cases, be solved by performing an impedance construction. For example, until the constructions for data set A were performed and the positive impedance contrast was identified, the polarity of the given data was not known. In fact, polarity was unknown for both data sets until preliminary constructions were performed. 104 REFERENCES Claerbout, J. F., and Muir, F., 1973, Robust Modelling with Erratic Data: Geophysics, V.38, P. 826-844. Deconvolut ion, Geophysics Reprint Series 1, Volumes I and II, ed. by G. M. Webster, 1978, Society of Exploration Geophysicists, Tulsa, Oklahoma. Galbraith, J. M., and G. F. Millington, 1979, Low Frequency Recovery in the Inversion of Seismograms: Journal of The CSEG, V.15, No.1, P. 30-39. Kramer, F. S., R. A. Peterson, and W. C. Walter, 1968, Seismic Energy Sources 1968 Handbook: Presented at the 38th Annual ISEG Meeting, Denver, Colorado, 1968 and Privately Published by United Geophysical Corp. Kramer, H. P., and M. V. Mathews, 1956, A Linear Coding for Transmitting a Set of Correlated Signals: IRE Trans. Inform. Theory, V.IT-2, P. 41-46. Lavergne, M., and C. Willm, 1977, Inversion of Seismograms and Pseudovelocity Logs: Geophys. Prosp., • V.25, P. 231-250. Levy, S., and P. K. Fullagar, 1981, Reconstruction of a Sparse Spike Train from a Portion of its Spectrum and Application to High Resolution Deconvolution: Geophysics, V.46, N.9, P. 1235-1243. Lindseth, R. 0., 1979, Synthetic Sonic Logs - A Process for Stratigraphic Interpretation: Geophysics, V.44, P.3-26. Lines, L. R., and R. W. Clayton, 1977, A New Approach to Vibroseis Deconvolution: Geophys. Prosp., V.25, P. 417-433. Lines, L. R., and T. J. Ulrych, 1977, The Old and the New in Seismic Deconvolution and Wavelet Estimation: Geophys. Prosp., V.25, P. 512-540. Mayne, W. H., and R. G. Quay, 1971, Seismic Signatures of Large Airguns: Geophysics, V.36, P. 1162-1173. Oldenburg, D. W., 1981, A Comprehensive Solution To The Linear Deconvolution Problem: Geophys. J. R. Astr. Soc. V.65, P. 331-357. Oldenburg, D. W., S. Levy, and K. P. Whittall, 1981, Wavelet Estimation and Deconvolution: Geophysics, (to be published in V.46, N.11, November, 1981). 105 Otis, R. M., and R. B. Smith, 1977, Homomorphic Deconvolution by Log Spectral Averaging: Geophysics, V.42, N.6, P.1146-1157. Parker, R. L.,1977, Understanding Inverse Theory: Ann. Rev. Earth Plant. Sci., V.5, P. 35-64. Peacock, K. L., 1979, Discrete Operators for Integration: Geophysics, V.44, N.4, P. 722-729. Peterson, R. A., W. R. Fillipone, and F. B. Coker, 1955, The Synthesis of Seismograms from Well Log Data: Geophysics, V.20, P. 516-538. Ready, R. J., and P. A. Wintz, 1973, Information Extraction, SNR Improvement, and Data Compression in Multispectral Imagery: IEEE Transactions On Communications, COM. 21(10), P. 1123-1130. Ricker, N. , 1940, The Form and Nature of Seismic Waves and the Structure of Seismograms: Geophysics, V.5, P. 348-366. Ricker, N., 1953, The Form and Laws of Propogation of Seismic Wavelets: Geophysics, V.18, P. 10-40. Robinson, E. A., 1967, Predictive Decomposition of Time Series with Application to Seismic Exploration: Geophysics, V.32, N.3, P. 418-484. Taylor, H. L., S. C. Banks, and J. F. McCoy, 1979, Deconvolution with the Li Norm: Geophysics, V.44, P. 39-52. Treitel, S., and E. A. Robinson, 1966, The Design of High Resolution Filters: IEEE Transactions, Geoscience Electronics, V.4, P. 25-38. Ulrych, T. J., 1971, Application of Homomorphic Deconvolution to Seismology: V.36, P. 650-660. Ulrych, T. J., 19-72, Maximum Entropy Power Spectrum of Truncated Sinusoids: J. Geophys. Res., V.77 P. 1396-1 400. Ulrych, T. J., and T. N. Bishop, 1975, Maximum Entropy Spectral Analysis and Autoregressive Decomposition: Rev. Geophys. and Space Physics, V.13, P. 183-200. Ulrych, T. J., and R. W. Clayton, 1976, Time Series Modelling and Maximum Entropy: Phys. Earth Planet. Inter., V.12, P. 188-200. 106 Appendix A: Derivation of the Common Trace Consider N traces of seismic data, where each trace s;(t) is composed of a common signal x(t) plus a difference signal nj(t); that is, If the nj (t) represent noise signals, then it is useful to find an x(t) which summarizes the most similar components of the N traces. In general, the common trace may be written as a linear combination of the N traces as follows u The most straight forward procedure is to find that common trace which minimizes the sum of integrated squared differences between x(t) and the Sj (t), i.e., N > z Obviously, this 4 is minimized with respect to x(t) minimization leads to the mean trace 107 where the P| in 2 are all l/N. There is no way to discriminate against an erroneous trace in this procedure because it is simply averaged in with equal importance. It is desireable to find that set of weights fii which emphasize very similar traces and rejects noisy or inconsistent traces. This can be done by formulating the common trace as a linear combination of orthornormal basis functions. Orthogonal Recomposition The s(- (t) span at most N dimensions in Hilbert space. If the s'-(t) are highly correlated, they will have a preferred orientation in that space. This suggests rewiting the common trace as a linear combination of basis functions which are orthonormal and contain one member which minimizes the projection distances between it and the N traces. The orthogonal recomposition of the common trace begins by forming the inner product or covariance matrix 108 Since P i is symmetric, it may De written as r = BA8 T where A-0 is a diagonal matrix of the eigenvalues of P arranged in descending order anc B -Kb,.. - bN is a matrix composed of the eigenvectors \jj of P . The eigenvectors of P define a set of principal orthogonal directions (analogous to the principal stress directions obtained from a stress matrix). The matrix 3 has the following property BST-BTB-I • lo^A.hki -^K J J~ where I is the identity matrix, Xj^is the Kronecker delta and bi'; =B. A set of orthonormal basis functions ^f^(t) may be 109 obtained by projecting the vector of seismic traces onto each eigenvector and normalizing by the corresponding eigenvalue as follows where and .b IT.-CO ^CO it = i 7 a The original traces may be reconstructed from the basis as follows u k=i This may be checked by substituting 6 into 8 as follows M IV U,. L - Si CO 110 The common trace Ms then written as a linear combination of M orthonormal basis- functions M The c^j may be found by minimizing Q-Z.) (xCt)-S,(0) <it Minimizing this, with respect to <^>CV£gives IV /b /- W IV 1 i^biK<fKco^} at 0 Thus N in Substituting 10 and 6 into 9 gives Mr- N If M=N, then 11 reduces to the mean trace If M<N, then 11 can be written tit) - £. Y^j bj • S where
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The recovery of subsurface reflectivity and impedance...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The recovery of subsurface reflectivity and impedance structure from reflection seismograms Scheuer, Tim Ellis 1981-03-26
pdf
Page Metadata
Item Metadata
Title | The recovery of subsurface reflectivity and impedance structure from reflection seismograms |
Creator |
Scheuer, Tim Ellis |
Date Issued | 1981 |
Description | This thesis is concerned with the problem of estimating broadband acoustic impedance from normal incidence reflection seismograms. This topic is covered by following the linear inverse formalisms described by Parker (1977) and Oldenburg (1980). The measured seismogram is modelled as a convolution of subsurface reflectivity with a source wavelet. Then an appraisal of the seismogram is performed to obtain unique bandlimited reflectivity information. This bandlimited reflecitivity information is then utilized in two different construction algorithms which provide a broadband estimate of reflectivity; from which a broadband impedance function may be computed. The first construction method is a maximum entropy method which uses an autoregressive representation of a small portion of the reflectivity spectrum to predict spectral values outside that small portion. The second and most versatile construction method is the linear programming approach of Levy and Fullagar (1981) which utilizes the unique bandlimited spectral information obtained from an appraisal and provides a broadband reflectivity function which has a minimum 1( norm. Both methods have been tested on synthetic and real seismic data and have shown good success at recovering interpretable broadband impedance models. Errors in the data and the uniqueness of constructed reflectivity models play important roles in estimating the impedance function and in assessing its uniqueness. The Karhunen-Loeve transformation is discussed and applied on real data to stabilize the construction results in the presence of noise. The generally accepted idea that low frequency impedance information must be supplied from well log or velocity analyses because of the bandlimited nature of seismic data has been challenged. When accurate, bandlimited reflectivity information can be recovered from the seismic trace, then an interpretable, broadband impedance model may be recovered using the two construction algorithms presented in this thesis. |
Subject |
Seismic reflection method Acoustic impedance |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-03-25 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052871 |
URI | http://hdl.handle.net/2429/22549 |
Degree |
Master of Science - MSc |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1981_A6_7 S38.pdf [ 5.07MB ]
- Metadata
- JSON: 831-1.0052871.json
- JSON-LD: 831-1.0052871-ld.json
- RDF/XML (Pretty): 831-1.0052871-rdf.xml
- RDF/JSON: 831-1.0052871-rdf.json
- Turtle: 831-1.0052871-turtle.txt
- N-Triples: 831-1.0052871-rdf-ntriples.txt
- Original Record: 831-1.0052871-source.json
- Full Text
- 831-1.0052871-fulltext.txt
- Citation
- 831-1.0052871.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052871/manifest