UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Deconvolution and wavelet estimation in exploration seismology Lines, Laurence Richard 1976

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1976_A1 L55.pdf [ 7.09MB ]
Metadata
JSON: 831-1.0052513.json
JSON-LD: 831-1.0052513-ld.json
RDF/XML (Pretty): 831-1.0052513-rdf.xml
RDF/JSON: 831-1.0052513-rdf.json
Turtle: 831-1.0052513-turtle.txt
N-Triples: 831-1.0052513-rdf-ntriples.txt
Original Record: 831-1.0052513-source.json
Full Text
831-1.0052513-fulltext.txt
Citation
831-1.0052513.ris

Full Text

DECONVOLUTION AND WAVELET E5 Tift A I I J N IM EXPLORATION SEISMOLOGY by Laurence R i c h a r d L i n e s B . S c , U n i v e r s i t y of A l b e r t a , 1371 M.Sc. , U n i v e r s i t y o f AlL>erta, 197J THESIS SUBMITTED IN PARTIAL FULFILMENT THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY i n the Department of G e o p h y s i c s and Astronomy We a c c e p t t h i s t h e s i s as c o n f o r m i n g t o r e g u i r e d s t a n d a r d The U n i v e r s i t y Of B r i t i s h Columoia August, 1976 (c) Laurence Richard Lines, 1976 In presenting th is thes is in p a r t i a l fu l f i lment of the requirements for an advanced degree at the Univers i ty of B r i t i s h Columbia, I agree that the L ibrary s h a l l make it f r e e l y ava i lab le for reference and study. I fur ther agree that permission for extensive copying of th is thes is for scho la r ly purposes may be granted by the Head of my Department or by his representa t ives . It is understood that copying or p u b l i c a t i o n of th is thes is for f i n a n c i a l gain shal l not be allowed without my wri t ten permission. Depa rtment The Univers i ty of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 ABSTRACT Seismic deconvolution i s investigated f o r signals generated by explosions and mechanical vibrators. In th i s i n v e s t i g a t i o n , the seismic trace i s modelled as the convolution of the earth's impulse response with a nonstationary wavelet. - By use of several deconvolution techniques, the geologically i n t e r e s t i n g impulse response i s estimated. Wavelet estimation i s often an esse n t i a l part of deconvolution. This study gives the f i r s t extensive comparison of the Wiener-Levinson, Wold-Kolmogorov, and homomorphie deconvolution wavelet estimators i n exploration seismology. In t h i s discussion, the Wold-Kolmogorov method i s shown to be equivalent to Hilbert transform wavelet estimation. wavelet comparisons on noiseless synthetics indicate that homomorphic deconvolution shows considerable promise as a wavelet estimator. However, homomorphic f i l t e r i n g encounters d i f f i c u l t y when additive noise i s present. Hence, cepstral stacking i s used to reduce the noise problem. This thesis also investigates multichannel Wiener f i l t e r s which exhibit i d e a l performance when the wavelets are known. Despite such i d e a l performance on synthetic data, the advantage of multichannel Wiener deconvolution over single channel methods becomes marginal or nonexistant when wavelet estimates are used. A case history of - deconvolving explosion-generated r e f l e c t i o n data i s shown. Comparisons of deconvolutions with well log synthetics are used in order to provide-an inte r p r e t a t i o n of the earth 1s impulse response i n the region of i i i n t e r e s t . Observed differences between synthetic seismoqrams and decon volutions are interpreted i n terras of txie c h a r a c t e r i s t i c s of seismic signals and v e l o c i t y logs. A new approach i s presented for deconvolving seismograms created by vibrator sources. Using c e p s t r a l f i l t e r i n g and spectral d i v i s i o n , the wavelet portion of the trace's amplitude spectrum i s removed. The undetermined portions of the impulse response's Fourier transform are f i l l e d i n by autoregressive prediction. Frequency domain prediction can s U D s t a u t i a l l y increase the time domain resolution of seismic a r r i v a l s . This deconvolution method i s p a r t i c u l a r l y well suitea to vibroseis data because of the phase c h a r a c t e r i s t i c s of the vibroseis wavelet and the known bandlimited spectrum of the vibroseis s i g n a l . Such an approach obviates r e s t r i c t i v e assumptions which are inherent i n conventional approaches to vibroseis deconvolution. The usefulness of t h i s new deconvolution method i s demonstrated for synthetic and r e a l data. i i i TABLE OF CONTENTS Page ABSTRACT i TABLE OF CONTENTS i i i LIST OF FIGURES iv ACKNOWLEDGEMENTS v i i PREFACE v i i i CHAPTER 1. INTRODUCTION TO DECONVOLUTION 1.1 The Impulse Response of a Layered E l a s t i c Earth 1 1.2 Wavelet Model of a Seismic Trace 9 1.3 Attanuation and Dispersion of a Seismic Pulse 13 1.4 Approaches to Deconvolution 14 CHAPTER 2. SEISMIC WAVELET ESTIMATION 2.1 The Wiener-Levinson Double Inverse Method 19 2.2 The Wold-Kolmogorov Factorization Method 23 2.3 Homomorphic Deconvolution 27 2.4 Comparisons of Wavelet Estimators 32 CHAPTER 3. THE DESIGN OF INVERSE FILTERS FOR WAVELETS 3.1 Single Channel Wiener Deconvolution 52 3.2 Multichannel Wiener Deconvolution 61 3.3 Time Adaptive Deconvolution 82 CHAPTER 4.. A CASE HISTORY OF DECONVOLVING REAL SEISMIC DATA 4.1 Deconvolutions Performed 103 4.2 Impulse Response Models Obtained from Synthetic Seismograms 111 4.3 Relating Deconvolutions and Synthetic Seismograms 119 CHAPTER 5. A NEW APPROACH TO VIBROSEIS DECONVOLUTION .5. 1 Introduction' to the Vibroseis Method 130 5.2 A Vibroseis Trace Model 131 5.3 The Vibroseis Deconvolution Problem 136 5.4 The Deconvolution Method 137 5,5' Applications to Synthetic and Actual Data 151 CHAPTER 6. SUMMARY 170 REFERENCES 173 iv LIST OF FIGURES Figure Page 1 P wave propogating at v e r t i c a l incidence 3 through a Goupillaud earth model 2a) Possible displacement waveform and recorded wavelet for a seismic disturbance. 2b) Variation i n the shape of the seismic wavelet as a function of propogation distance. 3 Flow chart outlining approaches to seismic deconvolution and wavelet estimation. 4 A comparison of source wavelet estimates for a synthetic trace. 5 An example demonstrating effects of additive noise on phase curves. 6 A set of f i v e traces from a given shot used i n comparison of wavelet estimators. 7 Use of cepstral averaging on r e a l data of Figure 3. 8 Use of log spectral and raw phase averaging to determine wavelet. 9 Comparison of wavelet obtained by log spectral averaging with minimum phase wavelets. 10 Comparisons of deconvolutions of real data using d i f f e r e n t wavelet estimates. 11 The eff e c t of prewhitening i n the design of single channel channel wiener deconvolution f i l t e r s . 12 A flow diagram of the Wiggins-Bobinson solution for the multichannel Wiener f i l t e r . 13a) A synthetic example exhibiting the ide a l performance of the multichannel wiener deconvolution f i l t e r . 13b) Sum of single channel Wiener deconvolutions using the actual wavelets. 12 12 18 35 38 40 43 47 48 50 60 66 69 70 14a) M u l t i c h a n n e l deconvolution using wavelet e s t i m a t e s . 74 V I 14b) Stacks of single channel deconvolutions designed using wavelet estimates 15 Real data used in a comparison between single channel and multichannel deconvolution. 16 Comparison of single and multichannel deconvolutions on rea l data. 17 Comparison of autocorrelations using the MEM and zero extension methods. 18a) Noisy synthetic trace used i n comparison of Wiener and HEM deconvolution. 18 b),c) A comparison of .Wiener and MEM prediction error f i l t e r i n g . 19 An application of Riley-Burg time adaptive deconvolution. 21b) Data displaying e f f e c t of spherical divergence correction. 22a) Wiener minimum phase deconvolution of the data in Figure 21a) using no prewhitening. 22b) Wiener minimum phase deconvolution of the same data with 10% prewhitening. 76 79 81 91 93 95 99 20 A bandpassed time adaptive MEM deconvolution. 102 21a) Bandpass f i l t e r e d data from an area near the 106 Gulf Coast. 106 108 108 23a) CDP stacks of data i n Figure 21 b). HO 23b) Wiener deconvolution of CDP stacks with no 110 prewhitening. 23c) Wiener deconvolutions of CDP stacks with '\Q% 110 prewhitening. 24 Check shot and CVL vel o c i t y models with 114 synthetic seismograms derived from the velocity curves. 2 5 Comparison of the autocorrelations of the CVL synthetic before and after bandpass f i l t e r i n g . 118 26 I t e r a t i v e deconvolution of common offset stacks for near v e r t i c a l incidence. 121 v i 27 Comparison of a u t o c o r r e l a t i o n a t various stages i n the i t e r a t i v e d e c o n v o l u t i o n process. 123 28 Deconvolution r e s u l t s obtained by Amoco 125 Research f o r the data s e t examined here. 29 Summary of best r e f l e c t i o n p i c k s from the 1 28 de c o n v o l u t i o n s and s y n t h e t i c seismograms obtained from the Gulf Coast data. 30 Example of an i d e a l i z e d model f o r a v i b r o s e i s 1 34 t r a c e . 31 a) Example of a Klauder wavelet, and i t s amplitude 143 spectrum. 31 b ) # c ) Deconvolutions and s p e c t r a r e s u l t i n g 145 from a u t o r e g r e s s i v e p r e d i c t i o n of the F o u r i e r - . transform of the Klauder wavelet. 32 a) Input s p i k e sequence and i t s spectrum. 147 32 b) Mo d i f i e d t r a c e produced by b a n d l i m i t i n g the 148 spectrum. 32 c) Deconvolution which r e s u l t s from p r e d i c t i n g the F o u r i e r transform . 33 d) Deconvolution of a t r a c e s y n t h e s i z e d by using a minimum phase " e a r t h f i l t e r " . 35 a) ,b) ,c) , Input v i b r o s e i s data shown with two subsequent stages of d e c o n v o l u t i o n . 35 d) Comparison of amplitude s p e c t r a at v a r i o u s stages of d e c o n v o l u t i o n . 36 A processed v i b r o s e i s s e c t i o n with subsequent d e c o n v o l u t i o n s . 149 33 a) The amplitude spectrum of a v i b r o s e i s wavelet 1 54 shown with the t r a c e model and the t r a c e spectrum. 33 b) Comparison between the t r a c e spectrum and a 1 5 5 c e p s t r a l estimate of the wavelet spectrum. ' 33 c) Deconvolution of t r a c e i n Fig u r e 33 a) and 1 5 7 r e s u l t i n g spectrum. 158 34 Minimum phase deconvolution of the v i b r o s e i s t r a c e model. 164 166 169 I wish to express my sincere thanks to my supervisor, Professor Tad Ulrych, whose constant encouragement, guidance, and generous assistance mads thi s project possible. I thank Rob Clayton for his invaluable assistance in applying homomorphic deconvolution and autoregressive prediction. I also thank Dr. Ralph Wiggins of Western Geophysical for providing several ideas on time sequence analysis. The assistance of other colleagues at UBC i s also acknowledged. I am grateful to Dr. Ron Clowes for providing a synthetic seismogram program and to Dr. Paul Soiaerville, B i l l Gumming, John B. Davies, and Mat Yedlin for their programming suggestions. I am indebted to Dr. Sven T r e i t e l of Amoco for many helpful discussions on wavelet estimation and deconvolution. I would also l i k e to thank Dr. Larry Wood, Dr. Paul Gutowski, Bob Lucas, Ken Peacock, and B i l l Thompson of Amoco for supplying me with seismic and well log data and data processing suggestions. Vibroseis data and s p e c i a l plotting f a c i l i t i e s sere kindly supplied by George Brinkworth, B i l l Davitt and Dr. Stan Jones of Chevron Standard. F i n a l l y , I would l i k e to thank the National Research Council of Canada for their generous f i n a n c i a l support. This thesis i s e s s e n t i a l l y a description of various approaches to deconvolution. Some of these approaches have been extensively discussed i n the l i t e r a t u r e , but many are applications of cepstral analysis and autoregressive prediction to exploration seismology that are presented for the f i r s t time. In the l a t t e r case, the research of my supervisor, Tad Olrych, had a large influence. To preserve continuity within the text, my o r i g i n a l contributions are emphasized only in the preface. In Chapter 1, I provide models of the seismic trace which are necessary i n l a t e r discussions. Also, a general outline of deconvolution processing techniques i s given as a prelude to the discussions i n Chapters 2 and 3. Homomorphic deconvolution and minimum phase wavelet estimators are compared i n Chapter 2. In that discussion, I show that Robinson's formulation of Wold-Kolmogorov f a c t o r i z a t i o n involves taking a Hilbert transform of the log amplitude spectrum. In estimating wavelets from noisy exploration seismic data, I also applied a technique of log spectral and raw phase averaging which was introduced by Clayton and Wiggins (1976) for the estimation of teleseismic wavelets. In Chapter 3, a close inspection of the multichannel Wiener f i l t e r algorithm showed that certain pathological cases of inverse f i l t e r i n g require prewhitening. Also, Sven T r e i t e l ' s c r i t e r i o n f o r a multichannel f i l t e r length allowed me to f i n d examples of i d e a l multichannel deconvolution. In discussing time adaptive deconvolutions i n Chapter 3, I show that ix differences between wiener and maximum entropy prediction error f i l t e r i n g e s s e n t i a l l y depend on the size of data gates. Following the discussions of inverse f i l t e r i n g , deconvolutions and well log synthetics from a complex geological region are compared. It was encouraging to note that a plausible interpretation could be extracted from such d i f f i c u l t data. In Chapter 5, a new vibroseis deconvolution method developed by myself and Hob Clayton i s presented. The approach exploits c e r t a i n properties of the vibroseis trace which allow cepstral analysis and autoregressive prediction to be used, Deconvolutions of both synthetic and r e a l vibroseis data exhibit the high resolution properties of our technique. 1 CHAPTEB 1 I I I IODOCTION TO DECOHVOLOTIO-N-1*1 i l J B i i i H Besfionse of a Layered E l a s t i c E-arth In recent years, deconvolution has received much attention as a means of enhancing information contained in r e f l e c t i o n seismic recordings. In exploration seismology, deconvolution has been used to resolve geological boundaries within the earth. As an introduction to deconvolution, the r e f l e c t i o n of a plane compressional wave by a series of horizontal rock layers i s considered. The case of v e r t i c a l incidence i s used so that shear waves need not be taken into account. This model of P wave r e f l e c t i o n has been described by several authors including T r e i t e l and Hobinson (1966), Kanasewich (1973), and Claerbout (1976). A diagram of the r e f l e c t i o n seismic method for a Goupillaud earth model i s given i n Figure 1. The modelconsists of a series of horizontal layers and i s mathematically convenient since the t r a v e l time for each layer equals a constant, T/2. As shown i n Figure 1, r e f l e c t i o n s and transmissions of the P wave occur at each interface. In t h i s multilayered medium, a downward propagating pulse causes a seguence of pulse a r r i v a l s to be recorded at the surface. A r r i v a l s resulting from single r e f l e c t i o n s of the seismic pulse are termed primaries, while a r r i v a l s resulting from two or more r e f l e c t i o n s are termed multiples. Figure 1. P wave propagating at v e r t i c a l incidence through a Goupillaud earth model. The rays are drawn at non-vertical incidence to simulate the change of the wave's v e r t i c a l position with time. 0"^ , denote upgoing and downgoing wave amplitudes at the top of each layer. 3 0 . T 2T The r e f l e c t i o n and transmission of a seismic wave at an interface i s described by two boundary conditions. These conditions state that the wave's displacement and the normal stress associated with t h i s displacement are continuous at a boundary. From these conditions, relationships between the incident, r e f l e c t e d and transmitted displacements can be deduced. The r e f l e c t i o n c o e f f i c i e n t for a wave t r a v e l l i n g from layer k to layer k + 1 i s given by: where , A r are the amplitudes of the incident and re f l e c t e d waves, v i s the P wave vel o c i t y , and i s the density of the layer. The transmission c o e f f i c i e n t for the wave i s given by: A 1 ftV«t f>kt, V«„ where A-,: i s the amplitude of the transmitted wave. Claerbout (1976) points out that the reflected wave amplitudes can be deduced from r e f l e c t i o n c o e f f i c i e n t s , and conversely, the r e f l e c t i o n c o e f f i c i e n t s can be deduced from wave amplitudes. To:elaborate on t h i s , the discussion of seismic wave propagation given by-Robinson (1967a) i s followed. Using the model of Figure 1, the z transform i s employed to relate wave amplitudes at the top and bottom of a layer. If the upgoing wave at the top of layer k i s represented by the time series ( u^(0) , u K(1),..., u^ {ti)) , i t s z transform i s given by If z=e , the z transform equals the discrete Fourier transform. The z transform i s a very convenient way of expressing a time series since mu l t i p l i c a t i o n by z gives a unit time delay and d i v i s i o n by z s i g n i f i e s a time advance. If the z transform of the downward propagating wave at the top of layer k'is denoted by D K (z) , a description of the r e f l e c t i o n s and transmissions occuring at the kth interface i s given by the following r e l a t i o n s i n z transform notation (Kanasewich (1973)) . N uH LO) + UKd) Z -f-(1.2) Repeated application of the above formulae for n+1 layers r e s u l t s i n the derivation of a relationship between the i n i t i a l pulse sent into the ground, D0 (z), and the upgoing and downgoing 6 saves i n any l a y e r . (1.4) i s the "communication matrix" f o r the i t h l a y e r . Assuming t h a t l a y e r n+1 i s an i n f i n i t e h a l f - s p a c e and that the i n i t i a l pulse sent i n t o the e a r t h i s a u n i t spike i m p l i e s t h a t U h + ) (z)=0 and D Q {z) = 1-rj 0 (z) . Hence, the values of D^, (z) and 0 o (z) can be determined from the values of the r e f l e c t i o n c o e f f i c i e n t s . C o n s i d e r a t i o n of the Kunetz equation, which r e l a t e s the i n p u t and output energy f l u x e s of t h i s model, allows the r e f l e c t i o n c o e f f i c i e n t s t o be determined from the waves (Claerbout (1976)). . Equating the energy f l u x of the top l a y e r to the f l u x escaping through the lower h a l f space y i e l d s 7 Y, [ (I - U0(^)(l-U.i''i) - U . W UoCV ' j - Y„ D^Wnj*-; (1.5) where Y r e p r e s e n t s the admittance constant f o r the l a y e r s , Equation (1.5) i m p l i e s t h a t the a u t o c o r r e l a t i o n of the escaping wave, D n t, (z) D n +, ( z ~ * ) , i s d i r e c t l y r e l a t e d t o the wave r e f l e c t e d to the s u r f a c e , U c ( z ) . The e s c a p i n g wave, D n +,j- (z) , can be used to determine the r e f l e c t i o n c o e f f i c i e n t s . Robinson (1967a) shows t h i s by use of the f o l l o w i n g polynomial. / V H ) = 1 + AK(l)z +M, AH(k) z * n.6) The polynomial i s r e l a t e d t o the escaping wave by: K U ) ° z " ? f o ±, ( 1 . 7 ) Claerbout (1976) shows that the c o e f f i c i e n t of z i n A K (z) i s egual to the r e f l e c t i o n c o e f f i c i e n t f o r the kth l a y e r . A l s o , A K (z) i s shown to be the z transform of a minimum delay time s e r i e s , which i m p l i e s that a l l the zeros of A k{z) are o u t s i d e the u n i t c i r c l e i n the complex z plane. Equating the c o e f f i c i e n t s of the non-neqative powers of z i n (1.5) a l l o w s the r e f l e c t e d p u l s e s t o be r e l a t e d t o the r e f l e c t i o n c o e f f i c i e n t s by 8 the following matrix agnation, Since the matrix containing the refl e c t e d pulse values i s in Toeplitz form, the algorithm formulated by Levinson (1947) can be used to solve (1.8). In applying the Levinson algorithm, the values of c n are determined at each step i n the recursive solution. As w i l l be seen l a t e r , the equations contained i n (1.8) are i d e n t i c a l to those defining a prediction error f i l t e r , where u Q (t) values are replaced by autocorrelation values, r (t). In f a c t , since the values i n the Toeplitz matrix of (1.8) equal the autocorrelation of the escaping wave multiplied by a constant, A n (z) represents the prediction error f i l t e r for the escaping wave. Hence, the l a s t c o e f f i c i e n t i n a prediction error f i l t e r i s often termed a r e f l e c t i o n c o e f f i c i e n t . Since the autocorrelation of a wave recorded at the surface i s used in the design of the prediction error f i l t e r , t h i s f i l t e r could be related to the r e f l e c t i o n coeffients only under the condition that the autocorrelation of the escaping wave is equal to the autocorrelation of the re f l e c t e d wave. This i s an u n r e a l i s t i c condition of s t a t i o n a r i t y for seismic waves which would not be s a t i s f i e d for the r e a l earth. Hence, i t would not be fe a s i b l e to determine the r e f l e c t i o n c o e f f i c i e n t s by use of prediction error f i l t e r c o e f f i c i e n t s . Although i t i s d i f f i c u l t to d i r e c t l y determine estimates of the r e f l e c t i o n c o e f f i c i e n t s , the u 0(t) time sequence i s of r e a l i n t e r e s t . This sequence, which w i l l be denoted as i ( t ) , i s termed the r e f l e c t i v i t y sequence or the impulse response of a layered earth. If the primary a r r i v a l s i n t h i s sequence can be determined, then a good interpretation of subsurface geology can be obtained. Fortunately, the task of determining primary a r r i v a l s in i ( t ) can be performed by suppressing multiples through use of predictive deconvolution (Peacock and T r e i t e l (1969) ) . _1.2 l a v e l e t Model of a Seismic Trace In the previous description of the earth's impulse response, the idealized case was used i n which the seismic si g n a l was a sequence of delta functions. Because of the nature of seismic sources, the seismic recording system, and the absorption properties of the earth, a r e a l i s t i c model replaces 10 the delta functions fay source generated wavelets. Savelets recorded by seismometers are proportional to the displacement ve l o c i t y of seismic waves and represent the manner in which the earth has f i l t e r e d a source generated s i g n a l . Figure 2a) displays a possible displacement waveform generated by an explosion as well as the recorded wavelet which i s i t s time derivative. A si m i l a r displacement waveform has been presented by Grant and West (1965) to describe the r a d i a l displacement generated by a pressure pulse within a spherical cavity. Many of the present concepts of seismic wavelets have evolved from the work of Ricker (1953). In Bicker's paper, seismic wavelets are considered to be the fundamental units of the seismogram. Using seismic wave theory as well as experiments performed over the Pierre shale i n Colorado, Ricker (1953) formulated fundamental laws governing wavelet propagation. With increase i n propagation time, the wavelet's amplitude decays while i t s width increases. A simulation of t h i s i s given i n Figure 2 (b). A s t a t i s t i c a l d e f i n i t i o n of the wavelet i s given by Robinson (1967b)., The wavelet i s described as a transient time function of f i n i t e energy which represents the predictable part of the seismic trace. F i g u r e 2a). P o s s i b l e d i s p l a c e m e n t waveform and r e c o r d e d wavelet f o r a s e i s m i c wave. F i g u r e 2b) V a r i a t i o n i n t h e shape o f the s e i s m i c w a v e l e t as a f u n c t i o n o f p r o p a g a t i o n d i s t a n c e when t r a v e l l i n g between the s h o t and the r e c o r d i n g geophone. Numbers 1 t o 4 r e l a t e t o f o u r d i f f e r e n t w a v elet shapes which are f u n c t i o n s of p r o p o g a t i o n time a f t e r t h e s h o t e x p l o s i o n . V. REFLECTING HORIZON 13 The seismic trace may be considered as a sum of time delayed wavelets which have been reflected from subsurface interfaces. The wavelet amplitudes are weighted by the values of the impulse response seguence. Hence, a stationary model of the seismic trace, x ( t ) , may be expressed as the convolution of a wavelet, w ( t ) , with the earth's impulse response, i ( t ) . 1*3 Attenuation and Dj^£ersion of a Seismic Pulse The time varying nature of the wavelet i s a r e s u l t of the attenuation and dispersion of seismic waves propagating through the earth. Assuming a linear wave theory and a l i n e a r variation of attenuation with frequency, Futterman (1962) showed that the dispersion of seismic waves i s a conseguence of attenuation. If the complex wavenumbar, k(oj-), i s a function of frequency, the displacement of a plane wave can be expressed as (1.9) 14 [or ix(Rd) = J u(o,*o e ' " di- ( 1 - 1 0 ) where k(w-) = K(ur) + i <^ {U)") . In t h i s case, the Fourier components, U . C O . ^ ' / e have dxfferent phase v e l o c i t i e s . The function K ( w ) determines the phase of each Fourier component and hence governs the wave's dispersion, while C*{6J-) determines the attenuation. Dispersion and attenuation are related by the Kramers - Kronig relationship which e s s e n t i a l l y expresses K ( u f ) as the Hilbert transform of <?C(ur) . The r e l a t i o n s h i p i s derived by using the pri n c i p l e of causality. A s u f f i c i e n t condition for u(R,t) to be causal i s that n(0,o^) has no zeros i n the upper half of the o& plane. This condition on u also implies that 'Futterman - s attenuation laws are minimum phase. Moreover, the attenuation and dispersion laws imply that the seismic wavelet i s nonstationary. 1"H lJE££°l£iiSS to Deconvolution Since wavelets contained i n the seismic trace are time varying and since additive noise i s present, a more r e a l i s t i c but l e s s tractable model of the seismic trace i s given by: 15 X(t) = c(ir) w(t~t} t) + (1.11) Here n (t) represents the additive noise term,-and w(2r,t) i s convolved with i ( t ) . The form w(-7j,t) i s used to indicate that the form of-the wavelet function varies with propagation time (Clarke (1968)). Excluding the noise term, expression (1.11) represents a discrete form of the Booton equation. Although the general solution of (1.11) i s not determined i n seismic processing, steps can be taken to reduce t h i s equation to the simpler convolution form given by (1.9). This can be done by using the seismic data processing techniques which are outlined by Wood and T r e i t e l (1975). Because of the data redundancy i n common depth point (CDP) recordings, normal moveout corrections (NMO) followed by stacking can decrease the noise to sign a l r a t i o i n the seismic trace. Also , the assumption of a stationary wavelet can be j u s t i f i e d i f r e s t r i c t e d time'gates of a trace are processed. Given x (t) as described by (1,9), the deconvolution problem involves estimating the geologically:interesting.impulse , response, i { t ) . To do t h i s , an operator must be designed to remove, or deconvolve, the wavelet from the trace. This i s a two f o l d problem which involves estimating seismic wavelets as "well as designing d i g i t a l f i l t e r s which w i l l shape the wavelet to a spike. • • Various approaches for deconvolving explosion generated 16 seismic data are outlined i n x F i g u r e 3, taken from Lines and Ulrych (1976). These deconvolution techniques are examined and compared i n the following chapters on wavelet estimation and inverse f i l t e r i n g . Figure 3. This flow chart outlines approaches to deconvolution and wavelet estimation that are examined. NMO = normal moveout correction and MEM = maximum entropy method. Input data should have s t a t i c s corrections applied. S O M E A P P R O A C H E S T O D E C O N V O L U T I O N I n p u t : C D P I r e f l e c t i o n ! d a t a O p t i o n a 1 P r e p r o c e s s i n g : b a n d p a s s f i l t e r , N M O , s t a c k , g a i n a p p i i c a t i o n D i v i s i o n of d a t a i n t o g a t e s 1 E s t i m a t i o n o f w a v e l e t s u s i n g : H o m o m o r p h i c d e c o n v o l u t i o n o r W o l d - K o l m o g o r o v f a c t o r i z a t i o n M E M t i m e a d a p t i v e d e c o n v o l u t i o n S i n g l e c h a n n e l W i e n e r o r M E M d e c o n v p l u t i o n 1 S i n g l e c h a n n e l W i e n e r s p i k i n g o f w a v e l e t s O p t i o n to a p p l y N M O a n d s I a c k M u l t i c h a n n e l W i e n e r o r K a l m a n d e e o n v o l n t i o n D E C O N V O L V E D T R A C E S W i e n e r • LevinsonJ d o u b l e i n v e r s e W A V E L E T E S T I M A T E S 19 CHAPTER 2 SEISMIC HAVEIET ESTIMATION Two general approaches for wavelet estimation are examined. One approach uses the assumptions of a minimum phase wavelet and a random impulse response. Various estimation technigues which use these assumptions have been examined by White and O'Brien (1975). A second approach, which does not require the assumptions of the minimum phase approach, i s homomorphic deconvolution, which has been successfully applxed by Ulrych (1971) and Ulrych et a l (1972) to the estimation of teleseismic wavelets. However, apart from the work of Buhl et a l (1974), Otis and Smith (1974), and Buttkus (1975), the applications of t h i s technique to r e f l e c t i o n seismology have not been explored. The Miener-Levinson and wold-Kolmogorov minimum phase estimation techniques are investigated here and are compared with the homomorphic deconvolution method for exploration , seismology problems. 2•I The Siener-LeVinson Double Inverse Method Wiener f i l t e r i n g i s a very common method for performing deconvolution, and i s occasionally used for wavelet estimation,. The theory concerning Wiener f i l t e r s and their design has been thoroughly described by Robinson and T r e i t e l (1967). The Wiener f i l t e r , f (t) , shapes a time sequence, x ( t ) , into an output which approximates the desired output, z ( t ) , in a least squares sense. In the f i l t e r design, the mean squared error given by -t 20 (2 .1) i s minimized by choice of f i l t e r values. That i s . 2B ^ 0 (2.2) for s=0,1,2,...,m. This yields t u t (2.3) Assuming s t a t i o n a r i t y , we can replace t-s by t to get q(s) - £ +(ol] rls-u) (2.4) In the above, r (s) = JJ x (t+s) x (t) i s the autocorrelation of the input, cj (s) = z (t+s) x(t) i s the •t crosscorralation between the input and the desired output, and N i s the length of the input sequence. The normal equations presented in (2,4) are usually written i n matrix form as 21 r r Cm) (I) (O) r(-D ... rC~m) r C O ) : r(0) (2.5) For our purposes, we design a Wiener inverse f i l t e r which t r i e s to shape the wavelet to a spike so that convolution of the f i l t e r with the seismic trace w i l l provide an estimate of the impulse response. If two assumptions are made, a Wiener deconvolution f i l t e r can be designed without e x p l i c i t l y evaluating the source wavelet. If the impulse response i s assumed to be a white noise sequence, the wavelet autocorrelation can be replaced by the trace autocorrelation. This can be seen by writing down the trace autocorrelation i n z transform notation and recognizing that the z transform of the convolution of two functions equals the product of t h e i r z transforms. Using the notation of Chapter 1, W(z) WW) I (z") 22 (2.6) Since I(z)I(z~')=1 for i ( t ) being a white noise seguence, the trace autocorrelation can be used in the f i l t e r design. Also, i f the desired output i s set to be a spike at zero delay, only the f i r s t member of the crosscorrelation vector i n the normal equations w i l l be nonzero. The zero delay spiking position w i l l be optimum only i n cases where the wavelet i s minimum delay. T r e i t e l and Sobinsoh (1966) show that f i l t e r performance w i l l improve i f delays i n the spike position of the desired output are allowed. However, to use a delay of n units i n the desired output requires that we have a good estimate of the f i r s t n points of the seismic wavelet to determine (9 (0) i 9 O) / • - • #9 (a) , 0,. . . ,0) i n the normal equations. To obtain an actual wavelet estimate, another Wiener f i l t e r i s applied to invert the o r i g i n a l minimum delay inverse f i l t e r . Consequently, t h i s wavelet estimation technique i s termed the Wiener-Levinson double inverse method. In using this approach, there exists the problem of using the r e s t r i c t i v e minimum phase and randomness assumptions as well as the problem of determining a suitable length for the Wiener f i l t e r . The Wiener inverse f i l t e r i s designed d i r e c t l y from the autocorrelation of the trace i n the f i r s t step of the wavelet estimation process. Therefore, the Wiener-Levinson double inverse method i s r a r e l y used i n practice, since i t i s the inverse f i l t e r , not the wavelet, which i s usually desired. 23 2.2 lold-Kolraogorov F a c t o r i z a t i o n The Wold-Kolmogorov f a c t o r i z a t i o n method, M h i c h i s d e s c r i b e d i n d e t a i l by Robinson (1967b), determines the minimum phase spectrum f o r a given amplitude spectrum. Given an amplitude spectrum f o r a wavelet of length n+1, there are 2n wavelets with d i f f e r e n t phase s p e c t r a which w i l l have t h i s amplitude spectrum. Only one of these wavelets i s minimum phase, and i t s spectrum may be deduced by the f o l l o w i n g approach. T h i s approach of determining the minimum phase wavelet i s e q u i v a l e n t to f a c t o r i n g the z transform of the auto-c o r r e l a t i o n i n t o minimum phase binomial f a c t o r s . F i r s t of a l l the l o g of the wavelet's F o u r i e r t r a n s f o r m , a ( o/) , i s giv e n by spectrum. For a c a u s a l wavelet, l o g (W ( ur) ) i s giv e n by a F o u r i e r expansion: (2.7) where |W(u>-) | i s the amplitude spectrum and O^iuf) i s the phase (2.8) As pointed out by Robinson, (2.8) i s d e r i v e d under the c o n d i t i o n t h a t l o g ( w (z)) be a n a l y t i c i n s i d e the u n i t c i r c l e i n 24 the complex plane. This requires that the wavelet's z transform, I (z), have a l l i t s roots outside the unit c i r c l e which i s the condition of minimum phase. The Fourier series for the even function, loglW(o-) I, i s given by oO /talWMl = oC(0) + 2 g COS (St (2.9) where Tr Jl-t M oc(±) = T T J / o j l WCA)| COS (2.10) Comparison of (2 .8) and (2.9) shows that # (t) = 2 ^ (t) for t> 1 , and # { 0 ) = c < ( 0 ) . Also, comparison of (2.7) and (2.8) then shows that o O &(<*) = ~ 2 Z ^ M (2 . 11 ) Hence, the minimum phase spectrum, &(kr), i s given in terms of the amplitude spectrum, J W ( ur) | , by the following r e l a t i o n . 25 IT &(ur) = ~2 J T sin est f cosM foSIW(JOldfl *2*12* 1 7 t=> v / 0 It can readily be shown that the above expression for O'(u-) derived by Robinson (1967b) i s the Hilbert transform of the log amplitude spectrum. F i r s t of a l l , since the Hilbert transform, 0 (&) , i s e s s e n t i a l l y the convolution of — with h j lw((*r/} , taking the Fourier transform of (f)(ur) y i e l d s the following, F(<J>) = - fe" ' M I'} f ^ M ( 2 . i 3 , Since i s an even function defined between —IT and TT , (2.13) can be rewritten as: TT F(0) = 2 i j C O S A * loj\VJLJl)\dJl > S3r>Lt) (2.14) Inverse Fourier transformation gives the o r i g i n a l value of the Hilbert transform. The transform i s discrete because we are dealing with d i g i t a l values in the time domain. Hence, (f)(<*) i s given by: 1*9 26 (2.15) Since sgn (t) i s an odd function of t, the above may be written i n the same form as (2.12). Hence, the Wold -Kolmogorov f a c t o r i z a t i o n method described by Robinson (1967b) i s often termed the H i l b e r t transform method. The wavelets obtained by use of the Hilbert transform are often quite s i m i l a r to those obtained by using the Wiener -Levinson double inverse method since both methods are determined by using an autocorrelation estimate and the minimum phase assumption. Galbraith (1971) points out that the methods would be equivalent i f an i n f i n i t e number of f i l t e r points could be used i n the Wiener f i l t e r . The analysis given-here i m p l i c i t l y uses the d e f i n i t i o n of a wavelet as a transient time sequence with f i n i t e energy, following Robinson (1967b) . In applyinq Wiener-Levinson wavelet estimation t h i s d e f i n i t i o n requires the choice of appropriate f i l t e r lengths, whereas the wavelet obtained by using the Hi l b e r t transform must be •truncated at some suitable point. The minimum phase wavelet estimators examined here have used a s t a t i s t i c a l d e f i n i t i o n of the impulse response as being a white noise sequence. Recently, Robinson (1975b) has used a dynamic modal for the impulse response fay applying the wave theory presented i n Chapter 1. 2.3 S2!°!9J.EhiLc: Deconvolution Homomorphic deconvolution , which was introduced by Oppenheim (1965), does not require the r e s t r i c t i v e conditions imposed by the minimum phase methods of wavelet estimation. In usinq homomorphic deconvolution, a c h a r a c t e r i s t i c transformation i s applied to a time sequence x ( t ) , i n order to y i e l d the complex cepstrum. The complex cepstrum of as time sequence i s defined as the inverse Fourier transform of the logarithm of the time sequence's Fourier transform. For a d i g i t a l l y sampled seguence, x ( t ) , the complex cepstrum, x(t) i s given by: The terms "quefrency" and , ,cepstrum , , were introduced by Bogert et a l (1963). The term "quefrency" describes the "frequencies" or r e p e t i t i o n rates of a frequency domain sequence and the complex cepstrum i s e s s e n t i a l l y a measure of the quefrency content of the log cf a sequence's Fourier transform. Calculation of the complex cepstrum transforms sequences 28 which are convolved i n the time domain into sequences which are additive i n a new domain termed the quefrency domain. The complex cepstrum of the convolutional model for the seismic trace w i l l be considered. The Fourier transform of x(t) = V w(s) i ( t - s ) yields X(^) ~ U/(w) J. Taking the logarithm gives laj X((JS) - / o j W(w) + /o j 1(UT) (2.17) The inverse Fourier transform of (2.17) then produces x(-t) = wit) ~h tit) (2.18) where w (t) and i (t) are the cepstra of the wavelet and the impulse response as defined by (2.16). The computation of the complex cepstrum involves certain essential steps which are outlined by Schafer (1969), Ulrych (1971), and S t o f f a et a l (1974), F i r s t , the expression for the log of the Fourier transform of x (t) i s considered. 29 (2.19) To make the phase spectrum, 0*(ux) , an analytic continuous function of freguency requires that the raw phase values computed by the Arctan function be unwrapped by adding multiples of 2TT. The l i n e a r component should also be removed from the phase since i t may obscure the intere s t i n g information i n the cepstrum. Removal of the li n e a r phase component, bur , from the phase, , simply s h i f t s the o r i g i n time of the trace by b units. Also, the mean of log X(or) i s usually removed. This operation sets x(0) to zero and i s equivalent to multiplying the o r i g i n a l time sequence by a constant exp(-a) where a i s the mean of log X(cj-). Hence, the t=0 value in the complex cepstrum i s only affected by the scaling in x (t) . As seen by (2.18), sequences which are convolved i n the time domain are additive i n the guefrency domain and an estimate of the wavelet cepstrum can be obtained by l i f t e r i n g ( i . e . l i n e a r f i l t e r i n g i n the guefrency domain). The wavelet cepstrum i s often estimated by low cut l i f t e r i n g . That i s , the wavelet cepstrum i s set egual to complex cepstrum values f o r J t | less than some cutoff value. The separation of the wavelet cepstrum and the impulse response cepstrum i s p a r t i c u l a r l y convenient when the impulse-response i s minimum phase. This can be seen by considering the form of i (t) . 30 The log of the Fourier transform of i ( t ) i s given by Lit) e it) e j (2.20) If \2L e J < 1, a Laurent expansion for log(I(ur)) can be used. Thus, / o , "t) e — ; = r ( 2 i « » e <2 .2 i , A The value of i ( t ) i s given by r2lT OO fc+/ o " l f r < » ^ T (2. 22) or (2 . 23) -IMS Hence, i f 1 ^  ^ s ) e | i s much less than one, the largest spikes i n the impulse response w i l l correspond to those i n the o r i g i n a l sequence, i (t) . This condition on I(cj) ensures minimum phase and can be imposed by exponential weighting. 31 Exponential weighting i s performed fay multiplying terms i n t sequence x{t) by °< , where 0<°<<1 .0 . This i s equivalent to replacing z=e x.n the z transform, X (z), by z=e J , and e f f e c t i v e l y moves the zeros of X (z) outward from the o r i g i n . Exponential weighting makes the impulse response appear to be minimum phase and x{t)=x(t). Therefore, i f the f i r s t and second a r r i v a l s contained i n a minimum phase impulse response are at t=0 and t=s respectively, the f i r s t spikes in the cepstrum of the impulse response are at t=0 and t=s also. An i d e a l deconvolution can then be performed i f the wavelet cepstrum i s l o c a l i z e d within t=s, since in such cases a complete separation of the wavelet cepstrum and the cepstrum of the impulse response can be made. As emphasized by Stoffa et a l (1974), exponential weighting can also be used to eliminate cepstral a l i a s i n g caused by the nonlinear operations of the logarithm, absolute value and the arctangent. After the cepstrum has bean calculated and l i f t e r e d to give an estimate of the wavelet cepstrum, the inverse c h a r a c t e r i s t i c transformation i s applied to y i e l d the source wavelet estimate. This inverse c h a r a c t e r i s t i c transformation i s performed by taking the Fourier transform, the exponential and the inverse Fourier transform. I f we denote the wavelet cepstrum by w(t), the wavelet estimate i s given by w ( t ) ^ { iu r t + 2 e W s ) } S=-oa 32 (2.24) Although the cepstrum of the impulse response could be transformed directly- to obtain a deconvolved trace, i t i s often d i f f i c u l t to obtain the actual cepstrum of the impulse response because the impulse response spikes contained in the cepstrum may be obscured by high guefrency noise. On the other hand, the cepstrum of the wavelet i s often l o c a l i z e d i n the low guefrency part of the complex cepstrum and can be obtained by applying a low cut l i f t e r to the complex cepstrum of the trace.. The wavelet obtained from t h i s cepstrum can then be used i n the design of an inverse f i l t e r to deconvolve the trace. 2.4 Comparisons of Wavelet Estimators The performance of the previously mentioned wavelet estimation techniques i s compared by using r e a l and synthetic data. To the author's knowledge, no other comparison of these wavelet estimators has appeared i n exploration seismology l i t e r a t u r e . Figure 4 compares wavelet estimation methods by using a synthetic trace. The trace i s formed by convolving a spike sequence with a wavelet. This wavelet, which was suggested to the author by G.L, Cumming, was synthesized by using a z transform of the form (-1,10 + z ) * (1.75+ z)'"4 and removing the average, The re s u l t i n g time sequence simulates a possible explosion generated wavelet, similar to those found by Ricker (1953). Because of removal of the average, the wavelet i s not quite minimum phase. Also, since the impulse response i s not a white noise sequence, the autocorrelation used i n the Wold-Kolmogorov and Wiener-Levinson wavelet estimates consists of the trace autocorrelation multiplied by a Parzen window of length 20. Although the windowed trace autocorrelation i s a good estimate of the source wavelet autocorrelation, the mini-mum phase assumption gives both wavelet estimates a d i f f e r e n t character from that of the actual wavelet. The homomorphic deconvolution estimate of the wavelet was obtained by low cut l i f t e r i n g the cepstrum. In this case the complex cepstrum of the trace was calculated after using an exponential weighting of 0.98 on the trace. As seen from the cepstra of Figure 4, the actual source wavelet cepstrum cannot be completely separated from the complex cepstrum of the trace by a low cut l i f t e r . Hence, the t h i r d lobe of the homomorphic estimate i s not quite the same as that of the actual wavelet although the f i r s t two lobes of the wavelet are v i r t u a l l y duplicated. Such tests on noiseless synthetics provide encouragement for the use of cepstral methods i n wavelet estimation. 34 F i g u r e 4 . A comparison of source wavelet e s t i m a t e s found by the Wold-Kolmogorov f a c t o r i z a t i o n method (WKFACT), the Wiener-Levinson double i n v e r s e method, and the method of homomorphic d e c o n v o l u t i o n . 35 COMPARISON OF WAVELET ESTIMATORS SOURCE WAVELET IMPULSE RESPONSE i — i 1 1 1 1 0.0 0.2 0.4 0.6 O.B 1.0 TIME SOURCE TRACE WKFflCT WIENER-LEVINSON AUTOCORRELATION . AUTOCORRELATION ESTIMATE DOUBLE INVERSE SOURCE CEPSTRUM TRACE CEPSTRUM HOMOMORPHIC ESTIMATE For the case of zero additive noise, the performance of homomorphic deconvolution as a wavelet estimator depends on how well the complex cepstrum of the wavelet can be separated from the complex cepstrum of the trace by low cut l i t t e r i n g . S i m i l a r l y , the performance of the minimum phase approaches depends cn how well the wavelet autocorrelation can be estimated by windowing the trace autocorrelation. Additive noise can cause major problems when using homomorphic deconvolution because i t i s d i f f i c u l t to determine i t s e ffect on the complex cepstrum. This i s because cepstral calculations involve nonlinear transformations. Furthermore, additive noise i n the seismic trace can cause severe i n s t a b i l i t i e s when unwrapping the phase curve. This problem i s demonstrated i n Figure 5. The noiseless trace shown i n the figure has a maximum amplitude of 1.0. Different r e a l i z a t i o n s of "white noise" with a standard deviation of 0.3 were added to t h i s trace. After the usual procedures of phase unwrapping and removal of the l i n e a r phase component, the phase curves show marked differences. Figure 5. These examples demonstrate that d i f f e r e n t noise r e a l i z a t i o n s added to a trace can produce s i g n i f i c a n t differences in the unwrapped phase curve, "white noise" r e a l i z a t i o n s , with standard deviations of 0.3, were added to a trace with a maximum amplitude of 1.0. 38 TRACE i.e NOISY TRACE WT=1.00 NOISY TRACE WT=1.00 NOISY TRACE WT=1.00 PHASE Clayton and Wiggins (1976) have investigated t h i s problem and have shown that deviations i n the phase caused by additive noise become severe when signal amplitudes are small. They also demonstrate that 2 Tf errors in the phase unwrapping can e a s i l y occur when the amplitude spectrum values approach zero. Therefore, from the work of Clayton and Wiggins (1976) and the author's experience, i t i s suggested that data redundancy be used to suppress the effects of noise on the phase curves when using homomorphic deconvolution. Recently, Otis and Smith (1974) and Clayton and Wiggins (1976) have devised methods of log spectral averaging which decrease the e f f e c t s of noise on the phase curve. These methods of c e p s t r a l averaging hopefully w i l l enhance the comp-lex cepstrum of the wavelet. Variations of cepstral stacking methods are applied to the windowed portion of the f i v e traces shown i n Figure 6. Traces from the same shot are used i n order to maintain redundancy i n wavelet shape. This example i s p a r t i c u l a r l y i n t e r e s t i n g because the f i f t h trace contains considerably more additive noise than the other traces. 40 Figure 6 . A set of f i v e traces which are used i n the comparison of wavelet estimators. The windowed portion of the trace was used in the cal c u l a t i o n of wavelets. 41 For c e p s t r a l computations, the windowed portions of the traces of Figure 6 were exponentially weighted using a weighting factor of 0.96. Figure 7 shows the phase curves of the exponentially weighted traces, the wavelets obtained by homomorphic deconvolution, the sum of the complex cepstra, and the wavelet obtained by low cut l i f t e r i n g of the c e p s t r a l sum . Due to the e f f e c t s of noise, the phase curve of the f i f t h trace i s considerably d i f f e r e n t from those of the other traces. Averaging the cepstra corresponds to averaging the log amplitude and phase spectra before taking the inverse Fourier transform. This can be seen from the expression for the average of the cepstra for N channels which i s given by ir W M I 2 - 2 5 ) 2.TT - TT where X K (ur- ) denotes the Fourier transform of x (t) . Figure 7. The use of cepstral averaging in wavelet estimation i s shown using the traces of Figure 6. This figure shows the phase curves of these traces, the wavelets obtained by l i f t e r i n g the cepstra of the traces, the cepstral average, and the wavelet obtained by applying a low cut l i f t e r of width 60 msec to the ceps t r a l average. Scale for cepstral average and average wavelet i s i n seconds. 43 PHASE CURVES -125 hz H 125 hz HOMOMORPHIC DECONVOLUTION WAVELETS fa F R E Q U E N C Y fa-fa V*r— 200 msec C E P S T R A L A V E R A G E r -0.4 0.0 0.4 S O U R C E W A V E L E T 0.2 44 Rewriting (2.25) in the following form shows that averaging the log amplitude spectra, i s equivalent to finding the log of the geometric mean of the amplitude spectra. (2.26) When finding the cepstral average, i t i s obviously less expensive to average the log amplitude and phase spectra before taking the inverse Fourier transform than to average the cepstra themselves. In taking a cepstral average, we are e s s e n t i a l l y using an average of the unwrapped phase. However, Clayton and Wiggins (1976) have proposed a method of averaging the raw phase values before unwrapping. This produces a phase curve with greater s t a b i l i t y than that obtained through averaging the unwrapped phase. Figure 8(a) displays the log amplitude spectrum with i t s mean removed as well as an unwrapping of the averaged raw phase values. The corresponding complex cepstrum and the low cut l i f t e r e d cepstrum used i n wavelet computation are also shown. Figure 8(b) shows that setting the log amplitude and phase spectra to zero when the signal's amplitude spectrum becomes less than a certain value has a negli g i b l e e f f e c t on the cepstrum. In Figure 9, the wavelet obtained by averaging the raw phase and log amplitude spectra i s compared with wavelets obtained by use of the Hilbert transform. From Figures 7 and 9, 45 i t can be seen that wavelets estimated by cepstral averaging methods are d e f i n i t e l y mixed phase, as compared to the minimum-phase wavelets obtained by use of the Hilbert transform. I t should be noted that the minimum phase wavelets shown in Figure 9 are delayed one point since the f i r s t sample of the wavelet i s not zero. This property, which i s common to a l l minimum phase wavelets, i s not physically r e a l i s t i c . Figure 10 displays deconvolutions which result from the application of Wiener f i l t e r s designed to spike the estimated wavelets. Although the obvious r e f l e c t i o n s are shown on a l l three sets of deconvolved traces, deconvolutions using the wavelet estimated from averages of the raw phase and log amplitude spectra give s l i g h t l y better resolution of overlapping re f l e c t e d pulses, especially between times of 1.5 and 1.7 seconds on the near v e r t i c a l r e f l e c t i o n s . Figure 8. a) Estimation of the wavelet using log spectral averaging for the traces of Figure 6. The average of the log amplitude spectra i s shown as well as the phase curve obtained by unwrapping the averages of the raw phase curve. An exponential weighting of 0.98 was used. The phase and log amplitude averages were used to compute the cepstrum. The l i t t e r e d version of the cepstrum (also shown) i s used to produce the wavelet estimate. Figure 8b) shows that setting the phase and amplitude spectra to zero after a certain freguency produces a very s i m i l a r cepstra to that obtained by using the entire frequency range. In b), the second cepstrum shown i s a l i t t e r e d version of the f i r s t cepstrum. Differences i n the wavelets obtained by using the spectra i n a) and b) could not be determined by visual inspection. The wavelet obtained i s shown i n Figure 9. a) LOG SPECTRAL" RVERRGING LOG RMP SPEC PHRSE CURVE o. •r. a o _ l o ft 0 0 6 ^ \ j 5 1 " UJ 125 hz COMPLEX CEPSTRUM TIME 0.4 COMPLEX CEPSTRUM 0.4 b> LOG SPECTRAL RVERRGING LOG RMP SPEC PHASE CURVE U J Vi <x X a. 125hz i 125 hz COMPLEX CEPSTRUM COMPLEX CEPSTRUM -0.4 0.4 HILBERT TRANSFORM WAVELETS i — j 200 msec Figure 9. Comparison 'of the wavelet obtained by log spectral averaging with wavelets obtained by use of the Hilbe transform method to the traces shown i n Figure 6. F i g u r e 10, Comparison of deconvolutions of the t r a c e s of F i g u r e 6 using a Wiener f i l t e r of l e n g t h 80 msec to supply an i n v e r s e f i l t e r f o r ; a) the wavelets obtained fay Wold-Kolmogorov f a c t o r i z a t i o n b) the wavelet obtained by c e p s t r a l averaging c) the wavelet obtained by l o g s p e c t r a l averaging, \ DECON WITH WKFRCT WAVELETS 50 — r 1 1 1 1 2.j a.< s.t J . i »•» ero Tt »~i 7* 1.4 i.t TIME. 2.0 k ) DECON WITH CEP RV WAVELETS —4' .5 H E K K £ n " " " " " DECON WITH LOG SPEC RV WRVELET 51 In discussing the r e l a t i v e merits of the estimation schemes, i t should be mentioned that despite r e s t r i c t i v e assumptions used i n i t s design, the Hilbert transform method requires only an estimate of the wavelet length i n order to produce a reasonable wavelet. On the other hand, the homomorphic deconvolution method may require more "parameter juggling" since one must choose an appropriate exponential weighting and a suitable low cut l i f t e r . Also, homomorphic deconvolution results are more dependent on the choice of data gates than the Hi l b e r t transform method, and the.additive noise presents a problem. These factors generally mean that homomorphic wavelet estimates are more expensive to obtain but can be more sensitive to the data when the noise l e v e l i s low and the source generated wavelet i s not minimum phase. Having discussed d i f f e r e n t approaches f o r estimating wavelets, attention i s now turned to the second half of the deconvolution problem, the design of inverse f i l t e r s which shape these wavelets to a spike. 52 CHAPTER 3 THE DESIGN OP INVERSE FILTERS FOR M&VELETS 3.1, Sinqle Channel Wiener Deconvolutioi]. Wiener f i l t e r i n g represents one of the most common approaches to deconvclution i n exploration geophysics. The equations defining the single channel Wiener f i l t e r vere outlined previously. The design and application of these f i l t e r s i s now discussed. Osing a wavelet estimate as input and choosing a spike as a desired output allows the Wiener deconvolution f i l t e r to be designed. Although shaping f i l t e r s are not discussed here, T r e i t e l and Robinson (1966) and DeVoogd (1974) have shown that Wiener f i l t e r s which shape the wavelet to some form other than a spike may also prove useful for deconvolution. Use of shaping f i l t e r s also presents the problem of deciding on a suitable desired output. Spiking f i l t e r s are considered not only because of their s i m p l i c i t y , but also because of t h e i r i n t e r e s t i n g relationship to autoregressive models and prediction error f i l t e r s . In assuming minimum phase, the Wiener spiking f i l t e r can be obtained d i r e c t l y from estimates of the wavelet autocorrelation. Peacock and T r e i t e l (1969) relate the minimum pnase Wiener inverse f i l t e r to a prediction error f i l t e r . Prediction error f i l t e r i n g i s examined here since i t provides the basis for a l a t e r discussion of the Burq algorithm. 53 I f t he u n i t d i s t a n c e p r e d i c t i o n f i l t e r o f l e n g t h m i s g i v e n by the seguence { ^ ( l ) , o ( m ( 2 ) , . » . , o£ m(m) ) , a p r e d i c t e d v a l u e of t h e s e g u e n c e , x ( t ) , i s deno ted by : - t - l I f t h e p r e d i c t i o n e r r o r s a r e r e g a r d e d a s a w h i t e n o i s e s e q u e n c e , x ( t ) , w r i t t e n as x (t)-= (1) x ( t -1) + <^ n > {2)x( t -2) + . . . + o^ m (in) x (t-m) + a / > 1 ( t ) , can be d e s c r i b e d a s an a u t o r e g r e s s i v e p r o c e s s of o r d e r m. The p r e d i c t i o n e r r o r seguence a s s o c i a t e d w i t h t h i s p r e d i c t i o n f i l t e r i s g i v e n by: am(i) = X U ) ~ £ <*«(•<:-*) ><M (3.2) o r (3.3) where ) f m { t ) i s t he m p o i n t p r e d i c t i o n e r r o r f i l t e r d e f i n e d by t h e f o l l o w i n g . 54 (3.4) Y j t ) = -°*M c-t) ) t> o By minimizing the prediction error power, Peacock and T r e i t e l (1969) have shown that t h i s prediction error f i l t e r s a t i s f i e s the eguations defining a minimum phase Wiener inverse f i l t e r to within a scale factor. That i s , where P^ i s the prediction error power associated with a f i l t e r of length m. The algorithm for solving the normal eguations presented in (1 . 8 ) , (2.5) and (3.5) was given by Levinson (1947). This algorithm i s a recursive solution which exploit s the Toeplitz form of the autocorrelation matrix. In t h i s scheme, the prediction error f i l t e r of length m+1 i s formed by combining the prediction error and hindsight error f i l t e r s of length m through use of the f i l t e r c o e f f i c i e n t o^m(m) . The f i l t e r i s chosen such that: (3.5) and where Therefore, °^_(m) i s chosen so that 56 (3.10) The recursive solution for the prediction f i l t e r i s then given by As seen from the previous discussion of v e r t i c a l l y r e f l e c t e d waves, oC (m) can , under certain conditions , be considered as a r e f l e c t i o n c o e f f i c i e n t . That i s , the l a s t f i l t e r c o e f f i c i e n t calculated by each step i n the Levinson recursion also has a physical interpretation. Before digressing from the discussion of prediction error f i l t e r s , i t should be pointed out that an important class of prediction error f i l t e r s e xists for which the prediction distance i s greater than one sample. As pointed out by Peacock and T r e i t e l (1969), such f i l t e r s can be very useful in eliminating water bottom multiples on marine seismic records. Having discussed the solution of the single channel Wiener f i l t e r , the question of determining a suitable f i l t e r length i s considered. f If a wavelet i s minimum phase, the z transform of i t s exact inverse f i l t e r i s given by: (3.11) • t = o 57 (3,12) In other words, a causal f i l t e r of i n f i n i t e length i s reguired i n the single channel case to perform an i d e a l deconvolution. Despite the fact that performance of the inverse f i l t e r should improve as i t s length increases, expense and roundoff errors are computational factors which r e s t r i c t f i l t e r length. Although the choice of a suitable length for the single channel Wiener f i l t e r i s subjective, i t i s i n s t r u c t i v e to monitor the normalized mean square error, E(m). This error function equals the previous expression given in (2.1) divided by a normalizing f a c t o r , r ^ 2 ( 0 ) , which i s the zero lag value of the autocorrelation of the desired output. For a f i l t e r of length m, E(m) reduces to The use of t h i s error function i n designing Wiener f i l t e r s i s described by T r e i t e l and Robinson (1966) and Wiggins (1968). Often an appropriate length for the Wiener inverse f x l t e r l i e s somewhere between the wavelet length and twice the wavelet length. In the design of inverse f i l t e r s , i t i s useful to apply prewhitening. Prewhitening involves adding a constant to the (3. 13) 58 spectrum of the wavelet. This eliminates zeros i n the wavelet spectrum so that the spectrum's reciprocal (the spectrum of the inverse f i l t e r ) i s well defined at a l l frequencies. In frequency ranges where the wavelet's spectral values are small, inverse f i l t e r s enhance noisy portions of the trace's spectrum. In Wiener f i l t e r design, prewhitening can also be accomplished by multiplying the wavelet autocorrelation's zero lag value by some number greater than 1 .0 . This i s equivalent to adding white noise to the wavelet. A demonstration of the e f f e c t s of prewhitening i s shown in Figure 1 1 . This figure displays a set of wavelet estimates from CDP stacks of data, as well as the Wiener inverse f i l t e r s which have been designed to shape these wavelets to a spike. F i l t e r s designed with no prewhitening are compared to those designed using 10% prewhitening. The output which results from the convolution of the f i l t e r s with the wavelets i s also shown. The f i l t e r e d wavelets obtained using no prewhitening give outputs more highly resolved than those obtained using prewhitening. However, these highly resolved f i l t e r outputs also contain high frequency noise. Furthermore, deconvolutions obtained from f i l t e r s designed with no prewhitening can display undesirable edge effects at the end of a time gate. Hence, in choosing a l e v e l of prewhitening, a tradeoff between resolution and the amount of high frequency noise must be considered. When a wavelet's spectrum i s limited to the lower frequencies, bandpass f i l t e r i n g i s obviously an alternative method of suppressing high frequency noise. Figure 11. The affect of prewhitening in the design of Wiener f i l t e r s . a) The wavelets shown were estimated from actual data. Inverse f i l t e r s designed with d i f f e r e n t l e v e l s of prewhitening are displayed. No prewhitening i s used in fa) and "\0% prewhitening i n c) . Outputs resu l t i n g from the convolutions of the inverse f i l t e r s with the wavelets are shown. a) S O U R C E W R V F L E T E S T I M A T E S b) W I E N E R F I L T E R S c) W I E N E R F I L T E R S F I L T E R O U T P U T 1*-+— ^ j u > - ^ V ~ l ~ A j ^ A - ^ - i ^ -F I L T E R O U T P U T 61 3. 2 Multichannel Wiener Deconvolution Thus f a r , only single channel deconvolution of seismic traces has been considered. Since modern seismic data have a great deal of redundancy due to multifold coverage, multichannel f i l t e r i n g can be used to exploit the information contained i n CDP recordings. The basic assumption in multichannel f i l t e r i n g i s that the output signals from the f i l t e r i n g system are related to information contained in more than one input s i g n a l . Hence, the i t h channel of the f i l t e r e d output i s given by Consequently, multichannel f i l t e r i n g involves operations with matrices instead of the scalar operations used i n single channel f i l t e r i n g . The following discussion , which i s very s i m i l a r to that given by Lines (1974), outlines the application of Wiener multichannel deconvolution. The extension of the Levinson algorithm to i t s multichannel form was o r i g i n a l l y qiven by Siqgins and Robinson (1965). Several Fortran programs for multichannel data analysis are given by Robinson (1967a), and an excellent review of the mathematics of multichannel f i l t e r i n g i s presented by T r e i t e l (1970). In that review, T r e i t e l derives the equations defining the multichannel f i l t e r . The f i l t e r matrices, f j t ) s a t i s f y the following system of normal equations for a f i l t e r of k J (3.14) 62 l e n g t h nu 1 ( f ( 0 ) f(/)... ffrOj j r(0) r ( 0 ... /g(o) 3O)... aft*); (3 .15) The matrix f (t) contains elements f,j (t) which operate on the j t h input channel to give a contribution to the i t h output channel. The matrices r(t) and g (t) represent the autocorrelation and crosscorrelation matrices respectively. Using the assumption of ergodicity, the values of r ( t ) and g (t) are given by the following r e l a t i o n s . r( t ) J( t ) = 1 s s 5 y 3 . i 7 ) 5 5" S I • *• : Here (x , (t) ,x 2 (t) ,..., (t) ) represents the input time sequences on k, channels and (Y( (t) , Y^  (t) ,..., Y^  (t) ) represents the desired output sequences on each of the J[ output channels. In the case of multichannel deconvolution, k i s usually the number of traces i n a CDP andj^ i s set equal to one since we desire a single impulse response for the redundant data. The f a c t that r(t) i s in the form of a Toeplitz matrix allowed Wiggins and Robinson (1965) to obtain a recursive solution to the equations defining a multichannel f i l t e r . A flow diagram of t h i s solution i s given i n Figure 12. This flow diagram corrects misprints in the previous diagram given in Wiggins and Robinson (1965) by following the description given by Robinson (1967a). One of the misprints was i n the description of the algorithm i n i t i a l i z a t i o n . This description f a i l e d to mention that the matrix _r (0) was inverted. By using the multichannel solution outlined in Figure 12, an inverse f i l t e r may be designed to convert the source wavelet to a spike. Multichannel convolution of f.(t) with the traces from a given CDP then gives an estimate of the impulse response. However, 64 the source wavelet estimates which are the multichannel inputs to t h i s program must not be i d e n t i c a l . As seen i n Figure 1 2 , i n i t i a l i z a t i o n of the recursive process for the f i l t e r design reguires that f(0)=g(0)r_ ( 0 ) , where r (0) i s the inverse of the autocorrelation matrix at zero lag. If a l l inputs to the f i l t e r are i d e n t i c a l then r_(0) i s a matrix with egual elements and i t s inverse does not e x i s t . This problem of the zero lag autocorr-e l a t i o n matrix being singular implies that prewhitening may be necessary i n the multichannel case. Prewhitening can e a s i l y be accomplished by multiplying the diagonals of r (0) by a number s l i g h t l y greater than one. Figure 12, This flow diagram of the recursive to the multichannel Wiener f i l t e r follows the diagram given by Wiggins and Robinson (1965) as well as a la t e r description given by Robinson (1967a). The polynomials a^Jz) and b.(z) represent z transforms of the prediction and hindsight error f i l t e r s . The sizes of the matrix c o e f f i c i e n t s of these z transforms are determined by the number of input and output channels. In the convention used here for the f i l t e r c o e f f i c i e n t s , f m(k) , k denotes the time index and the subscript m refe r s to the i t e r a t i o n number i n the recursive process. THE WIGGINS - ROBINSON RECURSIVE SOLUTION FOR THE MULTICHANNEL WIENER FILTER INPUTS r (0),r (I), . . . . r(N) g ( 0 ) , g (I),... g (N) INITIALIZATION Q 0 ( 0 ) = b 0 ( 0 ) = I m = 0 om(Z) = ( 0 ) + Q m (1) Z +.. . + Q'm (m) Z m b m (Z) = bm ( 0 ) + b m (1) Z + .. . + b m ( m ) Z m I r ( 0 ) Uam s a m (0) r (m + +... + o m (m) r (I) U bm = U a m T -A K am +1 = -^ bm -r i z ~ U b m ^ m ' 1 • 3. a m + i (Z) = bm+i (Z) = Qm (Z) + K a m + I Z">*' b m (Z ) b m (Z) + K b m + I Z m + ' a m (Z ) a m + , (Z), b m + |(Z) °<m*t = o t m + K o m t , U b m (?m + i : <?m+Kbnrn-iuom f 0 ( 0 ) = r"' ( 0 ) g (0) IN f m (Z) = f m (0) + f m (I) + . . . + f_ (m) Z m %m*-t : fm ( 0 ) r (m+11 + ^ ( 1 ) r (m)+..>fm(m)r(l) K fm +1 = (g (m +i) - 2f m + i) /9 m 11 fm + . (Z)= f m ( Z ) + K f m t l Z < " " b m + l ( Z - ) ' m * i (Z) S E T m E Q U A L T O ' m + i , . ' S T O P WHEN m = N 67 An interesting feature of the multichannel Wiener f i l t e r i s i t s p o s s i b i l i t y for i d e a l performance as pointed out by T r e i t e l (1974). By examining the number of eguations and unknowns for a system of k input channels and a single output channel, T r e i t e l (1974) demonstrates that the multichannel Wiener f i l t e r has an exact solution whenever the f i l t e r reaches a length egual to (N-k)/(k-1)| A , where N i s the length of the input sequence. As seen from the case where k=1, no si m i l a r condition of i d e a l f i l t e r performance exists for the single channel case since the f i l t e r must become i n f i n i t e i n length for an exact solution to be found. In the synthetic example shown in Figure 13(a), the use of Wiener multichannel deconvolution was studied under i d e a l conditions. In order to make r(0) nonsingular, d i f f e r e n t r e a l i z a t i o n s of white noise were added to the wavelet of Figure 4 to produce the wavelets of Figure 13. These were convolved with an impulse response to produce the synthetic traces shown. As required by the preceeding c r i t e r i o n for f i l t e r length, the length of the f i l t e r was chosen to be one less than the length of the source wavelet. When actual source wavelets are used i n the design of the inverse f i l t e r , the convolution of the f i l t e r with the synthetic trace yields a deconvolved trace which almost duplicates the impulse response. As shown i n Figure 13(b), single channel deconvolution followed by trace stacking does not produce id e a l r e s u l t s even when an i d e a l spiking position i s chosen and the actual wavelets are used. Figure 13a). T h i s s y n t h e t i c example e x h i b i t s the performance of the multichannel Wiener deconvolution f i l t e r under i d e a l c o n d i t i o n s . When convolved with the s y n t h e t i c t r a c e s , a multichannel f i l t e r of l e n g t h 14 g i v e s a deconvolved t r a c e which almost d u p l i c a t e s the impulse response. SOURCE WAVELETS I M P U L S E R E S P O N S E Y - A — S r i — 1 1— 0.0 0.2 0.4 TIME —i 1 0.6 C.3 S Y N T H E T I C T R A C E S i 1— 0.0 0.2 I 1— 0.4 0.6 TIME i 1 0.8 1.0 DECONVOLVED T R A C E I 1 1 1 1 1 1 0.0 0.2 0.4 0.6 0.6 1.0 1.2 TIME 70 Figure 13b). A sura of single deconvolutions using a f i l t e r of length 14 designed using the actual wavelets and an optimum spiking position. 71 Although the r e s u l t s of these idealized tests of the multi-channel Wiener f i l t e r are encouraging, the test was conducted under the u n r e a l i s t i c assumption that the source wavelets were known. As shown i n Figure 14(a), the deconvolution deteriorates from the i d e a l when wavelet estimates are used xn designing a multichannel f i l t e r . This example uses two traces from Figures 4 and 13(a) to demonstrate t h i s . The average wavelet estimate, obtained by homomorphic deconvolution, was used to produce the deconvolved traces. Deconvolutions were produced using spiking positions at the f i r s t peak and f i r s t trough of the wavelet with the l e v e l of prewhitening set at 5%. The obvious guesticn which arises at t h i s stage i s whether or not multichannel deconvolution o f f e r s a s i g n i f i c a n t improvement over combinations of single channel deconvolution and stacking. Davies and flercado (1968) pursued t h i s question for re a l data and found that the improvements re s u l t i n g from the use of multichannel deconvolution were i n s i g n i f i c a n t . When the impulse response i s known, a quantitative measure of the success of a deconvolution i s given by the energy normalized crosscorrelation, c ( t ) . If i ( t ) i s the actual impulse response and x(t) i s the deconvolved trace, c (t) i s The sums of single channel deconvolutions f o r the traces of Figure 14(a) are shown i n Figure 14(b) for a f i l t e r of the same K given by (3.18) 72 length and with the same l e v e l of prewhitening as used f o r the multichannel deconvolution of Figure 14(a). The maximum c(t) values for the multichannel deconvolutions were 0.666 and 0.700 for the spiking positions of 3 and 6 respectively. The maximum c(t) values f o r single channel deconvolutions with the same spiking positions were 0.653 and 0.691 respectively. When an optimum spiking position i s chosen f o r the single channel deconvolution of a trace stack, the maximum c (t) value was 0.694. Figure 14a). Multichannel deconvolution res u l t i n g from the use of source wavelet estimation i s shown f o r the cases where the desired outputs are spikes at delays of two units and f i v e units. 74 DELAY=5 Figure 14 b). stacks of single channel deconvolutions designed using the same wavelet estimates, the same spiking positions and f i l t e r lengths and the same prewhitening l e v e l as used i n the multichannel, case of Figure 14a). The thi r d trace shows the resulting deconvolution when the optimum spiking position was chosen, Delay = 2 Delay = 5 .0 Optimum Spike Therefore, as can be seen from Figures 13 and 14 and the values of the energy normalized crosscorrelation, the differences between single and multichannel deconvolution become less s i g n i f i c a n t when wavelet estimates are used instead of the actual wavelets. Let us consider an example which compares the r e s u l t s of single channel and multichannel deconvolution on noisy data. Figure 15 shows 15 seismic traces of near v e r t i c a l incidence r e f l e c t i o n s from 5 CDPs of a Gulf Coast seismic l i n e . Multichannel deconvolutions and stacks of single channel deconvolutions are computed for these bandpassad data, although stacking followed by deconvolution i s the least expensive of the "multichannel" methods, problems can develop when estimating wavelets from stacked data. as pointed out by Dunkin and Levin (1973), normal moveout corrections followed by stacking w i l l d i s t o r t the seismic pulse by enhancing the low freguency content. Figure 16 gives a comparison of the single and multi-channel deconvolution approaches. The deconvolution f i l t e r s have the same length and spiking positions, and are designed on the set of truncated minimum phase wavelets shown i n Figure 16(a). Figure 16(b) displays CDP stacks of single channel deconvolutions, while Figure 16(c) displays multichannel deconvolutions. Both approaches used a 1% l e v e l of prewhitening, the same f i l t e r lengths, and the same spiking position. F i l t e r lengths were sat at the optimum performance length f o r a system of two input channels, which i s the smallest number of channels per CDP i n t h i s set. The normalized mean 78 square error i s always monitored i n order to avoid the numerical d i f f i c u l t i e s which can arise i f the multichannel f i l t e r becomes too long. The r e s u l t s of the deconvolutions in Figure 16 are very s i m i l a r , and the picking of r e f l e c t i o n s on the multichannel deconvolutions i s no easier than on the stacks of sing l e channel deconvolutions. Although the performance of the multichannel approach i s better than the single channel approach for soma synthetic examples, the improvement found by using multichannel Wiener deconvolution on actual data i s often either marginal or nonexistent. This multichannel Wiener f i l t e r e x p l o i t s data redundancy by using wavelet crosscorrelations, whereas single channel methods employ stacking. This accounts for the greater expense required i n using multichannel deconvolution and the higher benefit to cost r a t i o which i s often obtained through use of single channel f i l t e r i n g . 79 flGC ON TRACES SI 7~. T. T. 1 ~i r I.I o i i.o I.J r.t ••• *•* COP 2 --^^\f\J\Af\^ — \ A / ^ -^y\/-^/\f\J\jJ\f\A^^^ CDP 3 -^^v«sf\j\jVv^^ — v A ^ V ^ \ / \ / \ / W \ / W v y \ A A A A ^ — \ / w - x A / COP CDP 5 Figure 15. Data used i n a comparison between si n g l e channel and multichannel deconvolution. Traces show r e f l e c t i o n s of near v e r t i c a l incidence from 5 CDPs. NMO corrections are less than 12 msec. h gain function has been applied to minimize attenuation e f f e c t s . The separation of the CDPs i s 120 feet. Figure 16 a) Truncations of minimum phase wavelets estimated from two egual gates of the traces of Figure 15. Figure 16b) shows CDP stacks of single channel Wiener deconvolutions using f i l t e r s of length 72 msec which were designed to spike the wavelets shown i n a) , while 16c) shows multichannel deconvolutions designed using the same f i l t e r length and spiking position as i n b) . Prewhitening was 1% i n both cases. 81 ESTIMATED WAVELETS i 1 0 . 2 0 . 3 0 . 2 TIME I I 1 0 . 2 0 . 0 _ 0 . 2 TINE b ) STACKED DECON i 1 1 1 1 1 1 1 i O.S 0.1 1 .0 1 .2 1.4 1 .6 l . l 2 . 0 2 . 2 TIME V-Y ^ y Y V V-/j^  C ) MULTICHANNEL "WIENER DECON Y Y w \|yyvvA^^ I 1 1 1 1 1 1 1 , 0 . 6 0 .1 1 .3 1 .2 1.4 1 .6 l . l 2 3 *-2 TIME I 1 1 1 1 1 1 1 1 0.1 0.1 1 .0 1 .2 1.4 l . l l . l 2 . 0 2 . 2 TIME 82 The Kalman deco n v o l u t i o n approach, as d e s c r i b e d by Crump (1974), a l s o a p p l i e s a type of multichannel f i l t e r . The use of the Kalman f i l t e r on the data analysed here produces r e s u l t s which are s i m i l a r to multichannel Wiener d e c o n v o l u t i o n ( Y e d l i n (1975)). The Kalman f i l t e r was a l s o designed on the wavelet estimates shown i n Fig u r e 16, and l i k e any of the methods of i n v e r s e f i l t e r i n g , i t s success depends on the g u a l i t y of data and the accuracy of the wavelet estimates. The time v a r y i n g nature of the wavelet l e a d s us to the next problem i n the design of i n v e r s e f i l t e r s , t h a t o f d e s i g n i n g a time a d a p t i v e f i l t e r . 3.3 Time Adaptive Deconvolutign In d e s i g n i n g f i l t e r s f o r s e i s m i c data, we are i n e v i t a b l y faced with the problem of n o n s t a t i o n a r i t y . To a p p r e c i a t e t h i s concept, we r e c a l l B i c k e r ' s 1953 t r e a t i s e which d e s c r i b e d how the shape of a s e i s m i c pulse changed during propagation. S t r i c t l y speaking, the property of s t a t i o n a r i t y r e q u i r e s t h a t a l l s t a t i s t i c a l moments of a time sequence be time i n v a r i a n t . I f the process examined i s assumed t o be Gaussian, s t a t i o n a r i t y r e q u i r e s that only the f i r s t and second order s t a t i s t i c a l moments ba time i n v a r i a n t . T h i s s t a t i s t i c a l p r o p e r t y , which i s termed weak or wide sense s t a t i o n a r i t y , r e q u i r e s t h a t the mean, v a r i a n c e , and autocovariance be independent of time (Jenkins and Watts (1969)). For a zero mean process such as a seismic trace, tests for. weak s t a t i o n a r i t y require examination of the'autocorrelation function. Hence, the autocorrelation function i s not only necessary i n the design of a Wiener f i l t e r but i s also a means of monitoring nonstationarity. Since s t a t i o n a r i t y i s assumed in the construction of a Wiener f i l t e r , such f i l t e r s are designed on certai n time qates of the data i n order to accommodate th i s assumption. Although the choice cf qate lengths i s often done in a subjective fashion, sophisticated techniques which u t i l i z e autocorrelation estimates may be used to determine a suitable gate length (Wang (1969)). Stationarity tests, however, are not the only c r i t e r i o n for choosing gate lengths. In choosing a gate, we should also r e f r a i n from s p l i t t i n g prominent r e f l e c t i o n events i n the traces. furthermore, although gate lengths must be s u f f i c i e n t l y short to insure some semblance of s t a t i o n a r i t y within the gate, the gates must also be s u f f i c i e n t l y long so that the assumption of a random impulse response i s not badly violated. In choosing gate lengths, i t i s useful to monitor the autocorrelation function. Conventional estimates of the autocorrelation, r (t) , are given by ' 84 N r ( t ) = 7777 ^ x(t-i-z)xc^) (3.19) T s 0 where (x (0) , x (1) ,.,. ,x (N) ) r e p r e s e n t s a time gate of the o r i g i n a l s e i s m i c t r a c e . T h i s v a l u e of r (t) r e p r e s e n t s the a u t o c o r r e l a t i o n o f a sequence which i s z e r o o u t s i d e t h e time g a t e . E x t e n s i o n o f a g i v e n time gate by the a d d i t i o n o f z e r o s i m p l i e s t h a t the a c t u a l a u t o c o r r e l a t i o n o f a time seguence i s m u l t i p l i e d by a B a r t l e t t window. An a l t e r n a t i v e e s t i m a t e of r ( t ) t o t h a t g i v e n by (3.19) r e p l a c e s t h e n o r m a l i z i n g f a c t o r , N+1, by N+1-t. However, t h i s e s t i m a t e can have the u n f o r t u n a t e p r o p e r t y t h a t r ( 0 ) i s not always the l a r g e s t v a l u e of r ( t ) , a p r o p e r t y which i n c o r r e c t l y i m p l i e s t h a t t h e energy of a time seguence can be n e g a t i v e . The i n c o r r e c t assumption o f z e r o e x t e n s i o n f o r the time gate may be remedied by u s i n g t h e maximum e n t r o p y method (MEM) proposed by Burg (1967). The t e c h n i q u e has been d e s c r i b e d i n d e t a i l by Burg (1975) and uMrych and Bishop (1975). The Burg a l g o r i t h m has been o u t l i n e d by Andersen (1974). Burg's method f o r c a l c u l a t i n g t h e a u t o c o r r e l a t i o n uses o n l y t h e a v a i l a b l e d a t a i n the time gate and m i n i m i z e s t h e assum p t i o n s made about the data b e h a v i o u r o u t s i d e the ti m e g a t e . T h i s i s e q u i v a l e n t t o s a y i n g t h a t t h e e s t i m a t e d i s p l a y s maximum e n t r o p y . The MEM c a l c u l a t i o n e s t i m a t e s t h e a u t o c o r r e l a t i o n i n a r e c u r s i v e f a s h i o n by u s i n g the p r e d i c t i o n e r r o r f i l t e r . The Burg a l g o r i t h m can be understood by r e f e r r i n g to the d i s c u s s i o n on p r e d i c t i o n f i l t e r i n g and a u t o r e g r e s s i v e p r o c e s s e s . For Burg's a l g o r i t h m o n l y v a l u e s o f the p r e d i c t i o n e r r o r 85 seguence between t=m and t=N are considered since i n t h i s range a m ( t ) i s defined only in terms of available values of x (t), That i s , (3.20) i s considered for values of t=m,m+1,...„ ,'N. As seen from the above, a m ( t ) i s a portion of the deconvolved trace since Ym{^) i s a minimum phase deconvolution f i l t e r . The hindsight error seguence, b ^ ( t ) , i s determined by convolving x(t) with the hindsight error f i l t e r , where the hindsight error f i l t e r i s defined as the time reverse of the prediction error f i l t e r . Again, only values of b m(t) which have been calculated from available data have been used. That i s , the following expression i s used for t=0,1,2,...,N-m. Considering values of a^ft) and b^ft) i n a r e s t r i c t e d range of t does not require that the data be appended with zeros, Consequently the p r i n c i p l e of maximum entropy i s not violated. 86 The Burg a l g o r i t h m f o r c a l c u l a t i n g % (t) r e q u i r e s the use of the L e v i n s o n a l g o r i t h m to s o l v e the normal e q u a t i o n s w i t h the c o n s t r a i n t t h a t f i l t e r c o e f f i c i e n t s be chosen to m i n i m i z e th e sum of the p r e d i c t i o n and h i n d s i g h t e r r o r powers. The e x t r a c o n s t r a i n t i s n e c e s s a r y s i n c e no a p r i o r i assumptions a r e made c o n c e r n i n g th e a u t o c o r r e l a t i o n e s t i m a t e s {Andersen ( 1 9 7 4 ) ) , I n t h e Burg a l g o r i t h m , the a u t o c o r r e l a t i o n i s c a l c u l a t e d i n a r e c u r s i v e f a s h i o n . As i n t h e L e v i n s o n scheme, the p r e d i c t i o n e r r o r f i l t e r of l e n g t h m+1 i s formed by combining th e p r e d i c t i o n e r r o r and h i n d s i g h t e r r o r f i l t e r s o f l e n g t h ra through use of t h e f i l t e r c o e f f i c i e n t °^m{®) • The updated p r e d i c t i o n f i l t e r c o e f f i c i e n t s a r e g i v e n by t h e r e c u r s i v e r e l a t i o n s h i p o f (3.11). I n u s i n g the above, r e c u r s i o n r e l a t i o n s between a ^ ( t ) and b m ( t ) can be deduced. (3.22) t=m, m + 1,. . , ,N. 87 (3.23) t=0,1 ,. N-m. In the Burg algorithm, the prediction f i l t e r minimizes the sum of the prediction and hindsight error powers. This requires Although the autocorrelation does not have to be e x p l i c i t l y evaluated to determine the Burg prediction error f i l t e r , using the prediction f i l t e r defined by (3.19) and (3.22) as well as the normal equations given i n (3.6) allows the autocorrelation function to be recursively determined. It i s demonstrated by Burg (1975) that r (t) can be calculated by use of the following r e l a t i o n s . r ( o ) = - i _ ^ A l l i * (3.25) a n d rim) = r ^ m - t ) j (3.26) As shown by Ulrych (1972), Burg's method for c a l c u l a t i n g power spectra (or i t s time domain equivalent, the autocorrelation) i s e s p e c i a l l y useful for short data lengths. This makes the Burg method i d e a l for calculating prediction-error f i l t e r s for short time i n t e r v a l s , which i s a necessary property of a time adaptive deconvolution scheme. For longer data lengths, differences between Burg's method and conventional estimates become less apparent. This i s shown for the autocorrelations computed f o r three gates of the trace shown in Figure 17. Moreover, i f only the f i r s t few lobes of the trace autocorrelation are r e l i a b l e estimates of the wavelet autocorrelation because of lack of randomness in the r e f l e c t i v i t y seguence, i t can be argued that the MEM approach does not offer s i g n i f i c a n t improvement over conventional approaches for the gate lengths used in exploration data processing. Figure 18 compares the use of MEM and minimum phase Wiener deconvolution on a noisy synthetic trace. Figure 18(a) displays 89 a s y n t h e t i c s e i s m i c t r a c e formed by adding white no i s e to the c o n v o l u t i o n of a minimum delay wavelet with a spike seguence. F i g u r e 18(b) shows t h a t the Wiener and MEM d e c o n v o l u t i o n s w i l l d i f f e r s i g n i f i c a n t l y f o r s h o r t data l e n g t h s , due to the d i f f e r e n c e s i n the a u t o c o r r e l a t i o n estimates. The l e n g t h of the p r e d i c t i o n e r r o r f i l t e r used i n both p a r t s o f Figure 18(b) i s chosen by use of the f i n a l p r e d i c t i o n e r r o r (FPE) c r i t e r i o n . T h i s c r i t e r i o n , which i s d i s c u s s e d i n Chapter 5, has been d e s c r i b e d by Akaike (1969) and U l r y c h and Bishop (1975). As seen from the values of c (0), the MEM deconvolution g i v e s a b e t t e r estimate of the impulse response than the Wiener d e c o n v o l u t i o n . The c (0) value f o r MEM deconvolution i s 0.348 while the c(0) f o r Wiener d e c o n v o l u t i o n i s 0.268. T h i s r e s u l t i s not unexpected because the MEM a u t o c o r r e l a t i o n estimate i s more accurate. Figure 17. Comparison of autocorrelations using MEM and zero extension methods. Nonstationary synthetic trace i s formed by joining convolutions of wavelets with the impulse responses shown. Autocorrelations are compared for the three gates shown. 91 W A V E L E T S I M P U L S E R E S P O N S E Figure 18. a) Noisy synthetic trace used i n the comparison of Wiener and Burg prediction error f i l t e r i n g . The trace i s formed by convolving a minimum delay wavelet with a spike seguence and then adding 101 noise. a) WAVELET Figure 18 b),c). A comparison between HS.H and Wiener deconvolutions i s shown. In b) the f i r s t 38 samples of the trace are used and i n c) the entire trace i s used. The minimum delay deconvolution f i l t e r s were designed from the autocorrelation estimates shown. F i l t e r lengths were chosen by use of the F.PE c r i t e r i o n , and the deconvolutions were bandpassed to reduce high freguency noise. 96 For c a l c u l a t i n g a u t o c o r r e l a t i o n s f o r s h o r t data l e n g t h s , the MEM approach i s e s s e n t i a l . Such an approach was a p p l i e d by R i l e y and Burg (1972) i n the d e s i g n of a time a d a p t i v e d e c o n v o l u t i o n scheme which t o seme e x t e n t o b v i a t e s t h e use of t i me g a t e s i n f i l t e r d e s i g n . T h i s approach c o n t i n u o u s l y updates t h e p r e d i c t i o n e r r o r f i l t e r t hroughout th e t r a c e . The p r e d i c -t i o n e r r o r f i l t e r i s c a l c u l a t e d from th e d a t a w i t h i n a s l i d i n g window whose l e n g t h i s s e t e g u a l t o the f i l t e r l e n g t h . The f i l t e r i s deduced from th e p r e d i c t i o n e r r o r and h i n d s i g h t e r r o r seguences c a l c u l a t e d from w i t h i n t h e window. The R i l e y - Burg approach a p p l i e s e x p o n e n t i a l w e i g h t i n g t o p r e d i c t i o n e r r o r and h i n d s i g h t e r r o r seguences i n the c a l c u l a t i o n of t h e v a l u e s so t h a t the p r e d i c t i o n e r r o r s f o r p o i n t s c l o s e s t t o t h e p r e d i c t e d v a l u e a r e weighted most h e a v i l y . The v a l u e of the p r e d i c t i o n e r r o r seguence, which has been c a l c u l a t e d i n the c o n t i n u o u s u p d a t i n g o f t h e p r e d i c t i o n e r r o r f i l t e r , i s e s s e n t i a l l y t h e d e c o n v o l v e d t r a c e . F i g u r e 19 shows an a p p l i c a t i o n of t h e R i l e y - B u r g approach which uses an a d a p t i v e f i l t e r c a l c u l a t e d o v e r a s h o r t t i m e gate of 40 msec. The r e s u l t a n t d e c o n v o l u t i o n has a very s p i k e d appearance but the s p i k e s are sometimes s h i f t e d and c f the wrong s i g n . For t h i s r e a s o n , the e n velope of the d e c o n v o l u t i o n s i s a l s o shown and compared t o t h e e n v e l o p e of t h e i m p u l s e response. Problems can a r i s e w i t h t h e R i l e y - B u r g approach when the gate l e n g t h becomes so s m a l l t h a t t h e assumption of a random impulse r e s p o n s e b r e a k s down and the t r a c e a u t o c o r r e l a t i o n does not g i v e a good e s t i m a t e of the wavelet a u t o c o r r e l a t i o n . I n t h a t c a s e , t h e p r e d i c t i o n -97 error f i l t e r i s not suitable for removing the wavelet from the trace. Also, to obtain a deconvolution that w i l l resemble the impulse response shown in Figure 19 requires that bandpass f i l t e r i n g be applied to the prediction error seguence. Although problems exist in t r y i n g to remove seismic wavelets by the Riley-Burg approach, the method i s appropriate f o r "deglitchinq" teleseismic data, (Davies (1976) ), since "glitches" are large i s o l a t e d prediction errors in the data which can be picked out well by short MEM prediction error f i l t e r s . Figure 19. An application of Riley-Burg time adaptive deconvolution. A continuously time adapting f i l t e r of length 20 msec (or 10 samples) was applied to a bandpassed input trace formed by the convolution of a wavelet with an impulse response deduced from a well log synthetic. The envelope of the deconvolution, or prediction error seguence, and the envelope of the impulse response are compared. SOURCE WAVELET I M P U L S E R E S P O N S E i 1 »•» I N P U T T R A C E BURG DECON it i . i . J . A Jill j! , J i J 1 .ll A i ll- I.i L.l 1 L i, li. 11 i. Ii. .L • , 1— p^'^j^N", i M .0 0.2 0.4 0.6 "] 1' 0 1 ' T TIME 2.2 IMPULSE RESPONSE ENVELOPE i r 0.0 0.2 T 0.4 i r 0.6 0.8 T — r 1.0 1.2 TIME 1.4 1.6 1.8 2.0 2.2 DECON E N V E L O P E 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 l.B 2.0 2.2 TIME VO VO 100 Gutowski (personal communication) has suggested an HEM approach which provides a good compromise between the time adaptive deconvolution suggested by Riley and Burg (1972) and the conventional deconvolution of larger data gates. This modified MEM approach also uses a s l i d i n g rectangular time gate i n which the prediction error f i l t e r i s constantly updated. Unlike the Riley-Burg approach, the f i l t e r i s not the same size as the window, but i s usually smaller. This modified time adaptive MEM approach was applied to a nonstationary trace model shown i n Figure 20, The CVL synthetic from 0.4 to 2.2 sec was used, as well as a wavelet which changed shape every 600 msec. Noise was added to the convolution of the time varying wavelet with the impulse response . A window size of 120 msec was used i n the design of a 50 msec prediction error operator. Bandpass f i l t e r i n g was applied to reduce high frequency noise i n the deconvolution. It i s often useful to calculate the envelope of time adaptive MEM deconvolutions by using the Hilbert transform ( Bracewell (1965)). Clayton et a l (1975) have demonstrated that the envelope can be very useful i n determining teleseismic a r r i v a l times. Figure 20 compares the envelope of the deconvolution with the envelope of the impulse response to demonstrate the usefulness of t h i s c a l c u l a t i o n for the determination of r e f l e c t e d seismic a r r i v a l s . Although the envelope c a l c u l a t i o n s a c r i f i c e s phase information, i t provides complementary information to the time adaptive deconvolutions of r e f l e c t i o n data. Figure 20. An application of time adaptive MEM deconvolution. wavelets of d i f f e r e n t lengths were convolved with sections of an impulse response deduced from a CVL. The position of the wavelets indicates whic portion of the impulse response they were convolved with A time adaptive f i l t e r of length 50 msec was applied to the trace and bandpassed to produce the deconvolution. The envelopes of the deconvolution and the bandpassed impulse response are also compared. 102 S O U R C E W A V E L E T S I M P U L S E R E S P O N S E I N P U T T R A C E i r 0.0 0.2 0.4 D E C O N V O L U T I O N J llhkM II lull JUJ1 ilk i l l f Mil/ (My * vv| I! |r|l IP i f flip Iff • 1 r 0.0 0.2 3.4 f 1 E N V E L O P E B A N D P f l S S E D I M P U L S E R E S P O N S E 2.0 2.2 D E C O N E N V E L O P E M l l r— • 1 1 " — l —i 1 1 : 1 r 0.0 0.2 0.4 0.6 0.0 1.0 1.2 1.4 1.6 l .B 2 0 2 2 TIME 103 CHAPTEE 4 A CASE HISTORY QF DECONVOLVING REAL SEISHIC DATA H'l Deconvolutions Performed Several deconvolution techniques were applied to seismic data from the Gulf Coast area. The deconvolution r e s u l t s were compared to synthetics obtained from well loqs i n order to int e r p r e t qeologic boundaries. The data set analysed i s from f i v e common depth points of dynamite generated land data. Parts of t h i s data set were shown e a r l i e r in Figures 6 and 15. In order to follow the assumptions used in the trace model of Chapter 1, most of the processing was done on traces of seismic waves reflected at near v e r t i c a l incidence. For these traces normal moveout corrections are small and the s t a t i c s corrections are 4 msec or l e s s . In most exploration work, large amounts of data are used and the assumption of v e r t i c a l incidence i s s a c r i f i c e d so that the noise l e v e l i n the traces can be decreased by stacking. The purpose of t h i s study was to compare the effects of several deconvolution techniques on a limi t e d amount of data rather than to process a large amount of data with any given method. Due to the noise in the s i g n a l , the data were bandpass f i l t e r e d using an Ormsby f i l t e r which passed freguencies between 10 and 55 hz. As shown i n Figure 21, the attenuation losses i n the data also had to be compensated for. This was accomplished by finding the least sguares f i t cf an exponential curve to the 104 a b s o l u t e a m p l i t u d e of the t r a c e and then m u l t i p l y i n g t h e t r a c e by the r e c i p r o c a l o f the e x p o n e n t i a l c u r v e . T h i s i s a form of g a i n c o n t r o l or s p h e r i c a l d i v e r g e n c e c o r r e c t i o n . I t s e f f e c t can be seen by comparing the t r a c e s o f F i g u r e 2 1 . Some of t h e d e c o n v o l u t i o n s f o r t h e s e t o f f i v e common depth p o i n t s were shown e a r l i e r . D e c o n v o l u t i o n s o f thes e d a t a were performed i n the comparison o f wavelet e s t i m a t o r s and i n the comparison o f s i n g l e c h a n n e l and m u l t i c h a n n e l d e c o n v o l u t i o n s . A comparison of Wiener d e c o n v o l u t i o n s on u n s t a c k e d d a t a i s g i v e n i n F i g u r e 22 . F i g u r e 22(a) shows a Wiener minimum phase d e c o n v o l u t i o n w i t h no p r e w h i t e n i n g w h i l e F i g u r e 22(b) g i v e s a Wiener minimum phase d e c o n v o l u t i o n w i t h 10% p r e w h i t e n i n g . These f i g u r e s demonstrate t h a t p r e w h i t e n i n g a l l o w s f o r a s t a b l e d e c o n v o l u t i o n w i t h l e s s h i g h f r e g u e n c y n o i s e w h i l e s a c r i f i c i n g r e s o l u t i o n . A l s o , s i n c e few r e f l e c t i o n s can be e a s i l y i d e n t i f i e d on t h e s e c t i o n i n F i g u r e 22 , s t a c k i n g of t h e data i s a d v i s a b l e . B e f o r e p r o c e e d i n g t o compare the d e c o n v o l u t i o n s w i t h s y n t h e t i c s , one more example of the e f f e c t s o f p r e w h i t e n i n g i s shown. F i g u r e 23 d i s p l a y s d e c o n v o l u t i o n s f o r the CDP d a t a s t a c k s of d a t a from F i g u r e 15 . As can be seen from t h e f i r s t t r a c e o f t h e d e c o n v o l u t i o n s i n F i g u r e 2 3 ( b ) , edge e f f e c t s can o c c u r when u s i n g no p r e w h i t e n i n g . On t h e o t h e r hand, t h e use o f 101 p r e w h i t e n i n g , as shown i n F i g u r e 2 3 ( c ) , e l i m i n a t e s t h i s problem. L e t us now examine s y n t h e t i c seismograms wtiich were o b t a i n e d from v e l o c i t y l o g s near t h e s e i s m i c l i n e . Figure 21a). Bandpass f i l t e r e d data from an area near the Gulf Coast. Distance between CDPs i s 120 f e e t . Figure 21b) shows the data after spherical divergence corrections were applied. Such corrections compensate for the energy losses in spherical wavefront amplitudes, which are proportional to the r e c i p r o c a l of the distance squared. These data were used in subsequent deconvolution t e s t s . 106 a) B A N D P A S S F I L T E R E D T R A C E S CDP 1 COP 2 — COP ^™*J\f^S\l^^^ 3 ^\/\A^AAA/^ V V v ^ v w v ~ n ^ v ' - w - ^ ~ ^ . COP 4 " n - - - = J . b > A G C O N T R A C E S i i 1 1 1 —r °/? lJf COP 1 I I i 1 O.t 0 . 9 t.O 1.2 1.4 1.6 1.0 2 . 0 1.1 X ^^^/\f\j\^ A / V -X /-^AJ\I\!\/YA!\^ • ^  -wvyyvwW\/\/\AMA/v\^ •^^\J\J\y-\j\ApJ^^ ^ - N / — / v -— • — • * s ^ J \ J \ J \ f \ J \ [ \ / ^ ^ — v / w w \ , COP 5 : s^\/W^r\f^^ ^ V V v v \ / \ / \ / W ^ 0-6 0.8 1.8 1.2 I 1 I ! 1 TI& 1,6 '•' " Figure 22 a). Wiener minimum phase deconvolution of the data i n Figure 21a) using no prewhitening and a f i l t e r length of 80 msec. Figure 22b) shows a minimum phase deconvolution of the data i n Figure 21a) using a f i l t e r of the same length with 10% prewhitening. 108 a> W I E N E R M I N P H A S E D E C O N I 1 1 1 : 1 1 1 1 0.8 0.8 1.0 1.2 1.4 1.8 l . l 2 .0 2.2 i 1 1 1 1 1 1 1 1 O.I 0 .9 1.0 1.2 1.4 1.6 1.0 2.0 2 .2 TIME to W I E N E R M I N P H A S E D E C O N i 1 1 1 1 1 1 1 1 0.0 0.8 1.0 1.2 1.4 1.6 1.6 - 2 .0 2.2 i 1 1 1 1 1 1 1 1 0.1 0 .8 1.0 1.2 1.4 1.6 1.9 2.0 2.2 TIME Figure 23 a). CDP stacks cf the data in Figure 21b). Wiener deconvolutions of stacks with no prewhitening are shown i n b), while similar deconvolutions of stacks with 10% prewhitening are shown in c ) . a , C D P STACKS i — 0.6 ~1 1.0 1.2 1 1 .* TIME - 1 — 1.6 i.e i'.o 2.? fo WIENER DECON ON SOURCE EST c) WIENER DECON ON SOURCE EST A A A J 111 4.2 Impulse Response Models Obtained From Synthetic S e i s a o j r a m s Estimates of the earth's impulse response i n the region of i n t e r e s t sere obtained by using two sets of seismic velocity information. These included a continuous v e l o c i t y log (CVL) and a velocity-depth interpretation obtained from check shot data. Synthetic seismograms of the earth's response to a spike impulse were obtained by using the approach of Berryman et a l (1960). a general program which applied t h i s approach was presented by Clowes (1969). A version of this program, revised by C. Chapman to include a source at depth and r e f l e c t i o n s from the free surface, was used to obtain synthetics f o r both sets of velocity information. Although density data ware also a v a i l a b l e , the method of Berryman et a l (1960) uses information pertaining to velocity variations i n the layers only. Also, the density logs from t h i s area are often unreliable. Furthermore, i n t h i s case, synthetics which used velocity logs were very similar to those obtained by using velocity and density (Lucas (personal communication)). This s i m i l a r i t y occurs because the percentage change i n velocity i s usually larger than the percentage change i n density so that the value of the r e f l e c t i o n c o e f f i c i e n t i s largely determined by velocity variations. Also, density variations i n t h i s case often coincided with large v e l o c i t y changes. The largest density variation for the well logs occured at a depth of 2200 feet where the density decreased from 2.08 gi/cc to 1.48 gm/cc in a distance of 20 feet. This large density change was not c l e a r l y shown in the seismic sections but 112 did coincide with a velocity change found in the check shot v e l o c i t y data. I t should be noted that using only v e l o c i t y variations to construct synthetics i s not va l i d i n the presence of gas deposits because the large density variation for these "bright spots" gives r i s e to large amplitude r e f l e c t i o n s . Also, i t should be pointed out that the synthetics were constructed such that attenuation was n e g l i g i b l e . For t h i s , a value of Q equal to 1000 was used. In addition, spherical divergence corrections were applied to the seismic traces i n Figure 21. The two models of i n t e r v a l v e l o c i t y versus depth are shown i n Figure 24 along with synthetics deduced from these models. The f i r s t velocity curve was obtained by recording surface shots using a borehole geophone placed at various depths. These ve l o c i t y data are termed check shot data since the velocity curve obtained by t h i s technique can be used to check the r e s u l t s of integrating a continuous velocity log. Figure 24. CVL and interpreted check shot velocity models with synthetic seismograms derived from these v e l o c i t y curves. The f i r s t synthetic of the check shot v e l o c i t i e s places the shot at the surface and contains n ghosts. While the second check shot synthetic contains multiples for the case where the shot i s buried at a depth of 100 feet. The CVL synthetic contains multiples for a shot at the same depth. INTERPRETED VELOCITY LU 0.0 1.0 _ J 2.0 _ j DEPTH IN KFT 3.0 4.0 5.0 _ l 1 L _ '.-CVL VELOCITY UJ 1.0 2.0 _ J 3.0 _ J DEPTH IN KFT 4.0 5.0 l i INTER VEL SYNTHETIC PRIMARIES ONLY o.o o.o o.o 115 The second velocity-depth model was obtained by using a continuous velocity log which was sampled at one foot i n t e r v a l s . In order to obtain a model which was computationally f e a s i b l e and geologically r e a l i s t i c , the velocity values of the CVL were averaged over 25 foot i n t e r v a l s . Since the logging sonde supplies r e l i a b l e velocity values only below the well casing, the i n i t i a l part of the CVL model was found by using the v e l o c i t y values from check shot data. Using the interpreted velocity from check shot data, synthetics were calculated to give only primary r e f l e c t i o n s in one case and primaries plus multiples i n the second case. In the f i r s t case the shot was placed at the surface and i n the second case the shot was buried at 100 feet. The CVL curve was used to calculate a synthetic which included primaries and multiples. The l a t t e r synthetic gives an i n d i c a t i o n of the geologic complexity in t h i s region. The delta function source used i n the synthetic was determined by using a f l a t spectrum between 0 and 250 hz. The synthetic obtained using a CVL supplies a model for the impulse response of a layered earth i n the v i c i n i t y of the seismic l i n e . In the minimum phase deconvolution approaches examined previously, the impulse response was assumed to be a random or white noise sequence. Therefore, i t s autocorrelation should be a delta function at zero delay. The autocorrelation shown in Figure 25 gives some support for t h i s assumption since the spike at t=0 dominates the autocorrelation. However, the e f f e c t of a ghost r e f l e c t o r i s also shown in the 116 autocorrelation. A ghost r e f l e c t o r i s also evident i n the second synthetic i n Figure 24 which includes the multiple r e f l e c t i o n s for the check shot velocity model. A ghost r e f l e c t i o n r e s u l t s when energy from the shot propagates to the surface and i s r e f l e c t e d down again. Ghost r e f l e c t i o n s have the opposite sign from the primary r e f l e c t i o n s and are separated in time frem the primary r e f l e c t i o n s by an amount equal to twice the shot to surface t r a v e l time. In t h i s case the shot depth was set at 100 feet i n the synthetic. As w i l l be seen l a t e r , the qhost r e f l e c t i o n i s not evident i n the autocorrelation of the seismic trace or the autocorrelations of subsequent deconvolutions. The ghost r e f l e c t o r i s present i n the synthetics because of the oversimplified velocity model which had to be used for the f i r s t 1630 feet of sediment. This should be kept i n mind when correlating the synthetics with seismic data or deconvolutions. Since the data ware bandpass f i l t e r e d to suppress low frequency surface, wave information and high frequency noise, the best that can be expected from conventional deconvolution i s a bandpassed version of the impulse response. The term "conventional" i s used here since Chapter 5 considers an unconventional method of extending the spectrum of the impulse response. The bandpassed synthetic and i t s autocorrelation are also shown i n Figure 25. As expected, bandpass f i l t e r i n g has the undesirable e f f e c t of decreasing the randomness of the impulse response while having the desirable effect of decreasing the noise contained i n the seismic trace. 117 Figure 25, Comparison of the autocorrelations of the CVL synthetic before and after bandpass f i l t e r i n g . Bandpass f i l t e r i n g was accomplished by using the same Ormsby f i l t e r as used on the seismic data. The corners of t h i s f i l t e r i n the frequency domain were at 10, 15, 45, and 55 hz. Note the.effect of the "ghost " r e f l e c t i o n on the autocorrelation. A U T O C O R R E L A T I O N 0.2* ' 0.4 0.6 0.8 1.0 1.2 1.4 1.6 . . . . TIME: n : — i 1 1.6 2.0 2.2 "ghost effect" AUTOCORRELATION AFTER BANDPASS 0.2V D.4 0.6 o . i ' ° T I M E 1 , 4 J " 6 n—' 1.8 ~1 2.2 1.4 1.  2.0 BANDPASS OF SYNTHETIC. i — 0.0 0.2 7$ i.2 Co 119 j 4 . 3 Relating Be convolutions and Synthetic Sei.siograms An int e r p r e t a t i o n of the tra v e l times for seismic r e f l e c t i o n s from various horizons i s obtained by using the synthetic seismograms and deconvolutions. F i r s t of a l l , the synthetic obtained from the CVL i s correlated with the most coherent set of deconvolved traces. Secondly, the best r e f l e c t o r s from a l l the deconvolutions and synthetics are picked and compared. The most coherent set of seismic traces was obtained by stacking traces with the same shot to receiver o f f s e t distance. These common offset trace gathers are mentioned by Wood and T r e i t e l (1S75). The stacking of traces with a common offset distance i s v a l i d over short distances in which the layers can be considered to be horizontal. The of f s e t stacks were deconvolved by using Wiener f i l t e r s which were designed on wavelets obtained by the Hilbert transform method. Successive deconvolutions were performed on the data i n order to increase the s i m i l a r i t y between the deconvolution and the impulse response. With each deconvolution, a new wavelet i s estimated and removed. This approach, which i s s i m i l a r to the i t e r a t i v e deconvolution described by Robinson (1975), represents a method of compressing seismic pulses in an i t e r a t i v e fashion. Figures 26(a) and 26(b) display the f i r s t and second deconvolutions of the offset stacks. The figures also show that the wavelet estimate i s compressed with each deconvolution. As shown by the autocorrelations of Figure 27, the randomness and high frequency content of the traces increases with each step in the i t e r a t i v e deconvolution process. 120 Figure 26. I t e r a t i v e deconvolution of the o f f s e t stacks for near v e r t i c a l incidence. The wavelet size i s compressed with each i t e r a t i v e step, and the high frequency content of the traces i s increased. 121 i 1 1 1 1 1 1 1 1 0.5 0.8 1.0 1.1 1.4 1.6 1.8 2.0 2.2 OFFSET STACKS Tlie O ft . « • . ! • WIENER DECONVOLUTIONS WIENER DECONVOLUTIONS BANDPRSS OF IMPULSE RESPONSE ^UTQCQRR WAVELETS V TIME 0 1 rir° 1 0,V 0.2 0.0 0.2 * TIME O^ tt 0.2 TIME ttJ QCORR WAVELETS rtr 1 •! Jo 0.2 i . « 0.2 TIME — i ' r 1 .4 I. TIME — i 1 2.0 2.2 I— 0.6 o'.S 1.0 - 1 — 1.2 6 1.8 Figure 27. Comparison of autocorrelations at various stages in the i t e r a t i v e deconvolution process. I t can be seen that the deconvolutions increase the high frequency content and randomness. No "ghost" r e f l e c t o r i s evident i n the autocorrelations of the deconvolutions. TRACE AUTOCORRELATIONS 3 tty 0.2 0.4 O.E 0.6 1.0 1.2 1.4 f TIME 1 W v YT O rt A rt r- ~ -l | | V W W! - 1 ' W V ^ ^ U <-» - | ~ „ ~ . r - . "> r * - - ^ ^ r i ^ °VH a' 2 °-4 0-6 O i 1.0 1.2 1.4 ^ T iME 0.2 0.4 0.6 0.8 l.o 1 2 TIME 1.4 JECQN AUTOCORRELATIONS °V* ° - 2 O'-B 0.6 l.O 1.2 1.4 TIME 0.2 0.4 0.6 0.8 1.0 1.2 TIME ^ • _ . ' ~ » L T . ' - v k y V , r » , — p 1.4 >-P 0.2 0.4 0.6 0.8 l.o 1.2 1.4 TIME ^ECOND DECON AUTOCORRELATIONS ? - 1 1 r -TI°M8E r |7yV \f t \ /vH^AV vA- • " i^/ v^^A^^v / V l " ^ 1 w i""-"' — ^ r - • T • 1 rH> 0.2 0.4 0.6 O B 1.0 r . 2 TIME 0.2 0.4 0.6 0.6 1.0 1 2 TIME I M P U L S E RESPONSE AUTOCORRELATION ~i— 1.4 124 As can ba seen from Figure 26, the agreement between the s e i s m i c data and tha s y n t h e t i c s i s not good. In f a c t , the energy normalized c r o s s c o r r e l a t i o n suggests that the s y n t h e t i c be s h i f t e d t o the l e f t by 20 msec to agree with the s e i s m i c r e s u l t s . Even so, the maximum value of the energy normalized c r o s s c o r r e l a t i o n f o r a s h i f t of 20 msec was only 0.15. Although the c r o s s c o r r e l a t i o n s between the s y n t h e t i c s and the deconvolutions were poor i n t h i s case, i t was encouraging to note t h a t s i m i l a r s e i s m i c d e c o n v c l u t i o n r e s u l t s f o r these data were independently obtained by Amoco Beseach. The Amoco r e s u l t s which are shown i n F i g u r e 28 d i s p l a y the prominent r e f l e c t i o n s at almost the same a r r i v a l times as those i n F i g u r e 26. The Amoco group experienced the same d i f f i c u l t y i n c o r r e l a t i n g s y n t h e t i c s with s e i s m i c d e c o n v o l u t i o n s . 125 Figure 28. Deconvolution r e s u l t s obtained by Amoco Research for the data set analysed i n t h i s chapter. 126 There are several reasons why impulse response synthetics and deconvolutions may not be compatible. F i r s t of a l l , the well's position i s not at the same position as the seismic l i n e . In the case examined here, the seismic data and the well log information were gathered at a distance of about 800 feet apart. Secondly the nature of a CVL i s such that no velocity information i s available for sediments above the bottom of the well casing. Hence, r e f l e c t i o n s and multiples r e s u l t i n g from near surface v e l o c i t y d i s c o n t i n u i t i e s w i l l not be included in the synthetic, although such a r r i v a l s can be very prominent on an actual trace. Also, well log measurements are taken i n the geologically disturbed environment of a borehole, and the problem of mud seepage can cause the velocity information to be inaccurate. Furthermore, the actual seismic trace w i l l not consist s t r i c t l y of a P wave r e f l e c t e d at v e r t i c a l incidence by horizontal layers. The trace w i l l include surface wave energy, S wave a r r i v a l s , as well as refracted, d i f f r a c t e d and converted a r r i v a l s . At this point, an interpretation i s obtained by merging the deconvolution and synthetic seisiogram information. Figure 29 summarizes the most probable r e f l e c t o r s which have been picked from the previous figures. Figure 29. Summary of the best r e f l e c t i o n picks from the deconvolutions and synthetic seismograms presented. The di f f e r e n t arrows denote r e l i a b i l i t y of the picks. R E F L E C T I O N P I C K S F R O M D E C O N V O L U T I O N S A N D S Y N T H E T I C S E I S M O G R A M S S O U R C E OF INFORMATION D E C O N O F FIG. IO ( o f r e r N M O ) D E C O N O F FIG. 16 D E C O N O F FIG. 2 2 D E C O N O F FIG 2 3 D E C O N O F FIG. 2 6 D E C O N O F FIG. 2 8 C H E C K S H O T S Y N T H E T I C (FIG. 2<*) C V L S Y N T H E T I C (FIG. 2 6 ) e . o 8 1. O 1. 2 1. <* IA 5 1. 3 2 O 2 2 1 i 1 t 1 t • 1 • t i J j • 1 i t I" I 4 i I t t t 1 t t t t t 1 1 1 1 t t t t * 1 1 11 : • 1 1 • 1 l • 1 1 t t ! i T I M E t .1 ( 3 0 5 0 Ft.) ( 4 ISO Ft.) t t ( 5 3 7 0 Ft.) ( 6 0 0 0 Ff . ) ( 8 3 0 0 Ft.) Q U E S T I O N A B L E I ho 00 129 From the well log synthetics and the seismic deconvolutions, cer t a i n r e f l e c t o r s are consistently shown in the i n t e r v a l between 0.600 and 2.200 seconds. These appear at times of 0.92, 1.20, 1.48, 1.63, and 2.12 seconds. For the time gate used, the depth range of r e f l e c t o r s was approximately between 1900 and 8700 feet. Osing the average v e l o c i t i e s deduced from check shot data, allows us to deterine the depth of the r e f l e c t i n g interfaces. These interfaces were situated at depths of approximately 3050, 4190, 5370, 6000 and 8300 f e e t . Both wall log synthetics also indicate an interface at a depth of 6700 feet but this i s not evident i n the seismic data. Although a plausible interpretation can be obtained from the information available, i t i s clear that no single deconvolution method w i l l provide miraculous results on such noisy data. Having considered deconvolutions on dynamite generated data, attention i s now focused on data generated by vibrator sources. 130 CH1FTJB 5 A JEW APPROACH TO VIBROSEIS DECONVOLUTION 5*1 I s i i o d u c t i o n t o t h e V i b r o s e i s Method The f o l l o w i n g d i s c u s s i o n , t a k e n from L i n e s and C l a y t o n (1976), o u t l i n e s a new approach f o r i m p r o v i n g t h e r e s o l u t i o n o f v i b r o s e i s r e c o r d s . C e r t a i n p r o p e r t i e s of the v i b r o s e i s s i g n a l make i t p a r t i c u l a r l y s u i t a b l e f o r u s i n g a d e c o n v o l u t i o n t e c h n i q u e which combines homomorphic d e c o n v o l u t i o n and a u t o r e g r e s s i v e p r e d i c t i o n . The t e c h n i q u e has a l s o r e c e n t l y been a p p l i e d t o o t h e r d e c o n v o l u t i o n problems by C l a y t o n and 'Wiggins (1976) and C l a y t o n and U l r y c h (1976) . In t h e v i b r o s e i s approach, a f r e g u e n c y modulated s i q n a l i s ge n e r a t e d u s i n g a m e c h a n i c a l v i b r a t o r . U n l i k e s i q n a l s produced by dynamite e x p l o s i o n s , the f r e q u e n c y range of t h e v i b r o s e i s s i g n a l can be c o n t r o l l e d . Since the a m p l i t u d e of a v i b r o s e i s sweep i s s u b s t a n t i a l l y l e s s than the a m p l i t u d e o f a dynamite g e n e r a t e d p u l s e , t h e time d u r a t i o n o f the s i q n a l i s i n c r e a s e d i n o r d e r t h a t t h e enerqy a s s o c i a t e d w i t h t he v i b r o s e i s t r a c e be c o m p a t i b l e w i t h t he energy o f t r a c e s g e n e r a t e d by dynamite e x p l o s i o n s . T h i s concept was implemented by Klauder e t a l . (1960) i n the r e s o l u t i o n o f c h i r p r a d a r s i g n a l s . A n a l y s i s o f the v i b r o s e i s s i g n a l r e q u i r e s t h a t t h e i n p u t f r e q u e n c y modulated s i g n a l be c r o s s c o r r e l a t e d w i t h t h e r e c o r d e d s i g n a l i n o r d e r t o "compress" t h e energy i n the v i b r o s e i s t r a c e so t h a t i t s appearance i s s i m i l a r t o a t r a c e g e n e r a t e d by a h i g h energy 131 dynamite s o u r c e . 5.2 A V i b r o s e i s Trace S o d e l A n o i s e l e s s u n c o r r e l a t e d v i b r o s e i s t r a c e , y ( i ) , can be c o n s i d e r e d as the c o n v o l u t i o n of the f r e q u e n c y modulated (FM) i n p u t sweep p ( t ) , w i t h the i m p u l s e response of a l a y e r e d e a r t h i ( t ) . y(t) = p(t) * c(t) (5.D The symbol * i s used throughout to denote c o n v o l u t i o n . The f r e q u e n c y modulated sweep p ( t ) , can be r e p r e s e n t e d as p ( t ) = A sin 2n-F0 ( t t | t * ) ( 5 - 2 ) The i n s t a n t a n e o u s f r e q u e n c y i s q i v e n by f 0 ( 1 + k t ) , which i s JTT -J— , where KJ xs t h e phase. I f T xs the d u r a t x o n of the sweep, t h e f r e q u e n c y of the i n p u t FM sweep ran g e s between f Q and f 0 (1+kT) . C r o s s c o r r e l a t i o n o f y (t) w i t h t h e i n p u t FM sweep y i e l d s an i d e a l i z e d model of a v i b r o s e i s t r a c e , x ( t ) . The f u n c t i o n x ( t ) may be r e p r e s e n t e d i n t h e f o l l o w i n g c o n v o l u t i o n a l f orm. x(t) = p(-t) * p U ) xid) 132 (5.3) or X(i) = Wit) $ (Ct) (5.4) In (5.4), w (t) i s the autocorrelation of the input FM sweep and represents a zero phase sequence termed the Klauder wavelet, A model of such an idealized trace i s shown i n Figure 30. This figure shows a vibrator sweep, p (t) , i n which f0=14 hz, T=1.0 sec, and k=3. Although the sweep i s generally much longer than one second, t h i s example i l l u s t r a t e s the properties of the trace model. Figure 30 shows the Klauder wavelet or autocorrelation of p (t ) , and an impulse response of random spikes. The recorded vibroseis trace, y (t) , and the crosscorrelated trace, x (t ) , are also given. The power spectrum of x(t) mainly shows the influence of the impulse response, i ( t ) . Figure 30. Example of an idea l i z e d model for a vibroseis trace. The figure shows the input FM sweep, the sweep's power spectrum, the Klauder wavelet or autocorrelation of the sweep, an impulse response model, the recorded vibroseis trace, the crosscorrelation of the sweep with the recorded trace, and the power spectrum of the crosscorrelated trace. 134 NPUT FM j.O POWER SPECTRUM AUTOCORRELATION -1 1 — 3.0 62.5 123.0 FREQUENCY i 1 1 r 1 1 1 0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 TIME IMPULSE RESPONSE 1 1 TIME 1.0 VIBROSEIS TRACE A 0.0 0.2 CROSSCORRELATED,TRACE POWER SPECTRUM J.D 0.2 TIM o.6 V V b.a i.o r*- " i 1 0.0 62.5 125.0 FREQUENCY 135 As shown by the figure, the power spectrum of the input sweep i s rather f l a t within the sweep's freguency l i m i t s (14-56 hz). The power spectrum of the input sweep i s the amplitude spectrum of the Klauder wavelet. Therefore, deconvolution of t h i s i d e a l i z e d v ibroseis trace within the freguency l i m i t s of the sweep i s f u t i l e , since a wavelet with a f l a t spectrum w i l l produce a trace which has the same amplitude spectrum as that of the bandpassed impulse response. Since any seismic signal transmitted through the earth i s dispersed and attenuated, a r e a l i s t i c model of y (t) considers the input FM signal to be f i l t e r e d by a function e{t) which describes the attenuation and dispersion properties of the earth. The fact that the input vibrator s i g n a l i s also not i d e a l may be incorporated into a (t), E s s e n t i a l l y e(t) represents the f i l t e r i n g effects of the earth on p(t) and w i l l tend to decrease the high freguency content i n y (t). In our model, the s i g n a l returning from a particular r e f l e c t o r w i l l be represented by p (t)*e (t) instead of p (t) and the new model for the trace i s given by (5.5) Hence, the crosscorrelated vibroseis trace, x ( t ) , i s then given by x L t ) = w ( t ) * ec-t) * i(-t) 136 (5.6) There i s evidence to suggest t h a t e (t) i s minimum phase. In t h i s r e s p e c t , the previous d i s c u s s i o n s of Futterman's minimum phase a t t e n u a t i o n laws are of i n t e r e s t . 5.3 The V i b r o s e i s Jjeconvglutign Probl-em Since the impulse response or r e f l e c t i v i t y seguence, i (t) , i s of g e o l o g i c a l i n t e r e s t , the seguence w(t)*e(t) i s the wavelet that we wish t o remove. As pointed out by Sistow and Jurczyk (1975), t h i s wavelet w i l l , i n g e n e r a l , be n e i t h e r minimum phase nor zero phase. T h i s makes the deconvolution d i f f i c u l t . T h e o r e t i c a l l y , the spectrum of the impulse response i s continuous up to the Nyguist freguency, but only a c e r t a i n bandwidth i s e x c i t e d because of the b a n d l i m i t e d source and the a t t e n u a t i o n and d i s p e r s i o n e f f e c t s of the earth. These e f f e c t s cause the v i b r o s e i s wavelet to have d i s t u r b i n g s i d e l o b e s , thereby l i m i t i n g the r e s o l u t i o n of the c r o s s c o r r e l a t e d v i b r o s e i s t r a c e . Various approaches have been proposed to s o l v e t h i s problem. Klauder et a l . (1960) weighted the spectrum i n order to r e s o l v e radar s i g n a l s . Gurbuz (1972) a p p l i e d t a p e r i n g t o reduce the s i d e lobes of v i b r o s e i s wavelets. While these approaches are e f f e c t i v e , only the v i b r o s e i s t r a c e spectrum w i t h i n a c e r t a i n freguency band i s d e a l t with. The proposed 1 3 7 d e c o n v o l u t i o n method deals not only with the spectrum of the t r a c e w i t h i n the bandpass on the input s i g n a l , but a l s o makes the optimum assumption about the spectrum o u t s i d e t h i s freguency band. The d e c o n v o l u t i o n method d e s c r i b e d here has the f o l l o w i n g o b j e c t i v e s , 1, To remove the e f f e c t s of a t t e n u a t i o n and d i s p e r s i o n on the v i b r o s e i s spectrum i n order to o b t a i n an estimate of the impulse response's spectrum wi t h i n the band of the FM sweep. 2 . To i n c r e a s e the r e s o l u t i o n of the v i b r o s e i s t r a c e by extending the F o u r i e r transform of the impulse response beyond the l i m i t s of the FM sweep. 5.4 The Decpnvplution Method-In d e a l i n g with the f i r s t o b j e c t i v e , c e p s t r a l a n a l y s i s i s employed. In doing so, we assume t h a t the spectrum of e (t) c o n t r i b u t e s to the low guefrency part of the spectrum of x ( t ) . Using the previous model f o r x (t) i n (5.6), and denoting the v i b r o s e i s wavelet as v(t) = w ( t ) * e ( t ) , we can write l o g ( X ( u r ) ) 138 ioa X(ur) = loa V(u) + loo Ifo) (5.7) where V(ar) anil I (to-) denote the Fourier transforms of v(f) and i ( t ) respectively. Hence, as shown before, the cepstrum of the trace, x (t) , can be rewritten as: where v ft) i s the cepstrum of the vibroseis wavelet and i (t) i s the cepstrum of the impulse response. In using the low guefrencies of x(t) to estimate the wavelet, log (X (<*r) ) i s e s s e n t i a l l y considered as a s i g n a l and the low "freguency" content of t h i s signal i s used to produce an estimate of log ( V ( u r ) ) . The low guefrency portion of log (X (or)) can be extracted by using x(t) values for | t|<T0 , where T Q i s the "cutoff quefrency". As shown by Clayton and Wiggins (1976), computation of the phase curve presents a major problem when computing the complex cepstrum. Since the Klauder wavelet i s zero phase and since phase unwrapping of noisy data i s unstable, the vibroseis wavelet i s assumed to be zero phase. The s e n s i t i v i t y of the method to t h i s phase assumption i s tested l a t e r . It should be noted that Wood (1976) has also used the zero phase assumption i n applying Wiener deconvolution f i l t e r s to symmetrical pulse (5.8) shapes. In the f i r s t step of the deconvolution process, the amplitude spectrum i s flattened by cepstral f i l t e r i n g of the zero phase cepstrum while the phase i s l e f t untouched. The zero phase cepstrum of the trace i s given by 2T I ^ T X 0 (t) = j /©j I X Cur) I e' dur (5>9) This cepstrum can also be written as Xo Ct) = ~± cfur (5.10) The f i r s t term i n (5.10) represents the log amplitude spectrum, or zero phase cepstrum, of the wavelet. Since the cepstrum of the wavelet generally contributes to the low guefrency content of the cepstrum, the f i r s t term i n (5.10) can be estimated by applying a low cut l i f t e r to the cepstrum. JI(u>)| can be estimated by applying the Fourier transform and exponential operations to the high cut portion of the cepstrum. A second approach for estimating l l f o r H requires that |V (or) | be estimated by applying the Fourier transform and exponential operations to the low cut portion of the cepstrum. Subsequent d i v i s i o n of \X(u)\ by | V (w) | then yields j I (<j) ]. Although these two approaches are mathematically equivalent and 1 4 0 give s i m i l a r r e s u l t s for synthetics, the l a t t e r approach was used when analyzing r e a l data since the problem of high guefrency noise i s avoided. Also, one deals d i r e c t l y with the amplitude spectrum of the vibroseis wavelet, which i s easier to deal with than the amplitude spectrum of the impulse response. Assuming the phase of the trace, &iUf), • - fo be the same as the phase of the impulse response allows the impulse response to be computed by the following inverse Fourier transform. Ut) - Tw J K w ) / e e dur ( 5 - 1 1 ) The impulse response which r e s u l t s from this type of cepstral f i l t e r i n g generally has a spectrum which i s a flattened version of the spectrum of the o r i g i n a l trace. Nevertheless, the fact that the spectrum i s bandlimited causes r e f l e c t e d pulses in i ( t ) to have disturbing side lobes. This leads to the second objective of the proposed deconvolution method, which i s that of extending the Fourier transform of the impulse response. Conventional data processing would leave zeros i n the spectrum for freguencies outside the range of tne input FM sweep. However, such a spectral model i s not i n accordance with the continuous nature of the impulse response spectrum, nor i s i t i n agreement with the pr i n c i p l e of maximum entropy. Furthermore, the Paley-fliener condition reguires that spectral values of stochastic processes be nonzero throughout the entire 141 frequency range (Smylie et a l . (1973)). Let us b r i e f l y review the maximum entropy p r i n c i p l e . This p r i n c i p l e reguires that the entropy or information content of a process be maximized when computing a power spectrum under the constraint that the spectrum be consistent with available data. The maximum entropy approach can be related to autoregressive modelling, and i t represents the optimum approach for c a l c u l a t i n g power spectra of f i n i t e data segments. In order to place nonzero values i n the impulse response spectrum between zero and the Nyquist frequency, an approach i s used which i s consistent with the p r i n c i p l e of maximum entropy. This approach, which has been successfully applied to teleseismic data by Clayton and Wiggins (1976), uses estimates of the impulse response's Fourier transform within the frequency range of the FM sweep i n order to predict Fourier transform values outside the range of the sweep. The prediction operator i s determined by use of the complex Burg algorithm, as described by Smylie et a l . (1973) and Burg (1975). Use of maximum entropy prediction provides a Fourier transform which i s consistent with known transform values within the frequency range of the input sweep while.remaining maximally noncommital with respect to unknown values. Ulrych and Bishop (1975) establish t h i s property of the maximum entropy method (MEM) and also point out that use of MEM i n designing a prediction error operator i s equivalent to using an autoregressive model for the predicted sequence. Use of autoregressive prediction allows us to f i l l in unknown parts of the spectrum by use of known 142 information instead of using the incorrect assumption of leaving zeros in the spectrum. The use of autoregressive prediction to resolve a single a r r i v a l i s shown i n Figure 31. In deconvolution, one wishes to remove t h i s wavelet and be l e f t with a delta function. In other words, id e a l s p e c tral prediction would give a completely f l a t spectrum between 0 hz and the Nyguist frequency. Figures 31 (b) and 31 (c) show spectra and deconvolutions re s u l t i n g from dif f e r e n t predictions of the Fourier transform of i ( t ) . Prediction using the l i m i t s of the FM sweep gives s l i g h t l y better resolution of the a r r i v a l while prediction using the transform values between 20 and 50 hz yields a re s u l t which i s e s s e n t i a l l y the desired delta function. The length of the autoregressive prediction f i l t e r had a length equal to half the length of the known transform. In our autoreqrassive modal, the impulse response of the earth i s represented as a discrete sum of delta functions. (5. 12) where ^ i s the sampling i n t e r v a l i n the time domain. 143 n i i 62.5 125.0 187.5 250.0 FREQUENCY 0.2 0.B 1.0 0.0 0.6 F i g u r e 31a) Example of a Klauder wavelet and i t s amplitude spectrum. The wavelet i s the a u t o c o r r e l a t i o n of an FM sweep g e n e r a t i n g f r e g u e n c i e s i n the 14-56 hz range. 144 Figure 31b). The spectrum aft e r applying an autoregressive predictor designed from the Fourier transform i n the 14-56 hz range. Figure 31c) shows the spectrum and wavelet resulting from the use of the predictor designed from within the 20-50 hz range. to SPECTRUM OF DECON JINIII 0.0 62.5 125~D 187.5 2S0.0 FREQUENCY DECONVOLUTION 0.0 0.2 0.4 T 0.6 0.8 o SPECTRUM OF DECON 1 1 0.0 62.5 125.0 FREQUENCY -T 1 187.5 250.0 DECONVOLUTION 146 The Fourier transform, I (<«>•) , i s then a sum of complex sinusoids In predicting values of I ( (*r ) beyond the frequency l i m i t s of the vibroseis sweep, the assumption i s made that a sum of sinusoids can be represented as an autoregressive process. Ulrych et a l (1973) and Ulrych and Clayton (1976) give examples which agree with t h i s assumption. A discussion on the use of autoregressive models to represent harmonic processes with additive noise has been presented by Ulrych and Clayton (1976). An example of a series of three delta functions with 1% additive noise i s given in Figure 32(a) along with i t s corresponding amplitude spectrum. Figure 32 (b) shows that bandlimiting the output of the delta functions produces side lobes on the a r r i v a l s . As shown in Figure 32(c), prediction of the bandlimited Fourier transform produces a Fourier transform very similar to the o r i g i n a l transform of Figure 32(a). This extension of the trace's Fourier transform provides an excellent deconvolution. (5.13) K 147 a) f— 0.0 AMPLITUDE SPECTRUM 1 w 1 1 33.25 62.5 FREQUENCY 93.75 ~1 125.0 INPUT SEQUENCE Figure 32a). The input spike seguence and i t s spectrum are shown. The spike sequence i s a series of three delta functions with 1% additive noise. I t s Fourier transform i s a sum of noisy sinusoids. b) M O D I F I E D T R A C E SPECTRUM T 0 . 0 31 .25 6 2 . 5 F R E Q U E N C Y 93 .75 "1 125.0 M O D I F I E D T R A C E 1 T I 2 U M E 0.4 Figure 32b). Spectrum which i s l i m i t e d to the 14-56 range produces the modified trace shown. 149 c i MEM PRED S P E C o.o ~i 1 31.25 62.5 FREGUENCY ~1 1 93.75 125.0 DECONVOLUT ION T 7 ~ TIME 0.4 Figure 32c). A prediction operator of length 10 designed on the bandlimited Fourier transform produces the spectrum shown and the re s u l t i n g deconvolution. 150 Let us consider the design of the prediction operator i t s e l f . If $ UT i s the frequency domain sampling i n t e r v a l , an rath order autoregressive representation of the discrete Fourier transform of the impulse response would have the following form. where I ( n S u r ) represents values of the discrete Fourier transform of i (t) , ot(r\$u>) i s the complex Burg prediction operator, and N ( n £w) 1 S a white noise seguence representing the prediction error or innovation of the model. In applying the Burg algorithm, the prediction error f i l t e r (1, - c < { $OJ) ,-c^{2 Sur) oC (m $(*>)) , and the complex conjugate of i t s time reverse are convolved with known values of I (n£« ) to y i e l d the prediction error and hindsight error sequences respectively. Values of o^(n £«•>-) are determined by minimizing the sum of the prediction and hindsight error powers and applying the Levinson algorithm. As pointed out by Claerbout (1976), t h i s produces a minimum phase prediction error operator, implying that predicted energy w i l l not increase outside the known passband. It i s important to note that the prediction f i l t e r i s designed over that portion of the spectrum f o r which the signal to noise r a t i o i s large. As shown by Figure 31, this portion w i l l be inside the l i m i t s of the FM sweep. Fortunately, the vibroseis technique allows us to know generally which portion of 151 the spectrum c o n t a i n s the s i g n a l . For p r e d i c t i o n , i t i s important to know what l e n g t h of the p r e d i c t i o n f i l t e r should be used, or e g u i v a l e n t l y , what the order of the a u t o r e g r e s s i v e process should be. A means of determining a s u i t a b l e order f o r the a u t o r e g r e s s i v e process has been presented by Akaike (1969). Akaike's f i n a l p r e d i c t i o n e r r o r (FPE) c r i t e r i o n determines the a u t o r e g r e s s i v e order by minimizing the sum of the p r e d i c t i o n e r r o r power and the e r r o r power i n v o l v e d i n e s t i m a t i n g the p r e d i c t i o n f i l t e r c o e f f i c i e n t s . For harmonic processes, however, the FPE c r i t e r i o n may not always be r e l i a b l e ( U l r y c h and C l a y t o n (1976)) and the l e n g t h of the p r e d i c t i o n operator may have to be chosen i n a s u b j e c t i v e manner. The p r e d i c t i o n f i l t e r l e n g t h determines the amount of p r e d i c t e d energy placed i n t o unknown p a r t s of the spectrum, and hence c o n t r o l s the degree of d e c o n v o l u t i o n . 5.5 Jk££iications to S y n t h e t i c and A c t u a l Data Let us c o n s i d e r s y n t h e t i c examples which i l l u s t r a t e the a p p l i c a t i o n of t h i s d e c o n v o l u t i o n method. F i g u r e 33(a) i l l u s t r a t e s a s e r i e s of s p i k e s whose spectrum i s m u l t i p l i e d by a model f o r |V(co^)'J. T h i s modal i s formed by the h a l f s of two d i f f e r e n t Gaussians and g i v e s a s i m i l a r t r a c e spectrum to sweep spectrum r a t i o as the example given by Ristow and Jurczyk (1975). A p e r f e c t l y f l a t spectrum i s used f o r the zero phase Klauder wavelet i n t h i s case. The e f f e c t of m u l t i p l y i n g 152 the trace spectrum of Figure 32 (a) by the wavelet spectrum, | V ( u r ) | , of Figure 33(a) i s shown by the modified trace spectrum and the modified trace i n Figure 33 (a). The problem i s to deconvolve the modified trace in order to recover the input spike series, since {V (or) \ i s a low guefrency phenomenon, a good estimate of j V ( w ) | can be obtained by using a low cut ceps t r a l f i l t e r . Low cut c e p s t r a l f i l t e r i n g estimates the wavelet spectrum from the amplitude spectrum of the trace as shown in Figure 33 (b). This wavelet spectrum contains the low guefrencies of the o r i g i n a l amplitude spectrum. As shown by Figure 33(c), the estimate of \V(ur)\ i s very similar to the actual value used i n construction of the synthetic shown in Figure 33(a). In t h i s case, a low cut f i l t e r of length 23, which was symmetrical about t=0, was applied to a cepstrum of length 256. The trace cepstrum i s then divided by the estimate of IV(itr)J to produce an estimate of \I(ur)\, Using the o r i g i n a l phase of the trace gives an estimate of the impulse response's Fourier transform, I (or), i n the frequency band from 18 to 50 hz. The unknown values of I (or) are f i l l e d i n by using an auto-regressive prediction f i l t e r of length 10. The spectrum r e s u l t i n g from cepstral f i l t e r i n g and autoregressive prediction i s also shown i n Figure 33 (c). As shown by Figure 33 (c), the ' re s u l t i n g deconvolution i s very s i m i l a r to the desired seguence of spikes. Figure 33a). The amplitude spectrum of a vibroseis wavelet i s shown along with the modified spectrum and trace. The amplitude spectrum of t h i s wavelet includes the e f f e c t of "earth f i l t e r i n g " . a) WAVELET SPECTRUM FREQUENCY 125.0 MODIFIED TRRCE SPECTRUM 0.0 1 r i 1 33.25 62.5 93.75 125,0 FREQUENCY MODIFIED TRRCE n I r~4 155 0.0 1 5 . 6 2 5 3 1 . 2 2 4 3 . 2 7 5 6 2 . 5 F R E Q U E N C Y Figure 33b). A comparison i s given between the spectrum of the modified trace and an estimate of the wavelet spectrum obtained by ce p s t r a l f i l t e r i n g . 156 Figure 33c). The spectrum r e s u l t i n g from c e p s t r a l f l a t t e n i n g and subseguent autoregressive prediction i s shown as well as the resulting deconvolution. The amplitude spectrum estimated by. cep s t r a l f i l t e r i n g i s very s i m i l a r to the o r i g i n a l wavelet spectrum shown i n Figure 33a). 157' C) RMP S P E C REMOVED 33.25 E2.5 93.75 FREQUENCY 125.0 MEM,PRED S P E C 0.0 T 1 r 3 1 . 2 5 P 2 . 5 9 3 . 7 5 FREQUENCY i 125.0 DECONVOLUTION 0.2 TIME • 9 1 0.4 0.6 158 d) . M O D I F I E D TRACE 0.0 D E C O N V O L U T I O N a w ^ t . . . 0.3 0 .4 TIME 0.6 Figure 33d), The association of a minimum phase spectrum with the wavelet amplitude spectrum gives the modified trace shown. Application of zero phase cepstral f i l t e r i n g and autoregressive prediction gives the resulting deconvolution. 159 The d e c o n v o l u t i o n t e s t i n F i g u r e 33 was performed using the zero phase assumption. Ristow and Jurczyk (1975) a p p l i e d wiener d e c o n v o l u t i o n and made minimum phase assumptions. Although the Klauder wavelet i s zero phase, e (t) i s o f t e n c o n s i d e r e d as a minimum phase f u n c t i o n . Let us examine the e f f e c t of having a minimum phase a s s o c i a t e d with the e f f e c t s of " e a r t h f i l t e r i n g " . The H i l b e r t transform of log(|E (w ) | ) g i v e s a minimum phase s h i f t a s s o c i a t e d with e ( t ) . T h i s produces the modified t r a c e shown i n F i g u r e 33(d). The a p p l i c a t i o n of the same d e c o n v o l u t i o n approach as i n F i g u r e 33 (a) s u p p l i e s a deconvolution which i s phase s h i f t e d from the r e s u l t of F i g u r e 33(c) but which d e f i n e s the d e s i r e d s p i k e s . F i g u r e 34 shows the r e s u l t of e s t i m a t i n g a t r u n c a t e d minimum phase wavelet from the modified t r a c e of F i g u r e 33(d), and then using a S i e n e r f i l t e r which shapes the wavelet to a s p i k e . Although the spectrum of the deconvolved t r a c e i s f l a t t e n e d , i t i s b a n d l i m i t e d . T h i s r e s u l t s i n a d e c o n v o l u t i o n which i s i n f e r i o r to that given i n F i g u r e 33 (d). There are two other reasons why c o n v e n t i o n a l minimum phase deconvolutions are r e s t r i c t e d i n t h e i r performance c a p a b i l i t i e s . F i r s t of a l l , the v i b r o s e i s wavelet i s e s s e n t i a l l y a f i l t e r e d Klauder wavelet, and t h e r e f o r e i t i s not minimum phase. Secondly, the d e s i r e d impulse response of the e a r t h may not approximate a white noise seguence. The l a t t e r assumption i s necessary i n minimum phase d e c o n v o l u t i o n approaches when e s t i m a t i n g the wavelet a u t o c o r r e l a t i o n from the t r a c e a u t o c o r r e l a t i o n . Figure 34 . Minimum phase deconvolution of the vibroseis trace model: a) Input trace i d e n t i c a l to that used i n Figure 33d). b) Truncation of minimum phase wavelet estimated using the Hilbert transform method. c) Spectrum of the deconvolution. d) Deconvolution r e s u l t i n g from the application of a Hiener inverse f i l t e r designed to spike the wavelet. WWRVELET ESTIMATE SPECTRUM OF DECONVOLUTION A 0.0 62.5 J25.0 FREQUENCY J67.S 250.0 DECONVOLUTION i 4 0.0 0 TJIH r •IA>«"»,''VW| 1 0.4 162 Having o b t a i n e d encouraging r e s u l t s from s y n t h e t i c data, the method of c e p s t r a l f i l t e r i n g and s p e c t r a l e x t e n s i o n was used with a s e t of v i b r o s e i s t r a c e s i n which a t l e a s t one ' b r i g h t spot' gas r e f l e c t o r i s l o c a t e d . The d e c o n v o l u t i o n was a p p l i e d at two stages i n the data processing seguence. F i g u r e 35 shows a d e c o n v o l u t i o n of the data a f t e r s p h e r i c a l divergence c o r r e c t i o n s , s t a t i c s removal, normal moveout c o r r e c t i o n s , and s t a c k i n g . F i g u r e 35(a) shows the i n p u t data which were obtained from and i n i t i a l l y processed by Chevron Standard. The d e c o n v o l u t i o n a t two d i f f e r e n t stages i s shown. In F i g u r e 35(b), c e p s t r a l f i l t e r i n g p rovides the deconvolution shown. Fi g u r e 35(c) d i s p l a y s the d e c o n v o l u t i o n a f t e r c e p s t r a l f i l t e r i n g and a u t o r e g r e s s i v e p r e d i c t i o n of the F o u r i e r transform. Each stage of the d e c o n v o l u t i o n i n c r e a s e s the r e s o l u t i o n of the v i b r o s e i s t r a c e at the expense of i n c r e a s i n g the high freguency n o i s e . Figure 35(d) compares the amplitude s p e c t r a of the t r a c e s at v a r i o u s stages of p r o c e s s i n g to show the frequency domain e f f e c t s of the d e c o n v o l u t i o n approach. A c e p s t r a l f i l t e r of l e n g t h 27 was used to f l a t t e n the spectrum, and a p r e d i c t i o n f i l t e r of l e n g t h 15 samples (or 60 msec)was designed using F o u r i e r transform values between 16 and 43 hz. I n s p e c t i o n of the amplitude s p e c t r a and deconvolved t r a c e s at each stage i n the d e c o n v o l u t i o n process provides a good means of determining s u i t a b l e p r o c e s s i n g parameters. Often, d i f f e r e n t s e t s of c e p s t r a l f i l t e r l e n g t h s and a u t o r e g r e s s i v e o r d e r s are used i n order to f i n d a s u i t a b l e d e c o n v o l u t i o n which provides a good compromise between r e s o l u t i o n and s t a b i l i t y . Figure 3 5 . a) Input stacked vibroseis data applied after NMO and s t a t i c s corrections. (b) F i r s t stage of deconvolution after cepstral f i l t e r i n g . (c) Second stage of deconvolution after cepstral f i l t e r i n g and autoregressive prediction of the Fourier transform. 165 Figure 35 d). Comparison of amplitude spectra; (i) before deconvolution, ( i i ) a f t e r cepstral f l a t t e n i n g , and ( i i i ) after autoregressive prediction. l ) S P E C T R A B E F O R E D E C O N V O L U T I O N i i ) A F T E R C E P S T R A L F I L T E R I N G 166 i i i ) A F T E R A U T O R E G R E S S I V E P R E D I C T I O N rREOUtXCT A A . A A A A A ^ A _ . A . A A J A . . A. 4 i J f r \ A j v l i ,»||.| ~\v 167 In Figure 36, the two step deconvolution process decreases the width of a r e f l e c t i o n event which was caused by the presence of natural gas. The location of the r e f l e c t i o n from the natural gas deposit i s shown by the arrows on the seismic sections of Figure 36. Figure 36(a) shows a section in which the 'bright spot' r e f l e c t i o n from a gas deposit was enhanced by Wiener deconvolution and amplitude sc a l i n g . The section in Figure 36(b) displays an additional•deconvolution of the section shown in Figure 36 (a). This deconvolution was applied to provide a better resolution of the "bright spot" r e f l e c t i o n . The deconvolution was obtained by using a cepstral f i l t e r of length 27 and a prediction f i l t e r of length 25* Deconvolution compresses the r e f l e c t i o n frcm the gas deposit to a si n g l e well defined event on the section to give a sharper d e f i n i t i o n of the desired horizon. Figure 36. a) F i n a l processed section with indicated "bright spot" emphasized by Wiener f i l t e r i n g and amplitude adjustments of data i n Figure 35. b) Subseguent deconvolution compresses pulse from "bright spot" r e f l e c t i o n shown i n Figure 36a). As i n Figure 35, the variable area - wiggle plo t has i t s p o l a r i t y reversed to emphasize the negative r e f l e c t i o n c o e f f i c i e n t r e s u l t i n g from the gas deposit. 170 CHAPTER 6 SOgflaSI This thesis investigates several methods of deconvolving r e f l e c t i o n recordings i n exploration seismology and makes recommendations as to the usage of these methods. Chapter 1 presents models for the impulse response of a layered earth and wavelets generated by explosive sources. The seismic trace i s considered to be,a convolution of the earth's impulse response with a source generated wavelet. The purpose of deconvolution i s to remove the wavelet and y i e l d the impulse-response of a layered earth. The problem of deconvolving explosion generated seismic data i s investigated in Chapters 2 to 4. Three wavelet estimation methods are considered. Wiener-Levinson wavelet estimation i s a minimum phase time domain approach which reguires an estimate of the wavelet length. The Wold-Kolmogorov estimation method i s a frequency domain approach which uses the same assumptions as the Wiener-Levinson method. Wold-Kolmogorov f a c t o r i z a t i o n e s s e n t i a l l y uses a Hilbert transform of the log amplitude spectrum to y i e l d a minimum phase spectrum. The Wiener-Levinson and Wcld-Kolmogorov methods provide s i m i l a r wavelet estimates for the synthetic example used. The homomorphic deconvolution method provides a wavelet estimate by separating the complex cepstrum of the wavelet from that of the trace. For noiseless synthetics, homomorphic deconvolution i s extremely promising when compared to the minimum phase approaches. Problems with phase curve unwrapping a r i s e with t h i s method when additive noise i s present i n the seismic trace. 171 Thus some form of cepstral stacking i s useful when estimating wavelets from noisy data. On actual explosion data, the wavelets obtained by homomorphic deconvolution are c l e a r l y not minimum phase, causing the minimum phase assumption to remain suspect. After estimation of a seismic wavelet, the problem of designing single channel and multichannel Wiener deconvolution f i l t e r s i s e s s e n t i a l l y one of choosing a suitable f i l t e r length and a l e v e l of prewhitening. Although the multichannel Wiener deconvolution f i l t e r i s superior to the single channel f i l t e r when the wavelet i s known, the differences i n the methods becomes much less s i g n i f i c a n t when estimated wavelets and noisy data are used. Differences i n the Hiener and MEM deconvolution f i l t e r s r e s u l t from the d i f f e r e n t autocorrelations being used i n their design. While the difference i n autocorrelation estimates becomes i n s i g n i f i c a n t for t y p i c a l gate lengths used i n seismic processing, the MEM approach for cal c u l a t i n g a time adaptive f i l t e r i s preferable when using short data segments. Deconvolution of noisy dynamite generated data i s presented in Chapter 4 . Although the deconvolution results were very s i m i l a r to those obtained by another independent research group, problems are encountered when attempting to correlate the impulse response synthetics and seismic data. The reasons for these problems l i e both in the data quality and the inherent differences between what the well logs and the seismic r e f l e c t i o n s are measuring. 172 In Chapter 5 , the properties of vibroseis recordings are exploited to develop a new approach to deconvolution which uses cepstral f i l t e r i n g and autoregressive prediction. This approach uses the assumption of a low guefrency, zero phase wavelet. The zero phase assumption avoids the problem of phase unwrapping which i s often encountered when using homomorphic deconvolution on r e a l data, and i s j u s t i f i e d by the f a c t that the vibroseis trace i s composed of f i l t e r e d Klauder wavelets. Even when the earth f i l t e r , e ( t ) , i s minimum phase, the deconvolution method i s successful. Cepstral f i l t e r i n g removes the low guefrency contribution of the wavelet to the spectrum within a r e s t r i c t e d passband. An autoregressive prediction f i l t e r i s then used to extend the Fourier transform of the impulse response beyond the l i m i t s of the FM sweep. This spectral extension scheme i s based on the idea that the Fourier transform of a sum of delta functions can be predicted by using an autoregressive model. Since the bandwidth of the vibroseis signal i s well known from a knowledge of sweep frequencies, i t i s not d i f f i c u l t to decide which portion of the spectrum i s to be used in designinq the prediction f i l t e r . As shown by the synthetic and r e a l data analysed, the deconvolution provides qood resolution of vibroseis traces. Since the method does not reguire the assumption of a random impulse response, i t can be used on short segments of the seismic trace. Also, because of the high resolution p o s s i b i l i t i e s of the method, i t can be p a r t i c u l a r l y useful in resolving vibroseis r e f l e c t i o n s from thin beds. 173 REFERENCES Akaike,H., 1969, F i t t i n g a u t o r e g r e s s i v e models f o r p r e d i c t i o n : Ann. I n s t . S t a t i s t . Math. ,v.21, p.243-247. Andersen, N. , 1974, On the c a l c u l a t i o n of f i l t e r c o e f f i c i e n t s f o r maximum entropy s p e c t r a l a n a l y s i s : Geophysics, v.39, p.69-72. Berryman, L.H., G o u p i l l a u d , P.L., and Waters, K.H., 1 958, R e f l e c t i o n s from m u l t i p l e t r a n s i t i o n l a y e r s : Geophysics, v.23, p.223-243. Bogert, R.P., Healy, M.J., and lukey,J.W., 1963, The guefrency a l a n y s i s of time s e r i e s f o r echoes: Proc. Symp. On Time S e r i e s A n a l y s i s , M. Rosenblatt, ed., New York, S i l e y , p.209-243. Bracewell,R., 1965, The F o u r i e r transform and i t s a p p l i c a t i o n s : New York, McGraw - H i l l . Buhl,P., S t o f f a , P . L . . And Bryan,G.M., 1974, The a p p l i c a t i o n of homomorphic deco n v o l u t i o n to shallow - water marine seismology - pa r t 2: r e a l data: Geophysics, v.39, p. 417-426. Burg,J.P., 1967, Maximum entropy s p e c t r a l a n a l y s i s : paper presented a t the 1967 SEG meeting i n Oklahoma C i t y , Okla. Burg,J.P., 1975, Maximum entropy s p e c t r a l a n a l y s i s : PhD t h e s i s , S t a n f o r d Univ., Palo A l t o , C a l i f . B uttkus, B., Homomorphic f i l t e r i n g : Geophys. Prosp., v.23, p. 723-748. C l a e r b o u t , J . , 1976, Fundamentals of g e o p h y s i c a l data p r o c e s s i n g : New York, McGraw-Hill. Clarke,G.K.C. 1968, Time v a r y i n g d e c o n v o l u t i o n f i l t e r s : Geophysics, v.33, p.936-944. C l a y t o n , R.W., and U l r y c h , T . J . , 1976, R e s t o r a t i o n of impu l s i v e f u n c t i o n s : IEEE Trans. I n f o . Theory, submitted f o r p u b l i c a t i o n . Clayton,R.W., McClary,B., and l i g g i n s , R . A . , 1975, Comments on •Phase d i s t o r t i o n and H i l b e r t t r a n s f o r m a t i o n i n m u l t i p l y r e f l e c t e d and r e f r a c t e d body waves' by G.L. Choy and P.G. R i c h a r d s : B u l l . Seism. Soc. Am., i n p r e s s . 174 Clayton,R.W., and Wiggins, R.A., 1976, Source shape estimation and deconvolution of teleseismic bodywaves: Geophys. J.R. Astr. S o c , i n press. Clowes, R.M., 1969, Seismic r e f l e c t i o n investigations of crustal structure i n Southern Alberta: PhD t h e s i s , Univ. Of Alberta, Edmonton, Aita. Crump,N., 1974 , A Kalman f i l t e r approach to the deconvolution of seismic signals: Geophysics, v.39 , p.1 - 1 3 . Davies,E,B., and Mercado, E.J., 1968, Multichannel deconvolution f i l t e r i n g of f i e l d recorded seismic data: Geophysics, v.33 , p.711-722 . Davies, J.C., 1976 , Maximum entropy spectral analysis of free o s c i l l a t i o n s of the earth: M.Sc. Thesis, University of B r i t i s h Columbia. DeVoogd, N., 1974, Wavelet shaping and noise reduction: Geophys. Prosp,, v . 2 2 , p.354 -369 . Fryer,G.J., Odegard,M.E., and Sutton,G.H., 1975, Deconvolution and spectral estimation using f i n a l prediction error: Geophysics, v.40, p.411-425 . Grant,F., and West,G., 1965 , Interpretation theory i n applied geophysics: Toronto, McGraw - H i l l . Futterman, W.I., 1962 , Dispersive body waves: J.Geoph.Res., v . 6 7 , p. 5 2 7 9 - 5 2 9 1 . Galbraith,J.N., 1971 , Prediction error as a c r i t e r i o n for operator length: Geophysics, v.35 , p.251-265 . Gurbuz, B.M., 1972 , Signal enhancement of vibratory source data i n the presence of attenuation: Geophys. Prosp., v.20 , p.421-438 . Jenkins,G.M. And Watts, D.G., 1969, Spectral Analysis and i t s applications: San Francisco, Holden-Day. Kanasewich, E.R., 1973 , Time sequence analysis i n geophysics: University of Alberta Press, Edmonton, Alberta. 175 Klauder,J,R., P r i c e , A.C., D a r l i n g t o n , S., and Albersheim, W.J.f 1960, The theory and design o f c h i r p radars; B a l l System Tech. Jou r . , v.39, p. 145-807. Levin s o n , N., 1947, The Wiener RMS (root mean square) e r r o r c r i t e r i o n i n f i l t e r design and p r e d i c t i o n : J o u r . Of Math and P h y s i c s , v.25, p.261-278. L i n e s , L.R., 1974, A note on the a p p l i c a t i o n of Wiener multichannel d e c o n v o l u t i o n : J . Can. Soc. E x p l . G e o p h y s i c i s t s , v.10, p.65-70. L i n e s , L.R. And U l r y c h , T.J., 1976, Technigues f o r s e i s m i c d e c o n v o l u t i o n and wavelet e s t i m a t i o n ; paper presented a t the 1975 SEG meeting i n Denver, Colorado and accepted f o r p u b l i c a t i o n by Geophys. Prosp. Pending r e v i s i o n . L i n e s , L.R. And C l a y t o n , R.W., 1976, A new approach f o r v i b r o s e i s d e c o n v o l u t i o n : paper presented at the 1976 CSEG convention i n C a l g a r y , A l b e r t a , and submitted f o r p u b l i c a t i o n t o Geophys. Prosp. Oppenheim, A.V., 1965, S u p e r p o s i t i o n i n a c l a s s of n o n - l i n e a r systems; Research Lab of E l e c t r o n i c s , MIT, Tech. Rep. 432. . Oppenheim, A.V., S c h a f e r , R.W., and Stockham, T.G., 1968, Nonlinear f i l t e r i n g of m u l t i p l i e d and convolved s i g n a l s : Proc. IEEE, v.65, p.1264-1291. O t i s , S.M. And Smith, R.B., 1974, Homomorphic d e c o n v o l u t i o n by l o g s p e c t r a l averaging: paper presented a t the 1974 SEG meeting i n D a l l a s , Texas. Peacock, K.L., and T r e i t e l , S., 1969, P r e d i c t i v e d e c o n v o l u t i o n : theory and p r a c t i c e : Geophysics, v.34, p.155-169. R i c k e r , N., 1953, The form and laws of propogation o f s e i s m i c wavelets: Geophysics, v.18, p.10-40. R i l e y , D. And Burg, J.P., 1972, Time and space ada p t i v e d e c o n v o l u t i o n f i l t e r s : paper presented a t the 42nd meeting of the SEG, Anaheim, C a l i f . , Nov. 28, 1972. Ristow, D., and J u r c z y k , D., 1974, V i b r o s e i s d e c o n v o l u t i o n : Geophys. Prosp., v.23, p.363-379. Robinson, E,A., 1967a, M u l t i c h a n n e l time s e r i e s a n a l y s i s with d i g i t a l computer programs: San F r a n c i s c o , Holden-Day. 176 Robinson, E.A., 1967b, Predictive decomposition of time series with application to seismic exploration: Geophysics, v.32, p.418-484. Robinson, E.A., and T r e i t e l , S., 1967, Princip l e s of d i g i t a l Wiener f i l t e r i n g : Geophys. Prosp., v.15, p.311-333. Robinson, E.A., 1975a, Deconvolution methods and seismic models: paper presented at the 1S75 SEG meeting i n Denver, Colorado. Robinson, E.A., 1975b, Dynamic predictive deconvolution: Geophys. Prosp., v.23, p.780-798. Schafer, R.W., 1969, Echo removal by discrete generalized l i n e a r f i l t e r i n g : Research Lab of Electronics, MIT, Tech Rep, 466. Smylie, D,E., Clarke, G.K.C., and Ulrych, T.J., 1973, Analysis of i r r e g u l a r i t i e s i n the earth's rotation: Methods of Comp. Physics, v.13, Acad. Press Inc., New York. Stoffa, P., Buhl, P., and Bryan, G., 1974, The application of homomorphic deconvolution to shallow water marine seismology - part 1: theory: Geophysics, v.39, p.417-436. T r e i t e l , S., and Robinson, E.A., 1966, The design of high resolution d i g i t a l f i l t e r s : IEEE Trans. On Geosci. Electronics, v.GE-4, p25-38. T r e i t e l , S., 1970, Princip l e s of d i g i t a l multichannel f i l t e r i n g : Geophysics, v.35, p.785-811. T r e i t e l , S., 1975, An ide a l performance property of the multichannel Wiener f i l t e r : J. Can. Soc. Of Expl. Gaophysicists, v.10, p.71-74. Ulrych,T.J., 1971, Application of homomorphic deconvolution to seismology: Geophysics, v,36, p.650-660. Ulrych, T.J., 1972, Maximum entropy power spectrum of truncated sinusoids: J. Geophys. Res., v.77, p.1396-1400. Ulrych, T.J., Jensen, O.J., E l l i s , R.M., and Somerville, P. G. , 1972, Homomorphic deconvclution of some teleseismic events: B u l l . Seism. Soc. Am., v.62, p.1269-1281. 177 Ulrych, T.J., Smylie, B,E., Jensen, O.G., and Clarke, G. K.C,, 1973, Predictive f i l t e r i n g and smoothing of short records by using maximum entropy: J.Geophys.fies., v.77, p.1396-1400. Ulrych, T.J. And Bishop, T., 1975, Maximum entropy spectral analysis and autoregressive decomposition: Reviews of Geophysics, v.13, p.183-200. Ulrych, T.J,, and Clayton, R.W., 1976, Time series modelling and maximum entropy: Phys. Earth Planet. Int., i n press. Wang, R., 1969, The determination of optimum gate lengths f o r time varying Wiener f i l t e r i n g : Geophysics, v.34, p.683-695, White, R.E., and O'Brien, P.U.S., 1974, Estimation of the primary seismic pulse: Geophys. Prosp., v.22, p.627-651. Wiggins, R.A., and Robinson, E.A., 1965, Recursive solution to the multichannel f i l t e r i n g problem: J. Geoph Res., v,70, p.1885-1891. Wiggins, R.A,, 1968, The use of expected error i n the design of least squares optimum f i l t e r s , Geophys. Prosp., v.15, p.288-296. Wood, L.C. And T r e i t e l , S,, 1975, Seismic s i q n a l processing: Proc. IEEE, v.€3, p.649-661. Wood, L.C, 1976, High resolution deconvolution of symmetrical pulses using Wiener f i l t e r s : paper presented at the 1976 CSEG convention in Calgary, Alberta. Yedlin, M., 1975, Kalman f i l t e r i n g and the deconvolution .problem: unpublished manuscript. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052513/manifest

Comment

Related Items