UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The deconvolution of teleseismic recordings Clayton, Robert W. 1975

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_1975_A6_7 C53.pdf [ 3.4MB ]
Metadata
JSON: 831-1.0052880.json
JSON-LD: 831-1.0052880-ld.json
RDF/XML (Pretty): 831-1.0052880-rdf.xml
RDF/JSON: 831-1.0052880-rdf.json
Turtle: 831-1.0052880-turtle.txt
N-Triples: 831-1.0052880-rdf-ntriples.txt
Original Record: 831-1.0052880-source.json
Full Text
831-1.0052880-fulltext.txt
Citation
831-1.0052880.ris

Full Text

THE DECONVOLUTION OF TELES £IS MIC RECORDINGS by Robert H. Clayton B.A.Sc, University cf Toronto, 1973 A THESIS SUBMITTED IK PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n the Department of Geophysics and Astronomy We accept t h i s thesis as conforming to the required standard The Unversity of B r i t i s h Columbia July, 1975 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t o f t h e r e q u i r e m e n t s f o r a n a d v a n c e d d e g r e e a t t h e U n i v e r s i t y o f B r i t i s h C o l u m b i a , I a g r e e t h a t t h e L i b r a r y s h a l l m a k e i t f r e e l y a v a i l a b l e f o r r e f e r e n c e a n d s t u d y . I f u r t h e r a g r e e t h a t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g o f t h i s t h e s i s f o r s c h o l a r l y p u r p o s e s may be g r a n t e d by t h e H e a d o f my D e p a r t m e n t o r by h i s r e p r e s e n t a t i v e s . I t i s u n d e r s t o o d t h a t c o p y i n g o r p u b l i c a t i o n o f t h i s t h e s i s f o r f i n a n c i a l g a i n s h a l l n o t be a l l o w e d w i t h o u t my w r i t t e n p e r m i s s i o n . D e p a r t m e n t T h e U n i v e r s i t y o f B r i t i s h C o l u m b i a 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Abstract A method i s presented for the deconvolution of a suite of teleseismic recordings of tha same event. The suite must be r e s t r i c t e d in azimuthal and delta ranges i n order that the assumption of source s t a t i o n a r i t y can be made. The source i s estimated by averaging the log amplitude and phase spectra of the recordings. This method of source estimation uses the redundant source information contained in secondary a r r i v a l s . The necessary condition for this estimator to resolve the source wavelet i s that the t r a v e l times of the various secondary a r r i v a l s be evenly d i s t r i b u t e d with respect to the i n i t i a l a r r i v a l s . The actual deconvolution of the seismograms i s done by spectral d i v i s i o n with two modifications. The f i r s t i s the introduction of a ioiniirum allowable source spectral amplitude termed the waterlevel (Helmberger and Wiggins, 1971). This parameter constrains the gain of tha deconvolution f i l t e r i n regions where the seismcgrarc has • l i t t l e or no infor nation, and also trades-off a r r i v a l time resolution with a r r i v a l amplitude resolution. The second modification, designed to increase the time domain resolution of the deconvolution, i s the extension of the transmission path impulse response spectrum beyond i t s optimal passband (the passband of the seismcgrams). The i i j u s t i f i c a t i o n for the extension l i e s in the fact that the impulse response i s "impulsive" by nature which means i t s spectrum i s not band-limited. Thus, the impulse response i s best represented by a continuous spectrum rather than one which i s set to zero outside the optimal passband. The actual prediction i s dene by an recursive a p p l i c a t i o n of a unit-step prediction operator determined by Burg's maximum entropy algorixhm (Burg, 1967). The envelopes of the deconvolution are used to detect the presence of phase s h i f t e d a r r i v a l s . The source estimator and deconvolution method are applied to three examples: Kern Ccuntv (16/10/62) , West Au s t r a l i a (24/3/70) , and Novaya Zemlya (14/10/69). Source estimation by hcmcmcrphic transforms (Ulrych, 197 1) i s also discussed i n this t h e s i s . This method cf source estimation was found to be cf limited use because of the i n v a l i d i t y of the low quafrency assumption that must be applied (Ulrych, 1971), and the phase i n s t a b i l i t i e s cf the transform i t s e l f . i i i Table of Contents I. Introduction .......1 II. Source Estimations ................................ 4 1. General Assumptions .........4 2. Source Estimation by Homomorphic Transforms ....8 3. Source Estimation by Averaged log Amplitude and Phase Spectra ..17 III. Deconvoluticn of the Seismcgratns 24 1. The Waterlevel Parameter 25 2. Extension cf the Impulse Response Spectrum .....29 3. Envelopes of Deconvolution 35 17. Examples Of Deconvolution ......................... 38 1. Kern County, 16/10/62 ,.40 2. West A u s t r a l i a 24/3/70 52 3. Novaya Zemlya 14/10/69 63 V. Summary 75 References ...78 i v L i s t of Figures 1. Example of a Source a r d Impulse Response C e p s t r u a .9 2. Average of A l l West A u s t r a l i a Amplitude Spectra ...10 3. E f f e c t of Additive Noise on Amplitude and Phase Spectra 12 4. E f f e c t of Changing a Zero Location on the Phase ...13 5. Prediction of a Single Sinusoid 31 6. Example of F i l t e r e d and U n f i l t e r e d Impulse Responses '. 33 7. Examples of the Envelopes of deconvolution ........36 8. Suite of Kern County Seismcgrams 41 9. Suite of West Au s t r a l i a Seismograros 53 10. Suite of Ncvaya Zemlya Seismcgrams 64 V L i s t c f Tables 1. Kern county Eeismogram Suite 40 2. West Aust r a l i a Seistnogram Suite 52 3. Novaya Zemlya Seismcgrair Suite ..63 v i I would l i k e to take t h i s opportunity to thank the people who contributed to the work on this thesis. My supervisor, Prof. Ralph A. Wiggins, provided the necessary academic advice, moral encouragement, and f i n a n c i a l support to sea t h i s thesis through tc c c i p H e t i c r . Several of his computer programs were used in t h i s thesis, and the data for the examples was supplied ty him. I was also greatly assisted by several discussions I had with Prof. Tad J. Ulrych on homomcrphic transfcirrs ard the maximum entropy method. F i n a l l y , the thesis was influenced by the work cf Paul Somerville and Larry Lines, who are a c t i v e l j using hcmo-morphic transforms in the i r f i e l d s of research. 1 I. INTRODUCTION Much of tha research in seismology has teen directed toward the f i t t i n g cf earth seismic velocity models from the body wave information contained in teleseismic recordings. In p a r t i c u l a r , the a r r i v a l times and amplitudes usually form the observation set for the inversion to the ve l o c i t y models. The observations, however, are often d i f f i c u l t to measure because they are obscured by source complications, such as long source coda lengths, which cause some a r r i v a l s to overlap i n time. Tha process of removing the source time function component from a seismcgram i s termed deconvolution because i n the forward case the source function i s combined with the earth's transmission path impulse response under the operation of convolution. The purpose cf t h i s thesis i s to examine the problem cf the deconvolution of teleseismic recordings. The f i r s t step i n any deconvolution process i s the estimation of the source time function. Chapter II cf t h i s t hesis deals with t h i s trcblem, A standard technique i n exploration seismology i s to decompose the source from the seismogram autocorrelation function. This method assumes that the source i s minimum phase and that the iirpulse sequence behaves l i k e white noise ( Robinson, 1967 ). Since neither of these assumptions appeared to be generally v a l i d for teleseismic recordings, the technique was not used. Source estimation by homomorphic transformation has also been suggested ( Ulrych, 1971 ), and this method i s 2 discussed in some d e t a i l in Chapter I I . However, the low guefrency assumption of t h i s method was found net tc be generally v a l i d , and there i s also a problem of phase i n s t a b i l i t i e s with the homomorphic transform i t s e l f . Both of the above mentioned methods attempt to estimate the source from a sing l e seismogram. I decided, that on the basis of a sing l e recording of an event, i t would be d i f f i c u l t , i f not impossible, to devise a general method cf estimating the source time function. For example, i f an earthquake source has more than one d i s t i n c t motion, then i t is d i f f i c u l t to decide whether the e f f e c t as observed on a single seismogram belongs to the source or to the transmission path. Fcr t h i s reason, I decided to r e s t r i c t the problem to the case where a suite of seisroc grams from the same event i s a v a i l a b l e . The redundant source information contained i n the suite w i l l allow one tc estimate the source with greater confidence than i f the estimation were done on the basis of a single seismogram. Re s t r i c t i o n s have to be placed on the delta and a2imuthal ranges of the suite of seismograms so that approximate source s t a t i o n a r i t y can be assumed. The actual deconvolution of the seismograms with a given source estimate i s discussed i n Chapter I I I . Techniques such as Wiener deconvolution, which convolve the seismogram with an inverse cf the estimated source, were not used because of the problems encountered when the source i s not minimum phase ( Robinson, 1967 ). In the frequency 3 domain where deconvolution i s the d i v i s i o n of the seismogram by the source, the problems with the above mentioned methods are easily seen as the spectral locations where the source amplitude becomes small. The s p e c t r a l d i v i s i o n a l deconvolution method can be e a s i l y modified to overcome t h i s problem by constraining the minimum allowable source amplitude l e v e l { Helmberger and Wiggins, 1971 ). A second modification which allows increased time domain resolution of the impulse response i s also outlined i n Chapter III. In Chapter IV, the deconvolution of three r e a l examples which employs the nethcds outlined in Chapters II and III i s presented. I I . SOURCE ESTIMATIONS 1. General Assum£tions lor t e l s e i s m i c recordings there are no independent observations of tha source time function, so i t i s necessary that i t be estimated d i r e c t l y from the data. This, i n general, w i l l present some problems, because the source wavelet w i l l usually have seme c h a r a c t e r i s t i c s which overlap with those of the transmission path impulse response. Seme assumptions w i l l be necessary tc arb i t r a t e this problem of non-uniqueness. The general model that w i l l be considered for the form of each seismogram recorded from the same event i s : XJ (t) = s(t) ® hj (t) + nj (t) [ 2 - 1 ] where xj(t) is. tha time seri e s recorded at the j " 1 s t a t i o n , s(t) i s the source wavelet, hj(t) i s the transmission path impulse response to the j " 1 s t a t i o n , nj(t) i s a white noise sequence, and ® denotes convolution. The assumption contained in thi s model i s that the source i s stationary i n both time and space. Thus, within the azimuthal-takeoff angle cone defined by the spread cf the suite of seismograms, and by the time window applied to them, the source function i s independent of azimuth and 5 takeoff angle. In general, t h i s assumption i s not true, as the kinematic d i s l o c a t i o n models of earthquake radiation i n d i c a t e ( Savage, 1966 ), but i t s v a l i d i t y can be improved by r e s t r i c t i n g the azimuthal, d e l t a , and time range of the suite of seismograms. A p r a c t i c a l view point may be considered for t h i s assumption. If a source can be estimated that in seme way best represents the majcr.features of a l l the source e f f e c t s over the suite of seismograms, then i t s removal w i l l leave only the impulse response and seme hopefully small source differences. The procedure usually produces seme v i s u a l enhancement of the recordings, which i s one of the purposes of deconvolution, even though the r e s u l t may not be a true representation of the earth's impulse response. The assumption mentioned above w i l l not be s u f f i c i e n t to solve equation [2-1] for the source wavelet. Additional assumptions about the features of either the impulse response or the source wavelet w i l l be necessary. Cne common assumption that i s made i n exploration seismology i s that the source i s minimum phase. This allows a unique decomposition of the source from i t s autocovariance function ( Robinson, 1967 ). I f thi s assumption i s coupled with the assumption that the impulse response behaves as a white noise sequence, then the source may be obtained by the decomposition of the estimated seismogram autcccvariance ( Robinson, 1967 ). The minimum phase assumption may be true f o r explosive sources, but i t i s not generally true for 6 earthquake sources. Two of the three examples used in t h i s thesis ( West A u s t r a l i a and Kern County ) appear to have non-minimum phase sources. The second assumption mentioned above i s also i n v a l i d because, within a p a r t i c u l a r time window of in t e r e s t on teleseismic recordings, there are usually only a few a r r i v a l s . For the case when a p a r t i c u l a r a r r i v a l on the record i s separated from i t s neighbours by a time that i s greater than the source length, then an obvious method of source estimation i s simply to pick t h i s a r r i v a l as the source. Seismograms recorded near a delta cf 30 degrees often have t h i s feature (De y-Sarkar, 1974). However, one has tc be confident that there i s only one a r r i v a l i n the p a r t i c u l a r time window, and that the window contains a l l cf the source. Also, i f the source wavelet i s non-stationary, then t h i s estimate may not be the " b e s t - f i t t i n g " one. This method was used to estimate the scurce time function for the Kovaya Zemlya example. A simple extension cf this method i s to average a p a r t i c u l a r a r r i v a l over several traces after suitable time alignment. This method does net use the redundant source information contained i n secondary a r r i v a l s such as the method proposed for t h i s thesis does. The simple time domain averaging method works well when the scurce i s strongly stationary or when i t i s desirable to base the source estimate on only one p a r t i c u l a r a r r i v a l . Another method" of source estimation which uses 7 homcmorphic transforms has been proposed ( Ulrych, 1971; Sto f f a et a l , 1974). This teethed replaces the minimum phase source and white noise impulse response assumptions by another which also appears to be peer, at least f c r complex earthquake sources. This method w i l l be discussed in some d e t a i l i n the next section because i t forms the basis cf the proposed estimator for this thesis, which i s presented in the l a s t section of t h i s chapter. 8 2. Source Estimation by Homomor^hic Iransfor The concept of a homomorphic transform from a convolutional space to an additive space was f i r s t proposed by Oppenheim (1967), and was elaborated upon by Schafer (1969). Olrych (1971) has suggested that i t i s a s u i t a b l e method for seismic source estimation. The basic idea of the method i s to transform the convolution operator cf eguaticn [2-1] into an additive operator so that the system may be treated by li n e a r f i l t e r i n g theory. The transformation that accomplishes t h i s and i t s inverse are given below ( Oppenheim, 1967 ). x = H(x) = F-»[ lcg[ F( x ) ] ] [2-2] x = H- l( x ) = F-»{ exp[ F( x ) ] } [2-3] where H(«) denotes homomorphic transformation, F(«) denotes Fourier transform, and x i s the transform of x, and by d e f i n i t i o n i s termed the cepstrum. The phase information cf x i s preserved under the transform because the complex logarithm i s used 1. The homomorphic transform of equation [2-1], for the case when there i s no additive noise i s : x = H ( x ) = H(s<S>h) = H ( s ) + H ( h ) = s + h [2-4] lThe cepstrum defined here i s usually termed the complex cepstrum ( denoting the usage cf the complex logarithm ), to di s t i n g u i s h i t from the cepstrum defined by Bogert et al (1963) which uses only the rea l f a r t cf the logarithm. Since the "complex cepstrum" i s the only one used in t h i s t h e s i s , i t w i l l be referred to as simply the cepstrum. 9 Additive noise causes s p e c i a l problems under this transformation because there i s no simple expression for the horaomorphic transform cf the addition operator. Additive noise predominately a f f e c t s the phase part of the transform. This aspect w i l l be discussed l a t e r . SOURCE IMPULSE Figure 1. ____£__ of a Source and I repulse iii;J2£c_s_ Ce_£tru_ To separate s and h in equation [2-4] an assumption i s needed about their c e p s t r a l r e l a t i o n s h i p . The usual A assumption i s that the source cepstrum s has a lower A quefrency content than the impulse cepstrum h ( Schafer, 1969; Ulrych, 1971 ). Quefrency i s the independent variable of the cepstral domain ( i . e . i t i s the "frequency" content of the log spectra >. If t h i s assumption i s true, then the 10 s o u r c e may be a p p r o x i m a t e l y r e c o v e r e d by f i l t e r i n g c u t t h e h i g h e r o r d e r q u e f r e n c y t e r m s . T h i s a s s u m p t i o n i s i l l u s t r a t e d by t h e e x a m p l e o f F i g . 1. AMPLITUDE SPECTRUM f 1 1 r O. \- 2. 3. F i g u r e 2. A y e r a a a o f A l l West A u s t r a l i a JLffi£liJ-i.i-.~ S£S2i£s. The f e a t u r e marked by t h e a r r o w a p p a a r s on a l l o f t h e a m p l i t u d e s p e c t r a r e c o r d e d f o r t h i s e v e n t . However, t h e l o w q u e f r e n c y a s s u m p t i o n i s n e t c o n s i d e r e d t o be s u f f i c i e n t l y v a l i d f o r s o u r c e e s t i m a t i o n o f t e l e s e i m i c e v e n t s , I f t h e z - t r a n s f o r m o f t h e s o u r c e f u n c t i o n h a p p e n s t o c o n t a i n a z e r o n e a r t h e u n i t c i r c l e , t h e n t h e a m p l i t u d e and p h a s e s p e c t r a w i l l c h a n g e q u i t e r a p i d l y i n t h e v i c i n i t y o f t h i s z e r o ( C l a a r b c u t , 1974 ) . I n t h e c e p s t r a l domain t h i s t r a n s l a t e s a s h i g h e r o r d e r q u e f r e n c y t e r m s f o r t h e s o u r c e . The W i s t A u s t r a l i a e v e n t i s g i v e n as an i l l u s t r a t i o n o f t h i s e f f e c t . The a m p l i t u d e s p e c t r a c f a l l 11 the recordings of t h i s event that are presented i n t h i s thesis were averaged together to produce the amplitude spectrum shown in F i g . 2. The presence of the dip ( marked by an arrow ) indicates that there i s a zero of the z-transform near the unit c i r c l e . Since t h i s feature appears in a l l of the amplitude spectra i t i s assumed tc be a scurce e f f e c t and not part of the impulse response. This dip makes contributions to the source cepstrum at guefrencies cf 1.1 seconds and higher. Since the impulse cepstrum i s expected to have contributions i n t h i s range, the source cepstrum cannot be e f f e c t i v e l y recovered by just f i l t e r i n g out higher guefrencies. The amplitude spectra of the Kern County recordings exhibited a s i m i l i a r behavior. Several p r a c t i c a l problems exist in the computation cf homomorphic transforms. The main problem i s the unwrapping of the phase curve. The complex logarithm i s net urigue with respect to i t s imaginary component ( the phase ) so a fac t o r of ±2nTf, n=0,1,,..., must be added at each point tc make the phase curve continuous ( Schafer, 1969 ) * The usual procedure for unwrapping the phase curve i s tc compare the difference between the phase at a p a r t i c u l a r lag with the phase at the preceding lag. If t h i s difference i s greater than rr , then a ±2 IT factor i s added to the phase to make the difference l e s s than TT ( Schafer, 1969 ). This procedure is started at zero freguency where the true phase 12 i s known t o be z e r o 1 . The l i n e a r p hase c o m p o n e n t , w h i c h c o r r e s p o n d s t o a s h i f t i n t i m e o f t h e s i g n a l , i s r e m o v e d f r o m t h e u n w r a p p e d p h a s e c u r v e t o p r e v e n t i t f r o m d o m i n a t i n g t h e c e p s t r u m . A d d i t i v e n o i s e c a n make t h e u n w r a p p i n g p r o c e d u r e u n s t a b l e f o r two r e a s o n s . F i g u r e 3. E f f e c t o f A d d i t i v e H o i s e on A m p l i t u d e and P h a s e S ^ ^ c t r a The d a s h e d c i r c l e a p p r o x i m a t e s one s t a n d a r d d e v i a t i o n o f a d d i t i v e n o i s e . The c o r r e s p o n d i n g p h a s e and a m p l i t u d e d e v i a t i o n s a r e shown. The f i r s t r e a s o n i s t h a t t h e p h a s e d e v i a t i o n due t o a d d i t i v e n o i s e c a n become l a r g e when t h e s i g n a l a m p l i t u d e i s l o w , To i l l u s t r a t a t h i s , s u p p o s e t h a t a p a r t i c u l a r s a m p l e o f t h e s i g n a l s p e c t r u m has G a u s s i a n n o i s e a d d e d t o i t s r e a l i l f t h e t i m e s s r i e s h a s a n e g a t i v e D.C. c o m p o n e n t , t h e n i t i s i n v e r t e d b e f o r e t h e t r a n s f o r m i s t a k e n ( S c h a f f e r , 1969 ) . 13 and i m a g i n a r y p a r t s . F i g u r e 3 i l l u s t r a t e s two c a s e s . When t h e s i g n a l a m p l i t u d e i s h i g h , t h e phase d e v i a t i o n s r e m a i n a t a r e a s o n a b l e l e v e l . H o w e v e r , when t h e s i g n a l a m p l i t u d e d r o p s , t h e n t h e pha s e d e v i a t i o n becomes q u i t e l a r g e . I n t h e c a s e where t h e s i g n a l a m p l i t u d e d r o p s t c z e r o , t h e a m p l i t u d e d e v i a t i o n s assume a R a y l e i g h d i s t r i b u t i o n and t h e phase h a s a u n i f o r m d i s t r i b u t i o n e v e r [-Tr»7T] ( B r a c e w e l l , 1965 , p.5 ) . T h u s , f o r l o w a m p l i t u d e r e g i o n s o f t h e s i g n a l s p e c t r u m , t h e p h a s e c a n have a l a r g e v a r i a n c e . The s e c o n d p r o b l e m a l s o o c c u r s w i t h low a m p l i t u d e r e g i o n s o f t h e s p e c t r u m . When t h e a m p l i t u d e s p e c t r u m has a d i p i n i t , i t i n d i c a t e s t h e p r e s e n c e o f a z e r o i n t h e z-t r a n s f o r m n e a r t h e u n i t c i r c l e . S u ppose t h a t t h e t r u e s i g n a l h a s a z e r o j u s t i n s i d e t h e u n i t c i r c l e and t h a t a d d i t i v e n o i s e f o r c e s i t j u s t o u t s i d e t h e u n i t c i r c l e . The ph a s e o f t h e t r u e s i g n a l w o u l d i n c r e a s e by V n e a r t h e z e r o , b u t when t h e n o i s e i s ad d e d t h e phase d e c r e a s e s by rt a t t h i s p o i n t . F i g u r e 4 shows t h i s e f f e c t on t h e u n w r a p p e d p h a s e as w e l l as on t h e u n w r a p p e d phase w i t h t h e l i n e a r t r e n d r emoved. I t i s o b v i o u s t h a t a c h a n g e i n l o c a t i o n c f z e r o s n e a r t h e u n i t c i r c l e c a n d r a s t i c a l l y a l t e r t h e s h a p e of t h e unwrapped p h a s e c u r v e w i t h t h e l i n e a r c omponent remov e d . S e v e r a l m e t h o d s were t r i e d i n an a t t e m p t t c remove t h e p h a s e i n s t a b i l i t y a t lew a m p l i t u d e s . A l t h o u g h none o f them worked s u f f i c i e n t l y w e l l t o be deemed u s e f u l , t h e y a r e l i s t e d h e r e b e c a u s e ' t h e p r o b l e m i s u s u a l l y i g n o r e d i n t h e 14 RMPLITUDE SPECTRUM • LQCRTIQN OP ZEROS IT UNWRRPPED PHASE L.P.C. REMOVED FREQ.(RRD./SEC) EREQ.(RRD./SEC) Figure 4. E f f e c t of ________ a Zero Location on t h e _____ When addir.ive noise changes the location of zero near the unit c i r c l e the deviations i n the unwrapped phase can become large. L.P.C. denotes l i n e a r phase component. l i t e r a t u r e . 1) The phase was estimated by in t e r p o l a t i o n from i t s neighbouring points when the amplitude dropped below a certain l e v e l . Two problems make thi s method impractical, F i r s t , i t i s d i f f i c u l t to define a cutoff l e v e l which r e j e c t s a l l poorly estimated phase terms while ret a i n i n g a s u f f i c i e n t number of well estimated 15 ones to make i n t e r p o l a t i o n worthwhile. Secondly, the method tends to remove a l l Tf jumps, even the ones that should be present. 2) Stoffa et a l . (1974) suggested computing the derivative cf the phase d i r e c t l y from the real and imaginary parts of the Fourier transform and then integrating tc obtain the unwrapped phase curve. If the Fourier transform i s X (u)) + i Y (u>), then the phase derivative i s : ( Stcf f a et a l , 1974 ) = [ Y(u>) X'(LJ) - X (w) Y'(U>) ] / A M * [2-5] where A (u>) i s the amplitude spectrum. The phase derivative also becomes unstable at low amplitudes because of the A 2(^) term i n the denominator. Also, integrating the phase derivative in regions cf rapid phase changes can be a problem, because of poor sampling of the phase de r i v a t i v e . Interpolating the phase derivative over low amplitude frequencies removed a l l ft" jumps, as with the f i r s t method. Es s e n t i a l l y , the problems with the phase unwrapping techniques reduce to the following. While i t i s possible to unwrap the phase in regions of high spectral amplitude, there appears to be no way of connecting the unwrapped sections of the phase curve across regions of low spectral ampli tude. I concluded that the f i l t e r e d cepstrum was not a good source estimator whan i t was used on a single recording. 16 The phase estimate i s unstable in the presence cf ncise and there i s doubt of the general v a l i d i t y cf the low quefrency assumption. However, more sucessful r e s u l t s were obtained by averaging the log amplitude and phase spectra of a number of records from the same event. This tends to make the phase more stable and i t frees the method from the low quefrency assumption. This estimator w i l l be discussed in the next section. 17 3. Source Estimation b_y_ Averac]ed L ccj Amplitude and Phase Spectra A method of averaged log amplitude and phase spectra i s proposed as a technique for freeing the single trace homomorphic transform from the low quefrency assumption and for improving the phase s t a b i l i t y . Since the expectation operator i s commutative with the Fourier transform, averaging the log amplitude and phase spectra i s equivalent to averaging the cepstra. When the phase i s averaged, i t w i l l be the raw phase and not the unwrapped phase which i s used, so that the i n s t a b i l i t i e s of the unwrapping procedure can usually be avoided. The object of the averaging w i l l be to enhance the source part of the cepstrum while diminishing the impulse response part. If the averaging succeeds i n eliminating most of the impulse cepstrum, then the lew quefrency assumption need net be applied to recover the source. If the source i s stationary, and tha earth did not disperse or phase-shift, then the impulse response i s a sum of Dirac impulses. _N_ h t = )_ c * t _ T " ) C 2~ 63 where C n i s the amplitude cf the n a r r i v a l and T* i s i t s a r r i v a l time. Dispersion of the a r r i v a l s reduces the higher frequency components of the source. Thus, the Dirac impulses w i l l be broadened out so they resemble Gaussian pulses. Phase s h i f t i n g w i l l phase s h i f t the Dirac impulses as w i l l be 1 8 shown i n the l a s t section of the chapter on deccnvcluticn. These effects w i l l be considered tc be minor perturbations on the impulsive nature of equation [2-4], The Fcurier transform of t h i s equation i s : • H = ^ C n exp(-i<0To ) [2-7] The log amplitude ( leg A {u>) ) and phase ( t (to) ) of H (u>) are: l o g & M = log{ £ C,? + 2 ^_ cftCmccs<o( T n- T m) } [2-8] M r c* sinoTrt i 0 (w) = tan-i | <gj | [2-9] L V~ C^cosOT* -» The log amplitude w i l l have i t s main quefrency contributions at a l l combinations of differences of the various t r a v e l times. The average of the leg amplitude spectra over a su i t e cf impulse responses w i l l tend to a constant i f the combinations of t r a v e l time differences are random from one impulse response tc the next. In the time domain, t h i s means that the various a r r i v a l s must have dif f e r e n t phase v e l o c i t i e s ( mcvecuts ). The constant that the average of the log amplitude snectra w i l l tend tc i s the N average of the energy in each impulse response ( the c * 2 t\-i term of equation [2-8] ). This w i l l introduce a m u l t i p l i c a t i v e constant into the source estimate which i s not considered s i g n i f i c a n t . The phase spectra w i l l also have contributions which correspond tc the various t r a v e l times so that an average of the phase over a suite cf impulse responses w i l l tend to zero under the same 19 conditions as with the amplitude spectra. The assumption that t h i s source estimator makes i s that the a r r i v a l times cf the seismograms are s u f f i c i e n t l y random that the average cepstra s i g n i f i c a n t l y reduces the impulse response e f f e c t . In general the e f f e c t w i l l not average to zero, but i f i t i s reduced s u f f i c i e n t l y , the scurce wavelet w i l l be enhanced when the cepstrum i s transformed to the time domain. It i s usually necessary to apply a tine windcw to the source wavelet to remove any remaining e f f e c t s of the impulse response. If there is a secondary a r r i v a l on the seismograms that does not moveout s i g n i f i c a n t l y with respect to the primary s i g n a l , then t h i s method would have trouble separating the source from the combination of the two a r r i v a l s . This brings up the basic question of the resolving c a p a b i l i t i e s of a p a r t i c u l a r suite of seismograms. If a scurce i s estimated and i t appears to have two separate parts, then i t has to be decided whether the source estimate actually contains two separate a r r i v a l s or whether both parts belong to the source. The easiest method of deciding t h i s question i s to check for the presence cf a dip ( c f . Fig. 2 ) in the amplitude spectrum of each recording. This dip corresponds to the time separation cf the two apparent components. If t h i s feature i s present on a l l records and i t s frequency p o s i t i o n does not change, then i t may be concluded that, within the resolving power of the s u i t e of seismcgrarcs, both parts belong to t h e s c u r c e . Examining the amplitude spectra 20 in t h i s manner i s mere r e l i a b l e than locking d i r e c t l y at the records themselves for apparent moveout because; f i r s t , the feature i s usually more obvious i n the amplitude spectra, and secondly, the amplitude spectra uses the redundant information of a l l the a r r i v a l s on the record. To form the scurce estimate, the leg amplitude and phase spectra of a number of records are averaged together. I found that the best re s u l t s were obtained, as would be expected, when noisy records or records with obviously d i f f e r e n t sources were excluded. A s u f f i c i e n t number of them must be retained to make the averaging process e f f e c t i v e . Some care i s necessary when the phase i s averaged because of the a r b i t r a r y branch cut in the arctangent function at ±TT . This break would only be suitable i f the expected value of the source phase i s zero. However, i f the expected value i s net zero, then the branch should be placed Tr radians above and below t h i s expected value. An i n i t i a l estimate of the phase i s necessary to e s t a b l i s h the correct i n t e r v a l to average over. Each phase value i s then placed in this i n t e r v a l by adding the appropriate ±2TV i f i t s difference from the i n i t i a l estimate is greater than TT . The i n i t i a l phase estimate was taken as the phase value whose amplitude spectrum i s closest tc the average amplitude spectrum at each frequency sampling point. This i s the phase that i s most l i k e l y tc be independent of the phase spectrum of the impulse response. It i s necessary that 21 scale factor differences amongst the amplitude spectra be removed before the averaging in order that the i n i t i a l phase estimate i s not biased by energy differences in the seismograms. Rather than normalizing by the peak trace height or the peak amplitude spectrum (toth of which could be strongly affected by the the impulse response) a le a s t squares f i t t i n g of a m u l t i p l i c a t i v e constant between amplitude spectra was dene. The f i r s t trace was picked as a reference and the following least squares problem was sclved for the { Cj }. e 2 = ^ { A, (w) - Cj Aj (co) }2 [ 2- 10 ] where Aj [uj) i s the amplitude spectrum of the j r h recording. The summation is taken ever the frequency passband of the recordings. The solution of Cj which minimizies the error is : Cj = £3 A. ^ hA <w> / L fl'2(uJ) [2-11} ^ CO The amplitude spectra cf a l l other traces are then normalized by the factors. Aj M = Aj M / Cj j=2...N [2-12] After the averaging of the log amplitude and phase spectra has been dene, the lew quefrency source assumption may be applied, i f i t i s desired. To do t h i s , the f i n a l inverse Fourier transform of the hemomerphic transform i s taken and the r e s u l t i n g cepstrum i s zeroed above some quefrency l e v e l . However, as mentioned e a r l i e r , the low quefrency assumption i s a dubious one for complex earthguake 22 sources. The assumption, i f i t i s used at a l l , must therefore be applied with caution. Furthermore, since the impulse response i s diminished by the averaging process, there w i l l generally be no clear i n d i c a t i o n of the correct upper quefrency l i m i t . The source estimation method outlined above was derived on the basis of homcmcrphic transforms. This method may also be viewed as just an application of a p a r t i c u l a r expectation operator in the Fourier domain. In t h i s case, the expectation operator i s w E(2} = exP{ )_ lcg{Zj) } [2-13] If Z; i s written as Z •= A;exp(i0,-), then the expectation operator reduces tc E{Z) - exp{ )_ In A, + i __ 0; ] [2-14] J 3 1 J*' It i s possible tc form scurce estimators based on other expectation operators such as the L 2 norm (the standard average), and the L, ncrm (the median) estimators. The I_ norm estimator i s commutative with the Fourier transform, and hence i s equivalent tc time domain averaging. This form of source estimation was discussed in the i n t r c d u c t i c n tc Chapter 2. The L, norm estimator has the advantage that i t tends to r e j e c t " e r r a t i c " data values (Claerbout and Puir, 1973). However, i t was found that source amplitude spectra estimated with the L, ncrm expectation operator, has a higher quefrency content than the estimator based on equation [2-14]. This means that the I, norm estimator i s 23 not as e f f e c t i v e at removing impulse response e f f e c t s . Another possible modification to equation [2-14] i s to use standard averaging for the amplitude spectra. That i s E{Z} = E( Aj) exp{i £_ 0- } [2-15] Tests with numerical examples indicated that source estimates based on equations [2-14] and [2-15] have s i m i l i a r performances, but that the estimator cf [2-14] came closer to representing the true source amplitude shapes. 24 I I I . DECONVOLUTION OF THE SEISHCGBMS A method w i l l now be outlined for the deccnvcluticn cf the i n d i v i d u a l seismcgrams with a given source estimate. The desirable properties of a deconvolution methcd are: (1) i t resolve the a r r i v a l s times sharply, and estimate the a r r i v a l amplitudes accurately, (2) i t should be stable with respect to small scurce estimation errors or source non-s t a t i o n a r i t y , and (3) i t should be robust with respect tc the r e j e c t i o n of random ncise. The basic method that w i l l be used i s that cf spectral d i v i s i o n a l deconvolution, because i t does not require assumptions about the form of the source wavelet < i . e . that i t be minimum phase ) . Two modifications w i l l be made to improve the s t a b i l i t y and resolution of t h i s methcd. The f i r s t modification w i l l be the introduction of a minimum allowable source amplitude l e v e l (Helmberger and Wiggins, 1971; Dey-Sarkar, 1974), termed the waterlevel, to reduce spurious noise components and the ef f e c t of small errors in source estimation. The second modification i s the bandwidth extension of the impulse response beycnd i t s optimum passband, tc improve the time domain resolution of the deconvolution. The l a s t section cf thi s chapter deals with the envelope of deconvolution, which provides a useful t c c l f c r the detection of phase-shifted a r r i v a l s . 25 1. The Water l e v e l Parameter The frequency domain form cf the seismogram model ( equation [2-1] ) i s : X = S <» H + N [3-1] where the c a p i t a l s denote the Fourier transform pairs of the quantities i n equation [2-1]. To obtain the estimate cf the impulse response H , the spectrum cf the seismogram i s divided by the estimated source S . ~* _ X r S S , r S H = — = I — I H + I - ~ > N t 3 " 2 ^ S | S | 2-1 L I S I 2J where * denotes conjugation. As the estimated source amplitude becomes small ( JS'J—> o ), the factor multiplying H i s 0(1) l, assuming that S* does not deviate too far from S. However, the factor multiplying the noise component N i s 0 ( Vis'l )• Therefore, i t i s e s s e n t i a l tc establish a minimum amplitude l e v e l for the source to prevent the noise term from becoming tec large. This l e v e l can be thought cf as a l i m i t on the gain of the f i l t e r 1/s" i n the s p e c t r a l regions where the seismogram has l i t t l e or no information. The minimum source amplitude i s termed tha waterlevel and is conveniently expressed as a f r a c t i o n of tha maximum source amplitude. With t h i s modification the estimator of H becomes: »0(-») i s the "order-of" symbol, ( Si r o v i c h , 1 972 , p.5 ) 26 ^ S S H + S N max{ | S( 2, (k|S /! m o x ) 2} where k i s the waterlevel parameter ( 0<k< 1 ) and |S| m a ) ( i s the maximum source amplitude. The deconvolution w i l l be stable with respect tc small errors in the source estimation i f the factor • [3-4 ] S S max{|S|2, ( k | S U a A ) 2 ] goes to unity independently of frequency. To show the effe c t of the waterlevel parameter on this type cf error l e t the r e a l source te: S = "s + c R [3-5] where R i s an arbit r a r y function cf unit maximum amplitude and c i s a scale factor which indicates the degree cf error in estimation. The factor [3-4] new becomes: (S + CR)S |S| 2 + c R S [3-6] max{ | S| 2 r k 2 | ? | ^ M } max{ |S| 2,k2|s |2 a x } As with the noise, the factor multiplying cB i s C( 1/|S| ). Therefore, the introduction of the waterlevel alsc prevents t h i s factor from becoming tec large, i f c i s s u f f i c i e n t l y small. An optimum chcice for k would be one which makes kjSJ greater than both the noise l e v e l and the deviations cf the estimated source from the re a l scurce. Since these factors cannot ba determined in practice, an alternate procedure i s suggested. Rather than attempt tc pick an optimum k, I suggest that the deconvolution te performed for a range cf 27 k 6[0,1 ], The s t a b i l i t y cf the deccnvcluticn can be checked by comparing the impulse response for the various waterlevels. This was the procedure that was used for the three examples of t h i s thesis. The waterlevel parameter has another i n t e r e s t i n g i n t e r p r e t a t i o n . As k approaches zero, the estimator becomes an unrestricted deconvolution of X by S. As k approaches unity, the estimator i s just a scale factor times the cros s c o r r e l a t i o n of X and s'. The unrestricted deconvolution attempts to remove a l l cf the source effects from the seismogram, thus making i t the best estimator cf the true impulse response. This fcrm cf the estimator w i l l be the best f o r resolving the t r a v e l times. The crossc o r r e l a t i o n i s the least sguares estimate cf the a r r i v a l amplitudes ( Helmberger and Wiggins, 1971 ). If a pa r t i c u l a r phase arrives at time I, then the squared error due to amplitude differences between the a r r i v a l and the estimated source i s : e2 = Yi i x ( t _ T ) " A s (*) ) 2 [3-7] t The amplitude factor A which minimizes this e r r c r i s : A = ]T s (t) x ( t-T ) / £1 s (t) * [3-8] t t which i s just a scale factor times the cro s s c o r r e l a t i o n of x(t) and s ( t ) , evaluated at time T. If s (t) i s normalized to unit maximum amplitude, then A i s the a r r i v a l amplitude. To obtain A from the spectral d i v i s i o n a l deconvcluticn with k= 1, the peaks are multiplied by the maximum source spectral amplitude squared and divided by the energy-of the scurce. 28 The waterlevel can therefore be a l s c interpreted as a parameter which trades off amplitude r e s o l u t i o n with a r r i v a l time resolution,, This i s added incentive tc perform the deconvolution for a range of waterlevels. 29 2. I ^ i ^ ^ s i o n of the Impulse Response Spectrum The second modification that was introduced in t c the d i v i s i o n a l deconvcluticn methcd was an impulse response spectral extension scheme designed tc increase the time domain resolution of the estimated impulse response. This scheme e s s e n t i a l l y consists of assuming a model f c r the impulse spectrum and predicting the unknown spectral components on the basis of t h i s model. In the section on source estimation i t was shown that the assumption cf source s t a t i o n a r i t y led tc an impulse response which was a sum of Dirac impulses for a non-dispersive, non-phase s h i f t i n g earth ( equation [2-4] ). The Fourier transform cf t h i s r e l a t i o n i s : which i s just a sum of N complex sinusoids. When d i v i s i o n a l deconvolution i s performed, i t i s necessary to f i l t e r the impulse response to a passband defined by the band of s i g n i f i c a n t energy of the p a r t i c u l a r seismogram. If th i s were not done, then spurious noise information could be included in the impulse estimate as i s shown i n the example of Figure 6. This approach, however, makes a zero extension assumption about the nature cf the impulse response outside the passband and, in view cf the model of the impulse spectrum ( eguaticn [3-7] ), this does not seem to be a p a r t i c u l a r l y good assumption. A better assumption would be one which f i l l s in the impulse spectrum i n accordance with the continuous nature of t h i s model. } [3-9] 30 A procedure f c r extending short r e a l i z a t i o n s of a continuous process by the maximum entropy method ( Eurg, 1967 ), has been given by Smylie et a l . (1973) and Ulrych (1973). A prediction f i l t e r i s found by f i t t i n g the following autoregressive process to the known data: x i = H a K X t - K + z t [3-10] where {a K} are the c o e f f i c i e n t s of the prediction operator, p i s the order cf the autoregressive process ( the length of the prediction operator ), and z t i s a white noise process termed the innovation. An alternate way of writing the autoregressive model i s in terms of the prediction error cperatcr: __ X* X t - K = z t [ 3-11 ] where &c=-a^, k=1...p, and = 1. For Burg's maximum entropy algorithm the process i s f i t t e d by minimizing the variance cf the innovation when the prediction error operator i s applied in both the forward and reverse directions ( Smylie et a l . , (1973 ). The method w i l l always produce a minimum phase prediction error operator ( Burg, 1967; Claerbout, 1974 ) which i s i d e a l for extending the impulse spectrum because the energy w i l l net increase away from the knewn passband. Since the impulse spectrum i s complex, the complex form of Burg's algorithm w i l l be required tc compute the prediction operator. The derivation of the complex Burg algorithm follows alcng the same l i n e s as the r e a l version, except that i t must be realized that the prediction operator i s conjugated f cr 31 reverse prediction (Smylie et a l . , 1973). I t could ba argued that, since the impulse spectrum model i s a sum of sinusoids, i t wculd be better to f i t a sinusoidal model rather than an autoregressive model tc the known data. The problem with this approach i s that the sinusoidal model can only be assumed to be approximately v a l i d , due to dispersion and phase s h i f t s of the a r r i v a l s . The autoregressive model makes the more general assumption that the given impulse passband i s part of a continuous process which covers the entire spectrum. If the spectrum is t r u l y a sum of sinusoids, the autcregessive model does an excellent job, as long as there i s some additive noise present ( Ulrych, 1973 ). Figure 5. Prediction of a ______ Sinuscid The s o l i d l i n e i s the assumed passband and the dashed l i n e s are the predicted portions of the spectrum. Once the p r e d i c t i c n cperatcr has been determined by the maximum entropy method, i t i s applied in a unit step prediction fashicn_, u n t i l the enti r e spectrum i s f i l l e d i n . 32 Figure 5 shows the r e s u l t s for the simple case cf extending a single sinusoid with additive noise. The energy cf the predicted spectrum always decreases away from the optimal impulse passfcand because; f i r s t , the prediction error operator i s minumum phase, and secondly, the innovation of the autoregressive model i s net included in the predic t i o n . The one parameter that has to be chosen for t h i s scheme i s the length of the prediction operator. Ihe same problem occurs i n the estimation of power spectra by the maximum entropy method (' Ulrych and Bishop, 1975 ). Cne c r i t e r i o n that has been suggested f o r optimally choosing t h i s parameter i s that of the f i n a l prediction e r r c r (Akaike, 1969; Fryer,1974; Ulrych and Bishop, 1975 ). This c r i t e r i o n constructs a trade-off curve between the minimum cf the innovation variance and the variance of the estimate of the prediction c o e f f i c i e n t s ( Akaike,1969 ). The minimum cf t h i s curve i s taken as the optimal autoregressive order. However, observational experience indicates that for sinusoidal type processes, s i g n i f i c a n t resolution can be obtained by increasing the autoregressive order beyond that indicated by the f i n a l prediction error c r i t e r i o n . The penalty that one pays f c r doing sc i s the pcssible introduction of spurious noise components into the estimate. Instead of attempting tc find an optimum autcregressive order, I decided to use the same procedure as with the waterlevel parameter. That i s , the deccnvoluticn was 33 performed for a number of prediction lengths and the s t a b i l i t y of the deconvolution i s checked by comparing the various r e s u l t s . SOURCE AMPLITUDE SPECTRUM IMPULSE AMPLITUDE SPECTRUM 0^0 lTo 2JD 375 FREQUENCY (HZ) UNFILTERED IMPULSE RESPONSE FILTERED IMPULSE RESPONSE 0^ 10. 20. 30. 40 . TIME (SEC) Figure 6. Example of F i l t e r e d and Unfiltered This example i s the Kern County event recorded at ft YSD. The passband that was used tc f i l t e r the impulse response i s shown on the impulse spectrum. The impulse response passband for the c a l c u l a t i o n of 3 4 the prediction c o e f f i c i e n t s was determined fcy f i r s t performing a preliminary deconvolution in a passband determined from the main spectral energy of the scurce estimate. If any large noise spikes wera found to be present, the passband was reduced to exclude them. Figure 6 shows an example of these noise spikes. 35 3. Env§Icj£e_s_ of Dsconvolution Several authors have demonstrated that phase s h i f t i n g of a r r i v a l s can introduce systematic errors into t r a v e l time and amplitude observations by changing the r e l a t i v e peak and trough positions of the a r r i v a l wavelets ( Choy and Richards, 1975; H i l l , 1974; Helmberger and Viiggins, 1971 ). A simple procedure to overcome t h i s problem i s tc compute the envelope of the deconvolution ( Helmberger and Wiggins, 1971 ). The envelope of a signal x (t) , i s the modulus cf the a n a l y t i c s i g n a l which i s defined to be ( Bracewell, 1965, pp.267-272 ): x (t) = x (t) - i Hi ( x (t) ) [ 3- 12 ] where Hi(<») denotes the H i l b e r t transform. The a n a l y t i c s i g n a l i s most e a s i l y computed i n the frequency domain. X(UJ) = X (oJ) « ( 1 + sgn<~> ) [3-13] where the c a p i t a l s denote the Fourier tranform pairs and sgn i s the usual signum, function. The a n a l y t i c form of a phase s h i f t e d wavelet x ( t ) , i s a r o t a t i o n of i t s an a l y t i c zero phase shifted form, x ( t ) . x'(t) = exp ( i t ) x (t) [ 3- 14 ] where t i s the phase s h i f t angle. If x(t) i s crosscorrelated with or deconvolved by x{t), i t can be simply shown with the aid of equation [3-13] that the r e s u l t w i l l have the same phase s h i f t properties as x ( t ) . For example, the analytic deccnvoluticn d (t) ( the deconvolution of x (t) by x (t) ) i s related to the zero phase sh i f t e d 3 6 deconvolution 6. (t) ( the deconvolution of x (t) by x (t) ) ty: d (t) = exp ( ifc } d (t) [ 3-15 ] Hence, the envelope of deconvolution |d(t)j i s independent S Hi(S) PHRSE SHIFTS: 0° 45° 90° 135° 180° SYNTHETIC SEISMOGRAM DECONVOLUTION WITH ENVELOPE CRQSSCORRELRTIQN WITH ENVELOPE Figure 7. _ f the 'Envelopes of The sourca wavelet S, and i t s Fiilbert transform Hi (S) , that were used to construct the synthetic seismogram ara shown. The angles shown refer to the phase s h i f t s i n the three following traces. of the phase s h i f t angle i , making i t s peaks suitable for 37 timing purposes. The r e l a t i o n between the envelope and the deconvolution w i l l accurately r e f l e c t any phase s h i f t s r e l a t i v e to the estimated source wavelet that occur i n the o r i g i n a l signal. An example cf the use cf the envelope i s given in Fig.7. A useful way of displaying the envelope tc show the phase s h i f t properties cf the signal i s to superimpose i t on the deconvolution, as was done i n t h i s example. A 9C degree phase s h i f t , for example, i s characterized by a trough followed by a peak, of equal amplitude under the peak cf the envelope. The envelope was computed and displayed i n t h i s manner for the three deconvolution examples of t h i s t h e s i s . 38 IV. EXAMPLES OF DJCCJVOLUTIOJS The source estimation and deconvolution techniques outlined in the previous two sections were applied to three r e a l examples. These examples are: 1) The Kern County, C a l i f o r n i a earthquake of Sept. 16/1962, recorded on the L. R. S. N. stations in the delta range cf 10° to 39". 2) The PKP phase of the West Australia earthquake cf March 24/1970 recorded on the Canadian network in the delta range of 1140 to 150°. 3) The nuclear detonation at Ncvaya Zemlya on Oct. 14/1969 recorded cn the Canadian network i n the delta range of 20<> to 58°. No i n t e r p r e t a t i o n cf these examples was attempted; they were merely used to demonstrate tha enhancement properties cf the deconvolution methcd. The d e t a i l s of the source estimation and deccnvcluticn w i l l be given under separate subsections f o r each example. The format of tha deconvolution plots i s the same f c r a l l cases. The s t a t i o n name, the primary phase on the record, the name of the event, and the delta separation cf the s t a t i o n are given in the t i t l e f c r each recording. The f i r s t and second traces on each plot are the source estimation and the actual seismic recording respectively. The following traces are various deconvolutions cf the 39 recordings by the estimated source. The numbers tc the rig h t of these deccnvcluticns indicate the waterlevel ( the decimal number ) and the autoregressive order ( the integer number ) parameters that were used i n the p a r t i c u l a r daconvolution. Each of the seismic traces was delayed so that the f i r s t obvious a r r i v a l occurs at the ten second mark on the a r b i t r a r y time scale given with the plots. The impulse that corresponds tc th i s a r r i v a l should appear at the zero second mark, but to make a clearer presentation, the impulse response was delayed by ten seconds. 40 1. Kern Count__ 16__0_6 2 This example has an in t e r e s t i n g source wavelet which i n v a l i d a t e s the assumption of low quefrency sources. It i s apparent on a l l the records beyond 29° that the scurce has a precursor feature about f i v e seconds before the main event. This feature i s v i s i b l e on most records of the s u i t e . If the low quefrency assumption was applied to t h i s source, then the precursor would be taken as part of the impulse response, and not part cf the scurce. It i s also apparent by examining HTHN ( delta = 21.16° ) and HNME ( delta = 38.60° ) that some r e l a t i v e phase s h i f t i n g takes place between these two parts of the source, or that the earth's impulse response contains a phase s h i f t . A source was estimated for t h i s event by averaging the log amplitude and phase spectra of the stations marked by ast e r i s k s in Table 1. These recordings showed the clearest P phases. The scurce wavelet was windowed by a 20 second wi ndow. The enhancement achieved by deconvolution was judged to be good f o r this event. For every recording, the i n i t i a l P wave onset i s ea s i l y discernable, and f o r the stations MPAR, HTHN, and HMTN, the second a r r i v a l i s very well spiked. Two problems with t h i s deconvolution method are pointed out by t h i s example. The f i r s t occurs with the waterlevel parameter near unity, i . e . the deconvolution i s approaching a c r o s s c o r r e l a t i o n . Fcr t h i s case, both parts of the source I Station | Delta | Azimuth jScurce 1 | F i g . | 1 LCNM | 10. 120 ! I t 253.750 | | 8a 2 PMWY | 11.340 | 302.210 | | 8fc 3 HKWY | 12.220 | 302.870 | | 8c 4 HBOK | 15.640 | 273.330 | * | 8d 5 WNSD | 15.670 I 303977° | | 8e 6 AYSD | 16.510 | 3 03.60 0 | I 8f 7 CTCK | 18.810 | 272.760 | I eg 8 SEMN | 19.810 | 302.860 | | 8h 9 MPAR | 20.410 | 273.970 | I 8i 10 HTtftf | 21.16 0 | 302 . 990 | * I 8j 11 ARWS | 24. 60o | 302.900 | * | 8k 12 MMTN | 26. 320 | 279.150 | * I 81 13 GDVA | 28.910 | 284.020 | * | 8 EG 14 BLWV | 29.370 | 2 84.98 0 | * | 8n 15 BBPA | 32.690 | 292.240 | * | 8c 16 BUQB | 33.340 | 3 CO.68 0 | * 1 8 F 17 H N ME | 38.60o | 301.400 | * 1 eg j . . j _ J._ i * t indicates the stati o n was used in the source estimate Table 1. Kern County Seismogram Suite tend to be "spiked" because the au t c c c r r a l a t i c n of the source i t s e l f contains three spikes. Thus, i t i s important to perform the deconvolution f c r a number of waterlevels in order that t h i s problem can be detected. The second problem i s cver- r e s c l u t i c n by the impulse spectrum extension .scheme. In pa r t i c u l a r , the deccnvclution of WNSE ( delta = 15.670 ) suffers from this problem. It i s not cl e a r what i s causing the resonance in th i s recording. The deconvolution indicates the a r r i v a l of two i n t e r f e r i n g phases, but I believe t h i s r e s u l t should be treated with suspicion. KERN 15/9/S2 LCNM HKWT HBOK J CTCK MPRR HT.'-iN RRVS ft 10.12 11.34 12.22 15.64 15.67 16.51 18.81 19.81 21 .16 24.69 Figure 8. KERN 16/9/62 MMTN —JW^ l 26 35 28 91 Wa^ w^ V 29 37 l/^/V^VV- 32 69 BUQB 33 34 HNMS \ 38 60 . I 1 1 1 1 — 0. JO. 35. 33. O . 43 44 45 CTOK P KERN 16/9/62' DELTR=18.81 F i g u r e 8g. SEMN P KERN 16/9/62 DELTR=19.81 F i g u r e 8h. L*7 48 49 50 cn CD ID in in a o a —* —i —4 CD CO CO —< a o o CM • o • in • a - in • a in CO o o o —< o —« o a O —i o o O —i o —4 HNME P KERN 1C/9/62 DELTR=38.60 Figure 8q. Lrt 52 2. West Au s t r a l i a 24</3^ 70 The source was al s c estimated f c r this suite of records by the averaged log amplitude and phase spectra method. As with the previous example the source appears tc have two components to i t . The j u s t i f i c a t i o n f c r including both i n the source estimation l i e s i n the fact that a l l of the amplitude spectra cf the recordings have a d i s t i n c t i v e dip that does not change with frequency ( c. f. F i g . 2 ). This dip corresponds tc the time separation cf the two components. The stations that were used in the source estimate are marked by asterisks in Table 2. | Station | Delta | Azimuth (Source 1! F i g . | H F- "+ -+- . +-I 1 1 INK | 113.790 336.280 1 1 1 9a I 2 MBC | 116.720 | 345.750 1 | 9t I 3 ALE j 118.870 | 358.720 1 * | 9c I *» VIC | 119.260 | 314.080 1 * | 9d i 5 PST | 121.750 | 3 15. 0 5O 1 * | 9e I 6 YKC | 122.560 | 331.020 1 * | 9f | 7 RES | 122.950 | 347.810 1 * | 9g I 8 SES | 127.130 j 317.090 1 | 9h I 9 BLC | 129. 240 | 337.460 1 | 9i I 10 FFC | 131.270 | 324.46 0 I | 9j I 11 FCC | 133.170 | 332.030 1 | 9k I 12 FBC | 137.050 | 350.100 1 | 91 I 13 LHC | 141.090 321. 230 i | Sir | 14 GWC | 142.050 | 337.200 ! | 9n ! 15 SCH | 145.670 | 346.040 1 * | 9c I 16 OTT | 15C.30O | 326.870 1 * | 9p I »* i i n d icates the stat i o n was used in the source estimate. Table 2. West Australia Seismogram Suite 1 The deconvolutions f c r this event were performed only at the 0.95 and 0.10 waterlevels because i t was fcund that 53 tha source did not drop s i g n i f i c a n t l y below the 0.05 l e v e l . For each record the main PKP phase a r r i v a l i s e a s i l y distinguishable on the deconvolved records. The deconvolutions at RES, YKC, PUT, VIC, ALE, KBC, IKK, and CTT are e s p e c i a l l y good i n t h i s respect. In p a r t i c u l a r , the deconvcluticns cf MEC, flIE, INT, YKS, SES, and FCC appear to reveal the presence of the PKP r e f l e c t e d branch. The interpretation cf this feature w i l l be reported i n a further study. 54 INK W RUSTRRLIR 24/3/70 PKP DELTR=113.79 MBC W RUSTRRLIR 24/3/70 PKP QELTR=i15. 72 .X . 33. F i g u r e 9b, 56 57 58 59 FCC W RUSTRRLIR 24/3/70 PKP 0ELTR=133. 17 0.95 i V o io F i g u r e 9 k . FBC W RUSTRRLIR 24/3/70 PKP DELTR=137.05 F i g u r e 91 . LHC W RUSTRRLIfi 24/3/70 PKP DELTR=141 09 n i'A'l -A Art A fill i ll 0.95 0 0.95 0.95 15 < ' T O W 1 ^ | j ^ 1 5 . ^ ^ w O . 1 0 0.10 Figure 9m. GWC W RUSTRRLIfi 24/3/70 PKP DELTR=142 05 23. 30. Figure 9n. SCH W RUSTRRLIR 24/3/70 PKP DELTR=145 67 23. 30. Figure 9o. OTT W RUSTRRLIR 24/3/70 PKP DELTR=150 30 Figure 9p. 63 3. K o v a j a _ _ _ _ _ _ _ _ _ _ _ _ _ _ The Novaya Zeralya example presented d i f f i c u l t i e s in source estimation. Several authors have noted that t h i s s i t e i s prone to Rayleigh tc P wave conversions due to the large t o p o l o g i c a l r e l i e f ( Greenfield, 1971 ). Most cf the energy present i n the seismograms a f t e r the i n i t i a l I wave i s probably due to t h i s cause. Thus, i t was d i f f i c u l t tc assume source s t a t i c n a r i t y f c r th i s example. For th i s Station I — f -Delta I Azimuth |Source 1 -+ 1 Fig. 1 ALE \ 20.80° 17.85° 1 1 | 10a 2 KBC N | 30.42° | 2. 76° 1 | 1Gb 3 RES | 30. 870 j 15.09° 1 | 10c 4 FBC | 37.920 | 37.110 1 | 1Cd 5 INK | 3 8,32° ) 354.84 ° 1 | 10e 6 CUC N| 38.740 | 6. 10° 1 | 1 Cf 7 BLC | 41 .05° | 18.81° 1 I 10g 8 YKC | 44.06° | 7. 1 1° 1 | 1 Ch 9 FCC | 46. 300 | 21 .80° 1 | 10i 10 GWC 47.45° | 34.760 1 I 10j 1 1 STJ | 49.62° j 57.73° 1 | 10k 12 FFC | 5 1.03° J 17.09O 1 * | 101 13 FSJ | 52.30° | 0.69° I | 10m 14 SES | 55.96° | 10.85° 1 | 10n 15 PHC N I 56.03° | 35 8.29° 1 | 10o 16 OTT | 56.630 j 39.78° 1 I 1 Cp 17 PNT 57.37° | 4.32° I I 10q 18 VIC | 58.22° j 1. 390 1 | 1Cr i . j X 1 indicates the station was used in the source estimate. Table 3. _2IIOil Ze_l_a __i__o__a_ Suite reason, I decided tc estimate the scurce only on the basis of the d i r e c t P wave phase. The recording at FFC ( delta = 56.03° ) was used as the source estimate because i t appeared to have the c l e a r e s t d i r e c t P wave a r r i v a l , that was 64 unemcumbered by a converted secondary a r r i v a l . This example was included tc demonstrate the increased enhancement of the seismograms by t h i s deconvolution methcd over the multi-bandpass f i l t e r i n g technique used by Greenfield (1971). .Again, the 0.01 waterlevel was net used because the source amplitude did net drop below that l e v e l . The P wave onset i s well resolved on a l l the recordings as well as most of the secondary converted a r r i v a l s . the EES ( delta = 30.87° ) recording indicates the "spiking" r e s o l u t i o n power cf the methcd. Tha seismogram contains f i v e obvious a r r i v a l s which are a l l well resolved in the f i n a l deconvolution. 65 a CD co en o a, o CO n CM in co CD co in m in in in in in MBC N P NVZM 14/10/69 DELTR-30.42 • i & r filftfl „ A 0.95 S K A ' A 0 95 ^M^^^ -X] 'jiV |f| ^ J '9 l> J , , r , , j U. JO. Lfl. 30. <]. 30 Figure 10b. Figure 10c. F B C P N V Z M 1 4 / 1 0 / 6 9 D E L T R = 3 7 . 9 2 f r !l 0 . 9 5 A\ p ncr i l l 0 9 5 I 2 0 4!^^^ — ^ • 1 0 0 . 1 0 10 I ^ l j B > o ^ o » ^ W ^ » « ^ ^ » M ^ - ^ > | > ^ 0 . 1 0 I 2 0 Figure 10d. 68 r-co cn II ex LlJ a to > -z. C J 2Z CJ ID in CD 01 • o o o a —c in • CO O CM o o o • a • o O —i O CM 4-1 o a) CM cn CO cn II cn UJ CD CO 19 cu o a) 3 00 •H BLC P NVZM 14/10/69 DELTfl=41.05 i ; 1 : 1 1 a. i s , n . 3a. «o. so. Figure lOg. YKC NVZM 14/10/69 DELTR=44.G6 FCC P NVZM 1 4 / 1 0 / 6 9 D E L T R = 4 S . 3 0 C. I'J. 30. <£. ^Q. Figure 101. GWC P NVZM 1 4 / 1 0 / 6 9 D E L T R = 4 7 . 4 5 0 isVw- 0 9 5 f v v < 1 0 J 0 . 1 0 0 0 . 1 0 30 0 . 1 0 20 Figure 10j STJ P NVZM 34/:0/69 DELTfl=49.62 few I A vW>>-0.95 "0 0 .95 ' 10 0.95 ' 20 0.10 0 0.10 • 50 0.10 '20 F i g u r e 1 0 k . FFC NVZM 14/10/69 DELTR=5J.03 41 f -ib 0.95 ' 0 0.95 10 0.35 '20 0.10 ' 0 0.1G 30 i.10 20 F i g i r e 1 0 1 . FSJ P NVZM 14/10/69 DELTR=52.30 pi h n qc ^ f ™ # ^ ^ n ' Ij'l ' J I f K . 0.95 ~'flk'''j-i'<!^^ ] 0 iji l l U , - . 0.95 ' i f o.io ^:K r y^n^^ 20 Figure 10m. SES P NVZM 14/10/69 DELTfi=55.96 Figure lOn, D HC N P NVZM 14/10/69 DELTR=56.03 A 0.95 10 4f&!>*®h^ 2o S S T 1' /.A 0.10 '0 0.10 10 0.10 '20 Figure lOo. QTT p NVZM 14/10/69 DELTR=56.63 f /•If! i Mil Jf\Ty I i'h 1 0.S5 0 0.95 10 0.95 20 0.10 0 0.10 10 0.10 20 Figure 10p. l 1 ! : r — a. JO. 2 3 . 3 0 . 41 . Figure lOq. VIC P NVZM 14/10/63 DELTR=58.22 Figure lOr. 75 V. SUMMARY In t h i s thesis I have presented a method for the deconvolution of teleseismic recordings. A summary cf t h i s method and two related topics that ware discussed i s given i n the following paragraphs. The f i r s t step in the deconvolution i s the source estimation. Since earthquakes tend to produce ccmplex source time functions, the a b i l i t y tc separate the source and transmission path e f f e c t s on the basis of a single seismogram i s l i m i t e d . Fcr example, i f the earthquake has two or more successive motions i t i s d i f f i c u l t , i f net impossible, to decide whether the observed e f f e c t on a p a r t i c u l a r recording i s due to the source or the transmission path. The obvious remedy f c r this s i t u a t i o n i s to use the redundant source information available frcm a s u i t e of recordings of the same event, whose de l t a , azimuthal, and time ranges have been r e s t r i c t e d s u f f i c i e n t l y to allow the assumption cf source s t a t i c n a r i t y . The s p e c i f i c method of source estimation that i s suggested i s tc average the log amplitude and phase spectra of the recordings of the suite together. This method has the advantage over simple time domain averaging that i t uses the redundant information of secondary a r r i v a l s on the recordings. The necessary condition for this method to resolve the source time function i s that the various a r r i v a l s on the recordings have d i f f e r e n t phase v e l o c i t i e s . The actual deconvolution of the seismograms i s dene by 76 spectral d i v i s i o n a l deconvolution because i t dees net require the usual minimum phase assumption cf predictive deconvolution ( Bobinson, 1967 )„ Two modifications tc t h i s method were made, The f i r s t i s the introduction of a waterlevel parameter which constrains the minimum allowable source s p e c t r a l amplitude l e v e l ( Helmberger and Wiggins, 1971 ). The waterlevel parameter i s defined as a f r a c t i o n of the maximum scurce spectral amplitude, and i s useful for l i m i t i n g the gain of the deconvolution in regions where the seismograms contain l i t t l e or no information. The parameter also trades o f f a r r i v a l time resolution with a r r i v a l amplitude r e s o l u t i o n . The second modification i s the extension cf transmission path impulse response beyond i t s optimal spectral passband to increase the time domain resolution cf the deconvolution method. This i s accomplished by predicting the unknown spectral regions with a prediction operator determined by Burg's maximum entropy algorithm from the estimated impulse response spectrum in the c p t i n a l passband. The length of the prediction operator and the waterlevel mentioned above are two parameters which have to be chosen for t h i s method. Since in practice neither can be chosen optimally, I suggested that the deconvolution be performed for a range cf these parameter values. The s t a b i l i t y of the deconvolution can be checked by comparing the r e s u l t s for the various parameter values. Two related topics were also discussed in t h i s t h e s i s . 77 The f i r s t is source estimation by homomorphic transforms of single soismograms ( Ulrych, 197 1 ). This method was found to be an unsatifact cry source estimator for earthquakes because of the i n v a l i d i t y of the low quefrency assumption and the phase i n s t a b i l i t i e s of the transform i t s e l f . The averaged log amplitude and phase spectra source estimator mentioned e a r l i e r i s equivalent tc averaging the cepstra cf the seismograms together. It i s not necessary to unwrap the average phase curve thus removing seme cf the inherent phase i n s t a b i l i t i e s of the homomorphic transform. Furthermore, the low quefrency assumption i s replaced by the assumption that the averaging s u f f i c i e n t l y reduces the c e p s t f a l contributions of the impulse response. The second related topic of the thesis i s that cf the envelope of deconvolution ( Helmberger and Wiggins, 1971 ). The appearance of a deconvolved seismogram can be s u b s t a n t i a l l y improved by superimposing the envelope of deconvolution on i t . This envelope i s independent cf the a r r i v a l phase s h i f t s . Therefore, the relationship between the envelope and the deconvolution indicates the presence of any a r r i v a l s phase s h i f t e d r e l a t i v e tc the estimated source. 78 REFERENCES Akaike, H. ( 1 9 6 9 ) . F i t t i n g autoregressive models for prediction, Ann. Inst. Stat. Math., 2 J , 2 6 1 - 2 6 5 . Bogert, B.P., Healy, H.J.R., and Tukey, J.H., ( 1963). The quefrency alanysis of time s e r i e s fcr echces: cepstrum, pseudc-autcccvariance, cross-cepstrum, and saphe cracking, Proceedings of the Symposium 2& l i i l S Series Analysis, ed. M. Rosenblatt, Wiley and Sons, Inc., New York, 206-243. Bracewell, R. , (1965). The Fourier Transform and Its A£_2i2__i2 ni# McGraw-Hill, New York, p.5 and pp.267-2727 Burg, J.P., (1 967) . Maximum entropy s p e c t r a l a n a l y s i s : Presented at the 37th Annual Int e r n a t i c a l SEG Meeting, Oklahoma C i t y , Get. 31, 1967 Choy, G. L. and Richards, P.G., ( 1975). Pulse d i s t o r t i o n and Hilfcert transforms in multiply r e f l e c t e d and refracted body waves, Bu l l . Seism. Sec, Air. 65, 55-70. Claerbout, J. P. , and Muir, F., ( 1973). Robust modeling with e r r a t i c data, __c_h__ics, 38, 826-844. Claerbout, J. F. , (1974). Geo_hyical Data Pr cces s in _ , manuscript of a textbook to be published in 1975. Dey-Sarkar, S. K., (1974). Upper mantle F-wave vel o c i t y d i s t r i b u t i o n s beneath Western Canada, EhD thesis, Unversity of Toronto. 79 9. Fryer, G. J. , Odegard, M. E. and Sutton, G. H., (1974). Deconvolution and s p e c t r a l estimation using f i n a l prediction error. Submitted to Geophysics, 1974. 10. Greenfield, R.J., (1971). Short-period F-wave gener-ation by Rayleigh-wave s c a t t e r i n g at Ncvaya Zemlya, J . Gecr.hyj=. Res. 76, 7988-8002. 11. Helmberger D., and Wiggins, R.A., (1971). Upper mantle structure of the Midwestern United States, J„ ££°2his° Res. 76, 3229-3271. 12. H i l l , D.P., (1974) . Phase s h i f t and pulse d i s t o r t i o n i n tody waves due to i n t e r n a l caustics, B u l l . Seism. Soc. Am. 64, 1733-1742. 13. Oppenheim, A.V., (1967). Generalized superposition, i S J S I J a t i c n and Ccntrcl, 11, 5 28- 536. 14. Robinson, E. A,, (1967). S t a t i s t i c a l Communication and Hfzi.£.£t.i£2r Hafnar Publishing Co., Kew York. 15. Savage, J a C., (1966). Radiation from a r e a l i s t i c modal of f a u l t i n g , E u l l . Seism. Soc. Am. 56, 577-592. 16. 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 E l e c t r o n i c s , M.I.T., Tech. Beport 466. 17, S i r o v i c h , L. (1971), Technigues of Asvniptctic Analysis, Springer-Verlag, New-York, p.5. 18. Smylie D.E., Clarke G.K.C., and Ulrych, T.J. , {1973). Analysis cf i r r e g u l a r i t i e s in the earth's rotation, Methods i n Computational Physics, Academic Press7~iew~Ycrk, vcl.~13, 39 1-430. 19. Stoffa, P.L., Buhl, P., and Bryan, G. M. , H974), The application cf hcmcmcrphic deconvolution to shallow w?ter marine seismology Part I: Models, Geophysics, 39, 40 1-416. 80 20. Ulrych, T.J., (1971). Application of homcrccrphic decon-volution to seismology., _ e c _ h j s i c _ , J6, 65C- 660o 21. Ulrych, T.J., Smylie, D, E. , Jensen, C.G., and Clarke, G.K.C., (1973). Predictive f i l t e r i n g and smoothing of shcrt records by using maximum entropy, J . _eo_h_s. Res, 18, 4959-4964. 22. Ulrych, T.J., and Cishop, T.N., (1975). Kaxiirum entropy sp e c t r a l analysis and autoregressive decom-posi t i o n , Rev. _eo£h_s., 13, 183-200. 

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-0052880/manifest

Comment

Related Items