Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Application of digital filtering techniques for reducing and analyzing in-situ seismic time series Baziw, Erick John 1988

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

Item Metadata

Download

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

Full Text

APPLICATION OF DIGITAL FILTERING TECHNIQUES FOR REDUCING AND ANALYZING IN-SITU SEISMIC TIME SERIES by Erick J . Baziw B . A . S c , The University of B r i t i s h Columbia, 1986  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTERS OF APPLIED SCIENCE  in  THE FACULTY OF GRADUATE STUDIES Department of C i v i l Engineering  We accept t h i s t h e s i s as confirming to the required standard  Dr. R. G. Campanella, C i v i l Engineering Dept.  Dr. K. W h i t t a l l , Geophysics & Astronomy Dept.  THE UNIVERSITY OF BRITISH COLUMBIA August, 1988  ©  Erick J . Baziw, 1388  In  presenting  degree  this  thesis in partial fulfilment of the  for  an  of  department  this thesis for scholarly or  by  his  or  her  I further agree that permission for  purposes  permission.  Department The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3  extensive  may be granted by the head of  representatives.  It  is  understood  that  publication of this thesis for financial gain shall not be allowed without  nF.firt/ft-n  advanced  at the University of British Columbia, I agree that the Library shall make it  freely available for reference and study. copying  requirements  copying  my or  my written  ABSTRACT The introduction of d i g i t a l f i l t e r i n g i s a new and e x c i t i n g approach in analyzing i n - s i t u seismic data.  D i g i t a l f i l t e r s are also in the same s p i r i t  as the e l e c t r i c cone which replaced the mechanical cone in CPT* t e s t i n g .  That  i s , i t i s desirable to automate CPT t e s t i n g in order to make i t less operator dependent and increase the r e l i a b i l i t y and accuracy. In CPT seismic cone t e s t i n g seismic waves are generated at the surface and recorded downhole with v e l o c i t y or ac c el er at i on transducers.  The seismic  receivers record the d i f f e r e n t seismic wavelets ( e . g . , SV-waves, P-waves) allowing one to determine shear and compression wave v e l o c i t i e s .  In order to  d i s t i n g u i s h the d i f f e r e n t seismic events, an instrument with f a s t response time i s desired ( i . e . , high natural frequency and low damping). instrument i s c h a r a c t e r i s t i c of an accelerometer.  This type of  The f a s t response time  (small time constant) of an accelerometer r e s u l t s in a very s e n s i t i v e i n s t r u ment with corresponding noisy time domain c h a r a c t e r i s t i c s .  One way to  separate events i s to characterize the signal frequencies and remove unwanted frequencies.  D i g i t a l f i l t e r i n g i s ideal for t h i s a p p l i c a t i o n .  The techniques of d i g i t a l f i l t e r i n g introduced in t h i s research are based on frequency domain f i l t e r i n g , where Fast F o u r i e r , Butterworth F i l t e r , and c r o s s c o r r e l a t i o n algorithms are implemented.  One based on time domain  techniques, where a Kalman F i l t e r i s designed to model'the instrument and the physical environment.  The c r o s s c o r r e l a t i o n method allows one to focus on a  s p e c i f i c wavelet and use a l l the information of the wavelets present averaging out any noises or i r r e g u l a r i t i e s and r e l y i n g upon dominant responses.  The K a l -  rnan F i l t e r was applied in a manner in which i t modelled the sensors used and the physical environment of the body waves and noise generation.  The KF was  investigated for i t s possible a p p l i c a t i o n to obtaining accurate estimates on the P-wave and S-wave amplitudes and a r r i v a l times.  The KF i s a very f l e x i b l e  tool which allows one to model the problem considered accurately.  In a d d i t i o n ,  the KF works in the time domain which removes many of the l i m i t a t i o n s of the frequency domain techniques.  * Cone Penetration Test i i  The c r o s s c o r r e l a t i o n f i l t e r concepts are applied by a program referred to as CROSSCOR. CROSSCOR i s a graphics i n t e r a c t i v e program which displays the frequency spectrums, u n f i l t e r e d and f i l t e r e d time s e r i e s and c r o s s c o r r e l a t i o n s on a mainframe graphics terminal which has been adapted to run on the IBM P.C. CROSSCOR was tested for performance by analyzing s y n t h e t i c and real data.  The  r e s u l t s from the analysis on both synthetic and real data i n d i c a t e that CROSSCOR i s an accurate and user f r i e n d l y tool which greatly a s s i s t s one in obtaining seismic v e l o c i t i e s . The performance of the Kalman F i l t e r was analyzed by generating a source wavelet and passing i t through the second order instrumentation.  The second  order response i s then fed i n t o the KF with the a r r i v a l time and maximum amplitude being determined.  The f i l t e r was found to perform well and i t has much  promise in respect that i f i t i s f i n e l y turned, i t would be possible to obtain a r r i v a l times and amplitudes on l i n e r e s u l t i n g in v e l o c i t i e s and damping characteristics,  respectively.  iii  TABLE OF CONTENTS Page xv  ACKNOWLEDGEMENTS ABSTRACT  1i  LIST OF TABLES  vi  LIST OF FIGURES  viii  I.  INTRODUCTION  1  II.  DESCRIPTION OF PHYSICAL PROBLEM  8  III.  APPLICATION OF FREQUENCY DOMAIN TECHNIQUES A.  Discussion of Crosscorrelation F i l t e r Formulation  22  B.  F i e l d Program and Discussion of Results  48  B.l B.2  Research Sites Cross-Over Method Vs. CROSSCOR with Analogue F i l t e r Applied and Different Seismic Probes Used  51 61  CROSSCOR with F i l t e r e d (Analogue) and U n f i l t e r e d ( i . e . , Removing P-Waves and Noise with CROSSCOR) Sei smi c Traces .  72  B.4  CROSSCOR with Nonpolarized Sources  80  B.5  CROSSCOR i n Obtaining Compression Wave  B.3  IV.  22  Velocities  89  B.6  Bandwith Considerations when Applying CROSSCOR...  98  B.7  Conclusions  107  KALMAN FILTER FORMULATION FOR ESTIMATING AMPLITUDES AND ARRIVAL TIMES  109  A.  Kalman F i l t e r Results  B.  Description of Continuous Version of the Kalman F i l t e r Description of Discrete Version of the Kalman F i l t e r Derivation of Transfer Function f o r Geophones  C. D.  Ill  iv  123 127 131  TABLE OF CONTENTS (Cont'd) Page E. F. G.  KF Equations f o r Estimating Amplitudes and A r r i v a l Times  136  Comments on Computing Estimation Error Covariance M a t r i x , P  147  Evaluation of P a r t i a l s with Respect to x2 and x4  153  V.  SUMMARY AND CONCLUSIONS  159  VI.  SUGGESTIONS FOR FUTURE RESEARCH  163  REFERENCES  164  APPENDICES A.  CROSSCOR Program  B.  Discrete Model of Second Order System  178  C.  Program L i s t i n g of the KF Fomulation  181  ;  v  165  LIST OF TABLES No.  Title  B.2.1  Calculated V e l o c i t y P r o f i l e s from McDonald Farm Comparing the Cross-Over Method with CROSSCOR. The Data was Obtained with a S e l f - F i l t e r i n g (Natural Frequency = 28 Hz) Geophone  63  Calculated V e l o c i t y P r o f i l e s from McDonald Farm Comparing the Cross-Over Method with CROSSCOR. The Seismic Traces were Acquired with a 300 Hz Low-Pass F i l t e r Applied and Using the SCMP with P i e z o e l e c t r i c Bender Accelerometers  65  Calculated V e l o c i t y P r o f i l e s from Laing Bridge Comparing the Cross-Over Method with CROSSCOR. The Seismic Traces were Acquired with a 300 Hz Low-Pass F i l t e r Applied. In t h i s Set of Data i t was not Possible to Obtain Accurate CrossOver Picks  68  Calculated V e l o c i t y P r o f i l e s from Lower Langley (232) Comparing the Cross-Over Method with CROSSCOR. The Seismic Traces were Acquired with a 300 Hz Low-Pass F i l t e r Applied  70  Calculated V e l o c i t y P r o f i l e s from Lower Langley (232) Comparing F i l t e r e d (300 Hz Low-Pass) and U n f i l t e r e d Seismic Data. The Data was Acquired with an Accelerometer  75  Calculated V e l o c i t y P r o f i l e s from T i l b u r y S i t e Comparing Shear Wave V e l o c i t i e s Obtained with CROSSCOR from Seismic Data Acquired by Two Difference Accelerometers. In t h i s Set of Data i t was not Possible to Obtain Accurate Cross-Over Picks  78  Calculated V e l o c i t y P r o f i l e s of Lower Langley (232) Comparing Shear Wave V e l o c i t i e s Obtained from the Hammer Shear Source and the Buffalo Gun Source. The Data was Acquired with an Accelerometer ( U n f i l t e r e d )  83  Calculated V e l o c i t y P r o f i l e s for Annacis P i l e S i t e , where V e l o c i t i e s Obtained from S t r i k i n g each Side of the Hydraulic Pads of the Hammer Shear are Compared. The Data was Acquired with an Accelerometer ( U n f i l t e r e d )  : 85  B.2.2  B.2.3  B.2.4  B.3.1  B.3.2  B.4.1  B.4.2  Page  vi  LIST OF TABLES (Cont'd) No.  Title  Page  B.4.3  Calculated V e l o c i t y P r o f i l e s for Annacis P i l e S i t e , where V e l o c i t i e s Obtained from S t r i k i n g each Side of the Hydraulic Pads of the Hammer Shear are Compared. The Data was Acquired with an Accelerometer ( U n f i l t e r e d )  87  B.5.1 B.5.2  B.5.3  Congressional and Shear V e l o c i t i e s i n Sediments and Rocks Calculated V e l o c i t y P r o f i l e from Lower Langley, where an Accelerometer was Used for Data A c q u i s i t i o n . The Compression Wave V e l o c i t i e s were Obtained from the Buffalo Gun Source and by Applying a Bandpass of 200 to 1000 Hz Calculated V e l o c i t y P r o f i l e from Lower Langley where an Accelerometer was Used for Data A c q u i s i t i o n . The Compression Wave V e l o c i t i e s were Obtained from the Buffalo Gun Source and by Applying a Bandpass of 750 to 1500 H z . . . .  vi i  9 2  9 3  96  LIST OF FIGURES No.  Title  1.1  Schematic of an Acceleration Measuring Device  3  1.2  Unit Step Response Curve Showing Transient Response S p e c i f i c a t i o n s t d , t r , t p , Mp, and ts  4  1.3  Accelerometer and Geophone Frequency Response Curves. The Geophone I l l u s t r a t e d has a Natural Frequency of 28 Hz, and the Accelerometer has a Natural Frequency of 500 Hz  5  1.4  Block Diagram of Source, Earth, Recorder, F i l t e r and Velocity Calculation  7  2.1  Schematic of Seismic Transmission P r o f i l e  8  2.2  CPT Seismic Equipment Layout f o r Mechanical Source  10  2.3  Schematic of Recorded Data  11  2.4  Ratio of Reflected or Transmitted Energy to Incident Energy, Eout/Erin, with no Change i n Wave Type Occurring....  13  2.5  Ratio of Reflected or Transmitted Energy to Incident Energy, Eout/Erin, when Waves Change T y p e . . . . .  14  2.6  Phase Response of F i r s t Order Low-Pass and High-Pass Filters  16  2.7  E f f e c t of D i f f e r e n t F i l t e r C u t - o f f Frequencies on Accelerometer Responses  19  2.8  Seismic Data Acquired from Annacis Island VibroCompaction S i t e on May 18, 1988. Data Recorded with Accelerometer having a 300 Hz Low-Pass F i l t e r Applied. Diagram I l l u s t r a t e s D i f f i c u l t y in Obtaining Cross-Overs Due to Signal being Masked by Many Low Frequencies  20  Seismic Data Acquired from McDonald Farm on May 14, 1985. Data Acquired with Geophone with Natural Frequency of 28 Hz, and 100 Hz Low-Pass F i l t e r Applied. The Diagram I l l u s t r a t e s Ideal Signals for Determining Cross-Overs  21  I l l u s t r a t i o n of F i l t e r i n g Desired Frequencies  22  2.9  3.1  Page  vi i i  LIST OF FIGURES (Cont'd) No.  Title  Page  3.2  Sampling and A l a i s i n g D i f f e r e n t Frequencies Sampled at 4ms Intervals (250 Times Per Second), (a) 75 Hz S i g n a l ; (b) 175 Hz Signal Yields Same Sample Values as 75 Hz; (c) 250 Hz Signal Y i e l d s Samples of Constant Value (0 Hz)  3.3  Changes in Phase and Amplitude Spectra Resulting from Interaction of F i l t e r and Input  3.4  Example of a Box Car Function, G(w), and i t s Transform, s i n c ( . 5 a t )  25  3.5  Composition of Rectangular Wave from Harmonics  27  3.6  I l l u s t r a t i n g Gibbs' Phenomenon (a)  Frequency  Spectrum and (b) Truncated Signal  28  3.7  Amplitude Response of a Butterworth F i l t e r  29  3.8  Fourier Transform of a Signal with Nyquist wn  31  3.9 3.10  I l l u s t r a t i o n of the A l a i s i n g Problem when F i l t e r i n g Relationship Between the Deformed Angular Frequencies wd, f o r the B i l i n e a r Z Transform and the Actual Frequencies, w, of the D i g i t a l F i l t e r  31  3.11  32  An I n f i n i t e l y Long Sinusoidal and i t s ' . Fourier Transform  34  A Truncated 60 Cycle Signal and i t s ' Fourier Transform. The Signal was Switched on f o r #/2 Seconds  35  3.13  A Data Window Used to Taper a Discrete Function  36  3.14  Autocorrelation of a Process  37  3.15  C a l c u l a t i n g the Crosscorrelation of Two Functi ons  37  Padding Signals with Zeros at the End of t h e i r Time Series  38  Autocorrelation of Box Car Function  39  3.12  3.17 3.18  ix  LIST OF FIGURES (Cont'd) No.  Title  3.19  Signal with D.C. Offset of - 0 . 2 0  40  3.20  Signal with D.C. Offset of  40  3.21  Crosscorrelation of Offset Signals  41  3.22  Signal 1 Corrected for D.C. Offset with the Number of Points Changed from 200 to 150 because Dominant Response Occurs Between 40 to 60 Time Units  41  3.23  Signal 2 Corrected for D.C. O f f s e t . .  42  3.24  Crosscorrelation Value for D.C. Corrected Signals  42  3.25 3.26  Page  -0.06  Shift  Superimposed Signals to give Visual Aid f o r Estimating Time Lag  43  Processing Synthetic Data with CROSSCOR. (a) Input Trace (b) Frequency Spectrums (c) Bandpassed Traces (40 to 60 Hz) (d) C r o s s c o r r e l a t i o n Function of the Fi l t e r e d Traces  47  B.l  General S i t e Location Map  B.l.l  Cone Penetration P r o f i l e f o r McDonald Farm S i t e  52  B.1.2  Cone Penetration P r o f i l e f o r Grant McConachie Way S i t e  54  Cone Penetration P r o f i l e f o r T i l b u r y Natural Gas Plant S i t e  56  Cone Penetration P r o f i l e f o r Annacis P i l e Research S i t e  58  Cone Penetration P r o f i l e f o r Lower Langley 232 S t . S i t e  60  Seismic Data Acquired from Annacis Island V i b r o Compaction S i t e on May 18, 1988. Data Recorded with Accelerometer Having a 300 Hz Low-Pass F i l t e r Applied. Diagram I l l u s t r a t e s D i f f i c u l t y in Obtaining Cross-Overs due to Signal being Masked by Many Dominant Low Frequencies.  52  B.l.3 B.1.4 B.l.5 B.2.1  x  : 49  LIST OF FIGURES (Cont'd) No.  Title  B.2.2  (a) V e l o c i t y P r o f i l e from McDonald Farm Comparing the Cross-Over Method and CROSSCOR. The Data was Obtained with a S e l f F i l t e r i n g (Natural Frequency = 28 Hz) Geophone. (b) Cone Bearing P r o f i l e from McDonald Farm  64  (a) V e l o c i t y P r o f i l e from McDonald Farm Comparing the Cross-Over Method and CROSSCOR. The Seismic Traces were Acquired from an Accelerometer with a 300 Hz Analogue F i l t e r Applied, (b) Cone Bearing P r o f i l e from McDonald Farm  67  (a) V e l o c i t y P r o f i l e from Laing Bridge Comparing the Cross-Over Method and CROSSCOR v The Seismic Traces were Acquired with an Accelerometer with a 300 Hz Analogue F i l t e r A p p l i e d . In t h i s Set of Data i t was not Possible to Obtain Accurate CrossOver P i c k s , (b) Cone Bearing P r o f i l e from the Laing Bridge S i t e . . .  69  (a) V e l o c i t y P r o f i l e from Lower Langley Comparing the Cross-Over Method and CROSSCOR. The Seismic Traces were Acquired with an Accelerometer with a 300 Hz Analogue F i l t e r Applied, (b) Cone Bearing P r o f i l e from Lower Langley S i t e  71  (a) Seismic Trace Recorded at Grant McConachie. Way with the P - P l a t e Source ( U n f i l t e r e d ) . (b) Frequency Spectrum of the above Trace I l l u s t r a t i n g Three Responses Recorded by the Accelerometer: S-Wave at 100 Hz, P-Wave at 1000 Hz, and Resonanting Accelerometer at 300 Hz  73  (a) Seismic Traces from Grant McConachie Way from Depths of 3.8 and 6.8 Metres. These Wavelets are then F i l t e r e d with a Bandpass of 80 to 120 Hz in order to Extract the Shear Wave, (b) Corresponding C r o s s c o r r e l a t i on Function  74  (a) V e l o c i t y P r o f i l e from Lower Langley Comparing F i l t e r e d (300 Hz Low-Pass) and U n f i l t e r e d Seismic Data. The Data was Acquired with an Accelerometer. (b) Cone Bearing P r o f i l e from Lower Langley  76  B.2.3  B.2.4  B.2.5  B.3.1  B.3.2  B.3.3  Page  XT  LIST OF FIGURES (Cont'd) No.  Title  B.3.4  (a) Calculated V e l o c i t y P r o f i l e s from T i l b u r y S i t e Comparing Shear Wave V e l o c i t i e s Obtained with CROSSCOR from Seismic Data Acquired by Two D i f f e r e n t Accelerometers. The V e l o c i t y P r o f i l e s are also I l l u s t r a t e d with the Cone Bearing Profile. In t h i s Seismic Data i t was not Possible to Obtain Accurate Cross-Over P i c k s , (b) Cone Bearing P r o f i l e from T i l b u r y Island S i t e  79  Seismic Section from Lower Langley where a Buffalo Gun was Used as the Source and an U n f i l t e r e d Accelerometer Recorded the Data. From the Figure one can Clearly Distinguished the Compression and Shear Waves  81  Seismic Traces Form Lower Langley I l l u s t r a t i n g Triggering Problems with the Buffalo Gun Source. The Data was Acquired with an Accelerometer at a Depth of 13 Metres. The Resulting Time Difference was 2.72 Msec  82  Calculated V e l o c i t y P r o f i l e s from the Annacis P i l e Research S i t e , where V e l o c i t i e s Obtained from Both Sides of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer ( U n f i l t e r e d ) . (b) Cone Bearing P r o f i l e of the Annacis Bridge P i l e Research S i t e . . .  86  (a) Seismic Trace from Grant McConachie Way where a P - P l a t e Source was Used and an Accelerometer ( U n f i l t e r e d ) Recorded the Data at 6.8 Metres, (b) Corresponding Frequency Spectrum of Seismic Trace I l l u s t r a t i n g the Three Dominant Responses: Shear Wave at 100 Hz, Compression Wave at 1000 Hz, and Resonanting Accelerometer at 3000 H z . . .  90  V e l o c i t y Versus Density for Compressional and Shear Waves in a l l Types of Sediments and Rocks. Poisson's Ration Versus Density at Top of Figure  91  Seismic Traces Obtained at Lower Langley where a Buffalo Gun Source was Used. The Traces were F i l t e r e d (200 to 1000 Hz) i n Order to Extract the Compression Wave Responses. The Resulting Compression Wave Responses are Highly Correlated  95  B.4.1  B.4.2  B.4.3  B.5.1  B.5.2  B.5.3  Page  xi i  LIST OF FIGURES (Cont'd) No.  Title  B.5.4  Seismic Traces Obtained at Lower Langley where a Hammer Shear Source was Used. The Traces were F i l t e r e d (1750 to 2000 Hz) in Order to Extract the Compression Wave Responses. The Resulting Compression Waves are Poorly Correlated  97  Seismic Trace Recorded at a Lower Langley where a Hammer Shear Source was Used, (a) Seismic Trace at a 2.7 Metres Acquired with an U n f i l t e r e d Accelerometer (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response at 50 Hz  99  Seismic Trace Recorded at Lower Langley where a Hammer Shear Source was Used, (a) Seismic Trace at 3.7 Metres Acquired with an U n f i l t e r e d Accelerometer (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response at 50 Hz  100  (a) Seismic Traces Recorded at Depths of 2.7 and 3.7 Metres, (b) Trace in (a) F i l t e r e d Between 20 to 200 Hz Resulting in a V e l o c i t y of 103 m/sec. (c) Traces i n (a) F i l t e r e d Between 30 to 70 Hz Resulting in a V e l o c i t y of 103 m/sec.  101  Seismic Trace Recorded at Lower Langley where a Buffalo Gun Source was Used, (a) Seismic Trace at 3.7 Metres Acquired with an U n f i l t e r e d Accelerometer (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Responses Situated between 50 to 200 Hz.  102  Seismic Trace Recorded at Lower Langley where a Buffalo Gun Source was Used, (a) Seismic Trace at 4.7 Metres Acquired with an U n f i l t e r e d Accelerometer (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Responses Situated between 50 to 200 Hz  103  (a) Seismic Traces Recorded at Depths of 3.7 and 4.7 Metres, (b) Trace in (a) F i l t e r e d between 20 to 200 Hz Resulting in a V e l o c i t y of 94 m/sec. (c) Traces in (a) F i l t e r e d between 60 to 180 Hz Resulting in a V e l o c i t y of 94 m/sec  105  Trace in Figure B.6.6 (a) F i l t e r e d between 100 to 200 Hz Resulting in a V e l o c i t y of 94 m/sec  106  B.6.1  B.6.2  B.6.3  B.6.4  B.6.5  B.6.6  B.6.7  Page  xi i i  LIST OF FIGURES (Cont'd) No.  Title  Page  4.1  Block Diagram of System, Measurements, and Kalman Filter  110  4.2  4.3 4.4 4.5 4.6 4.7 4.8 4.9.a 4.9.b  Schematic Illustrating the Interaction of Input Noise, Body Wave and Recording Instrument With Resulting Output  n  ?  The Effect An Increasing Time Constant Has On White Noise  113  Low-Pass Noise With Corresponding Geophone Response as the Time Constant Increases  114  Geophone Response Due to P and S-Wave Impulses and Noise With Increasing Variance  117  Illustrating the Synthetic (a) Input and (b) Geophone Response  119  Updating the States (a) x(2) and (b) x(l) (Amplitude)  121  (Arrival Time)  Mean Square Error for (a) State x(2) and (b) State x(l)  122  Block Diagram of State & Measurement Equations (4.1) & (4.2)  126  Block Diagram of Continuous Kalman Filter Equations (4.6) & (4.7)  126  4.10  Schematic of An Acceleration Measuring Device  131  4.11  Schematic of Forces on Magnet and Voltages in Measurement Circuit Block Diagram of Second Order Geophone Transfer Function  132  4.12  xiv  136  ACKNOWLEDGEMENTS  The work o u t l i n e d i n t h i s t h e s i s would n o t have been p o s s i b l e without the funding and guidance o f Dr. Carapanella.  Dr. Ken W h i t t a l l from t h e  Geophysics Department was fundamental i n h e l p i n g me w i t h t h e g e o p h y s i c a l aspects of t h i s t h e s i s .  In a d d i t i o n , Drs. Robertson, Dunbar and Slauson  provided considerable advice i n preparing t h i s The  thesis.  r e s t o f t h e i n - s i t u group's h e l p was i n v a l u a b l e .  Don G i l l e s p i e and  John S u l l y were always a v a i l a b l e f o r a s s i s t a n c e i n f o r m u l a t i n g my i d e a s . I would a l s o l i k e t o thank Jim G r e i g and M i c h a e l E h l i n g f o r t h e i r with the i n t e r a c t i v e  help  graphics.  I would l i k e t o express my thanks t o my buddies a t t h e Savage Mansion, Damon, Steve and Bear f o r keeping my s p i r i t s up a t a l l times. L a s t l y , my parents were i n v a l u a b l e i n h e l p i n g me prepare my r e s e a r c h and t y p i n g i t up (Thanks, Mom!).  xv  I.  INTRODUCTION  The objective of t h i s research i s to apply d i g i t a l f i l t e r i n g techniques, in i n - s i t u seismic cone penetration t e s t i n g , inorder to determine S (Shear) and P (Pressure) wave v e l o c i t i e s .  There i s considerable i n t e r e s t in these  s o i l c h a r a c t e r i s t i c s because they provide i n s i g h t i n t o the response of s o i l to imposed loads such as b u i l d i n g s , heavy equipment, and dynamic loads from earthquakes and explosions.  Richart et a l . (1970), Mooeny (1974) and Borm  (1977) have w r i t t e n extensive l i t e r a t u r e on the r e l a t i o n of seismic wave v e l o c i t i e s to dynamic e l a s t i c moduli. The seismic cone penetration t e s t was f i r s t investigated at UBC by Rice (1984) where he found the a p p l i c a t i o n of the technique could provide a rapid and accurate method f o r carrying out a downhole shear wave v e l o c i t y survey. In his work, Rice determined that a u n i a x i a l geophone was s u f f i c i e n t for i d e n t i f y i n g shear waves, i f the instrument was orientated properly.  Rice also  investigated several d i f f e r e n t types of sources and the dominant waves they generated.  Rice worked predominantly with geophones which have slower  response time compared to an accelerometer.  From his work with geophones,  which were not amplified at the p i c k - u p , Rice stated "Below 30 metres the strength of i n d i v i d u a l energy impulses was attenuated such, that c l e a r i n d i v i dual shear wave excursions were not always i d e n t i f i a b l e . " The next to investigate the a p p l i c a t i o n of the seismic cone in i n - s i t u t e s t i n g was Laing (1985).  Laing focused her a t t e n t i o n on sources and receivers  with the seismic cone f o r both on-shore and o f f - s h o r e a p p l i c a t i o n s .  Laing came  up with the same conclusion as Rice with respect to the use of the geophone as a seismic r e c e i v e r , she s t a t e s , "accelerometers tend to show a response which more c l o s e l y represents the response of the s o i l .  Thus, damping c h a r a c t e r i s t i c s  of s o i l should t h e o r e t i c a l l y be a t t a i n a b l e from accelerometer responses." I d e a l l y one would l i k e an instrument to t e l l the whole seismic s t o r y . Figure 1.1 i l l u s t r a t e s a t y p i c a l a c c e l e r a t i o n and/or v e l o c i t y measuring device. The dynamics behind the V e l o c i t y and ac c el er at i on measurements can be expressed by the following two second order equations:  1  acceleration: X ( )/s 0  S  2  X i  = X (s)/X 0  = k/((s /a)^)+2^s/u) +l) 2  i  n  (1.1)  velocity: x (s)/sx Q  i  = x ( s ) / x . = s/(s +2^ s+a) ) 2  Q  (1.2)  2  n  r e l a t i v e displacement of mass M x0 x. = absolute displacement of mass M u  n  =.natural frequency  X = damping k = proportionality  constant  s = Laplace transform v a r i a b l e The response of the instrument i s described by w and c n  the higher the value of x, and smaller value of w in the s e n s i t i v i t y of the instrument.  n  In g e n e r a l ,  the greater the reduction  This corresponds to an instrument with  a large time constant r e s u l t i n g in slower t r a n s i e n t response s p e c i f i c a t i o n s , t h i s i s i l l u s t r a t e d in Figure 1.2.  In t h i s f i g u r e , an instrument with a  large time constant associated with i t would have comparatively longer response time s p e c i f i c a t i o n s .  Figure 1.3 shows t y p i c a l frequency response curves f o r  geophones and accelerometers.  As i s i l l u s t r a t e d , an accelerometer has f l a t t e r  response curves than a geophone f o r lower frequency i n p u t s , and becomes unstable f o r higher frequency inputs near i t s natural frequency.  Since P waves t y p i c a l l y  consists of frequencies between 200 to 2000 Hz, and S waves are t y p i c a l l y between 50 to 300 Hz, i t i s desirable to have an instrument which has a f l a t r e s ponse curve below 2000 Hz and i s l i g h t in mass, due to lack of room in the seismic cone.  The above conditions would be met with a l i g h t weight a c c e l e r o -  meter with natural frequency at approximately 4000 Hz. The desire of using s e n s i t i v e accelerometers in i n - s i t u t e s t i n g c a l l s f o r the necessity of using d i g i t a l f i l t e r s in order to separate out noise and d i f ferent seismic events.  The technique of data reduction used by Rice and Laing  r e l i e d upon the p o l a r i t y of shear body waves and the a p p l i c a t i o n of analogue f i l ters.  In Section I I ,  t h i s technique i s o u t l i n e d with some of i t s l i m i t a t i o n s .  2  Recorder  Figure 1.-1. . Schematic of an Acceleration Measuring Device (Close & Fredrick, 1978)" In Figure 1 . 1 , the following parametres are defined: M  -  mass of suspended magnet  B  -  flux  K  -  stiffness of spring which supports the magnet  D  -  viscous damping between magnet and the case  N  -  number of turns in coil  d  -  diameter of coil  l  -  length of c o i l , ir-d-N  R  -  resistance of coil  -  inductance of coil  -  resistance of recorder  -  acceleration of case relative to a fixed reference  -  equilibrium position of M relative to the case  -  dynamic displacement of M relative to the case  a (t) c  density between magnets  3  TRANSIENT RESPONSE SPECIFICATIONS  1. Delay time t 2. Rise time t  d  r  3. Peak time t  p  4. Maximum overshoot M 5. Settling time t,  4  p  a>  u o •a: 0) c o o  d •rt C  e Q  8 10  025  6  8  15 2 ? 25-3040 5060 80 100  100  150200  300  403fj  200  300  400  Frequency,  Figure 1.3.  600  500  1000  1000  2000  3000 4000  2000 3000  Hz  Accelerometer and Geophone Frequency Response Curves. The Geophone I l l u s t r a t e d has a Natural Frequency of 28 Hz, and the Accelerometer has a Natural Frequency.of 500 Hz.  4000  The dynamic behind source wavelet generation and f i n a l v e l o c i t y mination i s as f o l l o w s .  deter-  A source, such as a hammer blow or seismic cap,  generates p o t e n t i a l s at the surface which correspond to SV*, SH**, and P*** waves.  These i n i t i a l p o t e n t i a l s can be considered to be Dirac d e l t a functions  ( i . e . , impulses) having a broad frequency spectrum.  As these p o t e n t i a l s propa-  gate through the earth as seismic waves, t h e i r frequency spectra are bandlimited ( i . e . , higher frequencies are attenuated). a low-pass analogue f i l t e r .  This i s due to the earth acting as  At a c e r t a i n depth, the generated s i g n a l s are  recorded with the three waves ( i . e . , SV, SH, and P) and noise being superimposed. The recorded signals are d i f f i c u l t to analyze due to the superposition of many s i g n a l s .  Therefore, a d i g i t a l f i l t e r is. applied so one can focus on a  s p e c i f i c wavelet ( e . g . , SH).  Once f i l t e r i n g i s a p p l i e d , the time o f f s e t between  two s i m i l a r waves ( e . g . , SH) recorded at successive depths i s determined and interval  v e l o c i t i e s are c a l c u l a t e d .  The following diagram (Figure 1.4)  illus-  t r a t e s the e s s e n t i a l r e l a t i o n between source, e a r t h , recorder, f i l t e r , and interval velocity calculation. Referring to Figure 1.4, the major goal These f i l t e r s should  give  was to develop d i g i t a l  filters.  c l e a r signals without d i s t o r t i n g the o r i g i n a l  signal in terms of time s h i f t s and amplitudes. The two f i l t e r i n g methods considered in t h i s thesis are the one based on the frequency domain f i l t e r i n g (Section I I I ) ,  where Fast F o u r i e r ,  Butterworth,  and c r o s s c o r r e l a t i o n algorithms are implemented, and one based on the time domain techniques (Section IV), where a Kalman F i l t e r i s designed to model, both the instrument and the physical environment.  *  A h o r i z o n t a l l y t r a v e l i n g shear wave so polarized, that the p a r t i c l e motion is a l l v e r t i c a l .  ** A h o r i z o n t a l l y t r a v e l i n g shear wave so p o l a r i z e d that the p a r t i c l e motion i s all horizontal. *** Compression wave c o n s i s t i n g of p a r t i c l e motion with a l t e r n a t i n g condensations and rarefactions during which adjacent p a r t i c l e s of the s o l i d are c l o s e r together and f a r t h e r apart during successive half c y c l e s . The motion of the p a r t i c l e s are always in the d i r e c t i o n of wave propagation. 6  Depth 2  Depth 1  SV Source Potentials  SH  5^1  I  I  Noise P-Wave S-Wave  Noise P-Wave S-Wave  Earth Model  I  I Recording Equipment  1 Band-Pass Filter  <5> Superimposed Signal V  AV =6d  Figure 1.4.  Block Diagram of S o u r j ^ . J i a r t h , and V e l o c i t y C a l c u l a t i o n . 7  Velocity Calculations  Recorder,  Filter  II.  DESCRIPTION OF PHYSICAL PROBLEM  A seismic cone investigation is conducted in order to obtain i n - s i t u values of shear wave v e l o c i t i e s of a stratigraphic s e c t i o n .  The shear wave v e l o c i t i e s  are d i r e c t l y related to the maximum shear modulus (Rice 1984) of the s o i l by the following mathematical expression. f i  G_...  =  ^  -  max  •  2 V 's  where  V G  s  max  i s t h e mass d e n s i t y  o f the s o i l  - i s the shear wave v e l o c i t y J  - i s t h e maximum s h e a r m o d u l u s  10  at strains  less  than  percent.  The significance of the above s t r a i n level i s that at strains of 10"  per-  cent and l e s s , the shear modulus i s maximum and l i t t l e effected by s t r a i n amplitude.  The maximum shear parameter i s an important s o i l parameter; shear  modulus reduction curves can be generated from i t . The shear modulus reduction curves are very important in dynamic s o i l analysis,because they relate the change in the shear modulus with induced s t r a i n s . In seismic cone testing a response signal i s recorded at one metre i n t e r v a l depths.  The wavelet i s generated at the surface by a sledge hammer ^low or a  seismic cap and received at the cone with either a geophone or an accelerometer. The following diagram i l l u s t r a t e s the seismic cone transmission p r o f i l e . sledge hammer bl ow  ! j  <*iv  s  i  cone  1  •  — i  d  _ i  , 2 d  V  sl  s2  »S3  3  d'  -  d"  -  T  l  "  T  2  -  recording depth to < recording  *  1 i  inserted here  ! ' '  1 —rr -  d  1  "  d  "  *  T  arrive to f i r s t seismic recei ver time f o r shear to a r r i v e to second seismic receiver  ... Figure 2 . 1 .  Schematic of Seismic Transmission P r o f i l e 8  The geophone and accelerometer responses were recorded on a N i c o l e t 4094 d i g i t a l o s c i l l o s c o p e with a CRT screen and a floppy disk storage c a p a b i l i t y . A complete d e s c r i p t i o n of t h i s o s c i l l o s c o p e i s given by Rice (1984). The seismic cone i s pushed into the s o i l using the U.B.C. ing v e h i c l e .  I n - s i t u Test-  The truck i l l u s t r a t e d in Figure 2.2 weighs approximately 11 tons  and has two hydraulic c y l i n d e r s with which the cone i s pushed into the s o i l . To provide the reaction force needed to push the cone, the truck i s raised on two hydraulic pads as shown in Figure 2 . 2 . Wave sources.  The hydraulic pads also act as SH-  SH-Wave sources are desired because i t i s easy to reverse the  p o l a r i t y of the SH-Wave.  This energy reversal r e s u l t s in two oppositely  p o l a r i z e d wavelets which are recorded on the N i c o l e t o s c i l l o s c o p e as i s shown in Figure 2 . 3 .  The f i r s t major crossover times i l l u s t r a t e d in Figure 2.3  allow the c a l c u l a t i o n of the r e l a t i v e t r a v e l time d i f f e r e n c e of shear waves between successive depths which allows the c a l c u l a t i o n of shear wave v e l o c i t i e s . The above technique i s thoroughly discussed by Rice (1984) and Laing (1985), and w i l l be referred to as the reverse p o l a r i t y method and/or crossover method. The crossover times discussed above are r e l a t e d to the v e l o c i t y p r o f i l e by the f o l l o w i n g equations:  V  -d' d II corr corr  s3  (2.2a)  where d ii corr  (2.2b)  d corr  (2.2c)  and d , d " , T, , and T  9  are defined i n F i g u r e . 2 . 1 .  9  F i g u r e 2.2  CPT S e i s m i c Equipment Layout f o r M e c h a n i c a l Source, ( a f t e r R i c e , 1984)  10  one-way travel time (msec)  f i r s t arrivals  Figure 2 . 3 .  Schematic of Recorded Data  11  These equations and Figure 2.1 apply only i f i t i s assumed that the angle of incidence approaches zero, and the r a t i o of the i n t e r v a l v e l o c i t i e s i s of the order of u n i t y .  shear  These conditions w i l l generally hold  for i n - s i t u t e s t i n g , because v e l o c i t i e s f o r ajacently s t r a t i f i e d s o i l s do not d i f f e r (Reference 4) s i g n i f i c a n t l y ( i . e . , generally less than 1 . 5 ) , and the source i s not separated much from the r e c e i v e r , such that we have near normal incidence. Referring to Figures 2.4 and 2 . 5 , i t i s seen that for these conditions an incident SV-Wave w i l l have-a n e g l i g i b l e amount of energy converted to a P-Wave as r e f l e c t i o n or transmission (Figures 2.5b and d ) .  The majority of the i n c i -  dent SV energy w i l l be converted to an SV transmission with the angle of i n c i dence approximately equal to the angle of r e f r a c t i o n (Figure 2.4b). can be said for the P-Wave.  The same  An incident P-Wave subjected to these boundary  conditions w i l l transmit most of i t s energy, at an i n t e r f a c e , to a P-Wave in the lower medium with the angle of incidence equal to the angle of r e f r a c t i o n . Thus, the assumption of s t r a i g h t ray propagation, with l i t t l e r e f l e c t i o n and energy conversion, holds.  The above considerations are more throughly assessed  in Rice (1984), where he s t a t e s , "For most p r a c t i c a l purposes a s t r a i g h t  line  t r a v e l path assumption would be quite s a t i s f a c t o r y . " As was previously discussed, the i n t e r v a l times, T^ and J^,  are determined  by taking the d i f f e r e n c e of the crossover times f o r the successive wavelets. When only noise free shear waves e x i s t , the above procedure i s trouble f r e e . In many cases, one i s not dealing with an ideal wavelet, but has a seismogram with corrupting P-Waves and noise.  These p a r t i c u l a r signals make i t  to pick a r r i v a l and crossover times.  difficult  The method which i s currently used to  obtain c l e a r shear wave signals applys an analogue low-pass f i l t e r to remove : the higher frequencies of the recorded s i g n a l .  Stokbe and Hoar, (1978) raised  an important point when considering analogue and/or e l e c t r i c a l f i l t e r s .  They  suggest that f i l t e r i n g should be minimized because i t can s i g n i f i c a n t l y d i s t o r t the signal and erroneously a l t e r a r r i v a l times.  Stokoe and Hoar also  state that e l e c t r o n i c f i l t e r s usually introduce time delays which vary with the frequency of the input s i g n a l ; hence, u n f i l t e r e d signals were generally preferred and were used in t h e i r study.  A numerical example of the possible  phase s h i f t or corresponding time delay, f o r an analogue f i l t e r ,  12  follows.  a  7  Pi  l  1.036 1.07 1.07 1.0'J 1.10 1.09 1.00 1.14 1.30 1.103 1.286 1.2tl6 0 « t a i l i doubtful  2 3 4  ' in E  •7T  a  • {*>) S K .  ijji  •  i .1  <*l  °2  0.24 0.25 0.28 0.25  0.24 0.26 0.16 0.25  p - density of s o i l a - shear wave v e l o c i t y B - compression wave v e l o c i t y a - attenuation  II ' ]i /I' s  •  ^  <H  ct/  v  -  s  20°  \  </  Ii  i  L  40°  J  '  60°  —  i  —  80'  Angle of Incidence  Figure 2.4. Ratio of Reflected or Transmitted Energy to Incident Energy, " V Change in Wave Type Occurring (Ewing, J a r d e t z , and Press, 1957)  E  out  / -j > E  n  M  i  t  h  N o  Figure 2. 5.  Ratio of Reflected or Transmitted Energy to Incident E n e r g y , V E Waves Change Type (Ewihg, Jardetz, and Press, 1957)  QUt  / E , When in  The accelerometers used i n t h i s i n v e s t i g a t i o n were bimorph transducers, which were 0.5 inches square i n dimension and made of G1195 piezoceramic bender m a t e r i a l .  When mounted, these receivers had a natural frequency of  3 kHz, and underdamped response of 1.3 of c r i t i c a l .  More information of t h i s  accelerometer can be obtained i n technical papers from Piezo E l e c t r i c Product's, I n c . , New Jersey.  Figure 2.6 i l l u s t r a t e s the phase response curve f o r a  typical first-order N i c o l e t recorder.  low-pass analogue f i l t e r s i m i l a r to f i l t e r s w i t h i n the Equation (2.3a) i l l u s t r a t e s the i n t e r a c t i o n of the input  wavelet, f i l t e r , and output wavelet. |AU)|e  1  |Y(«o)|e  input  '  |A( )||Y( )|e  filter  M  j(* +<j> ) 1  2  u  (2.3a)  output  define A(UJ) = |A(w)|e  Y(co) = |Y(o.)|e  1  J4><  G(u>) = |A(to) | | Y(ui) | e  Input  (2.3b)  Filter  (2.3c)  Output  (2.3d)  angular frequency If i t i s assumed that the f i l t e r has unity g a i n , then Equation (2.3d)  reduces  to, G(u>) = |A(u)|e  j(+ +<|) ) 1  2  (2.3e)  The Fourier transform, f o r an aperiodic wavelet, of Equation (2.3'e) i s CO  g(t)  i r = ~- J  Jwt G(o))e dto ....  In terms of the frequency,•f, we have  15  (2.3f)  FREQUENCY  ,.25 .3  -90  .6 .7  .8.91.0  2.0  3.0  4.0  u  (a)  Low-Pass Response  K - Time constant S - Laplace transform v a r i a b l e  FREQUENCY 25  .3  .4  .5  .6 .7 . 8 . 9  1.0  2.0  3.0  4.0  0°  • out - in  3  •  KS—— -> • *  -45°-  <  X a.  -90°-  (b)  Figure 2.6.  High-Pass Response  Phase Response of f i r s t Order Low-Pass and High-Pass F i l t e r s (Lancaster, 1982) 16  g(t)  r = /  G(f)e  j2irft  df  (2.3g)  S u b s t i t u t i n g Equation (2.3e) into Equation ( 2 . 3 g ) , gives  g(t)  /  |A(f)|e  j(27Tft+Y(f))  df  (2.4a)  -co  where <l' (f)+'t' (f)  Y(f)  1  2  For an actual waveform, g(t)  i s real and hence Equation (2.4a) reduces to  (2.4b) -co  If in Equation (2.3a) i t i s assumed that <j>^=0 ( i . e . , there i s no phase lag in the input s i g n a l ) , then y(f)  (2.4c)  = <> l2  We next consider varying the frequency of the input wavelets and determine actual time s h i f t s by r e f e r r i n g to Figure 2 . 6 , which i s normalized to a c u t o f f frequency of 100 Hz.  Since P-Waves are found to Tie t y p i c a l l y between 200 to  2000 Hz, and S-Waves are t y p i c a l l y between 50 to 300 hz in s o i l s , we w i l l work w i t h i n these ranges. If our input wavelet has a frequency of 50 hz, then from Figure 2 . 6 , we have <t> =-25°.  From Equation ( 2 . 4 c ) , we also have  ?  Y  (f)  = -25° = -0.436 rad  At maximum amplitude, the cosine term in Equation (2.4b) should equal 1 . 0 ; thus 2-rrft + y ( f )  =0  17  which implies that a time s h i f t of  t  2^m  =  =  l  -  3 9  msec  -  If we consider a 150 hz s i g n a l , Figure 2 . 6 , indicates that y(f)  = *  2  = -55° = - . 9 6 0 rad  which corresponds to a time s h i f t of t = 0.960/(2*(250) = 1.02 msec Thus, between 50 and 150 hz, we have a time s h i f t d i f f e r e n c e of 370 usee. Figure 2.7 i l l u s t r a t e s the delay in the shear wave a r r i v a l f o r applied c u t - o f f  different  frequencies.  An a d d i t i o n a l consideration i s that the reverse p o l a r i t y method demands that in picking a crossover, a consistent method must be used.  For instance,  one may decide that the second crossover w i l l be used throughout the seismic investigation.  But i t has been found that in numerous cases t h i s may prove to  be a d i f f i c u l t c r i t e r i o n to implement due to the signal being masked by low frequency noise as i l l u s t r a t e d in Figure 2.8 (where a 300 hz low-pass was a p p l i e d ) . Figure 2 . 9 .  filter  Ideal wavelets for determining crossovers are i l l u s t r a t e d in Therefore,  since d i g i t a l f i l t e r i n g allows one to band-pass a  desired frequency range ( i . e . , S or P-Wave) without d i s t o r t i o n in and phase s h i f t , and the c r o s s c o r r e l a t i o n function makes use of a l l the information in the signals (averaging out i r r e g u l a r i t i e s and putting s i g n i f i c a n c e on dominant responses), t h i s method should provide better and more r e l i a b l e r e s u l t s than the reverse p o l a r i t y method.  Details of t h i s new approach w i l l be discussed  in the next chapter.  18  Low Pass F i l t e r McDonald's Farm 8.0 m Depth  s h e a r wave a r r i v a l s "  (delay  due  Co p h a s e  Cut-off. frequency  •I  = 20 Hz  shift)j  1  1  ^ %  I  20  40  .IWVIAAA.  GO  30  Time  i g u r e 2.7.  100  120  500  H:  600  Hz  1000  Hz  140  (msec)  E f f e c t of D i f f e r e n t F i l t e r C u t - o f f F r e q u e n c i e s on A c c e l e r o m e t e r Responses (Laing, 1985)  19  I  1  1 30.0  Figure 2.8.  1  1 60.0  1  1  I  L_  90.0  Time (msec) Seismic Data Acquired from Annacis Island Vibro-Compaction Site on May 18, 1988. Data Recorded with Accelerometer having a 300 Hz Low Pass Filter Applied. Diagram Illustrates Difficulty in Obtaining Crossovers due to Signal being Masked bv Manv Dominant low F r e a u e n c i e s .  Oppositely Polarized SH-Waves at 2.1m  44.61  20 Figure 2 . 9 .  48.89  Oppositely Polarized SH-Waves at 3.1m  53.77  Oppositely Polarized SH-Waves at 4.1m  40  60  80 120 140 100 110 Time (msec) Seismic Data Acquired from McDonald Farm on May 14, 1985. Data Acquired with Geophone with Natural Frequency of 28 Hz, and 100 Hz Low Pass F i l t e r Applied. The Diagram I l l u s t r a t e s Ideal Signals for Determining Crossovers  III. A.  APPLICATION OF FREQUENCY DOMAIN FILTERING  Discussion of the Crosscorrelation F i l t e r Formulation In order to obtain an appreciation of the discussed low pass frequency  f i l t e r , phase l a g s , and time s h i f t s , one should have a fundamental idea of what i s meant by d i g i t a l signal f i l t e r i n g and the corresponding possible side e f f e c t s .  Frequency f i l t e r i n g refers to obtaining a desired frequency  spectrum from a signal which i s corrupted with many d i f f e r e n t components. Most f i l t e r s permit the s e l e c t i o n of the upper and the lower l i m i t s of the passband.  This concept i s i l l u s t r a t e d in the following diagram: A(w)  X(u))  YD  OJ +> «=c  f f o l T  Figure  T  3.1 .  r  f f 2 3  Passband  T  f  o  f  l  Frequency  I l l u s t r a t i o n of F i l t e r i n g Desired Frequencies  It i s u s u a l l y possible to select the sharpness of the c u t o f f (the rate 9 at which the g a i n , |Y(")| , decreases as i t leaves the passband). Typical f i l t e r response curves are s p e c i f i e d by t h e i r cutoff frequencies, that i s , the frequency at which the gain has dropped by 3 dB (30% of amplitude, 50% of power). F i l t e r design i s c r i t i c a l l y dependent on the sampling rate of data acquisition.  The sampling theorem defines the minimum sampling i n t e r v a l one  can apply data a c q u i s i t i o n without l o s i n g a wavelet's i d e n t i t y  (aliasing).  In general, no information i s l o s t by regular sampling provided that the sampling frequency i s greater than twice the highest frequency component in the waveform being sampled.  The above concept i s i l l u s t r a t e d in Figure 3 . 2 ,  where we have a 4 ms sampling i n t e r v a l f o r wavelets of 75 hz, 175 hz and 250 hz frequency components.  Note that only (.a) in Figure 3.2 gives a true  response since the signal of 75 hz i s being sampled at 250 hz or about 3 1/2 times f a s t e r than the s i g n a l . 22  (c) Figure  3.2.  Sampling and A l a i s i n g D i f f e r e n t Frequencies Sampled "at 4 ms Intervals (250 Times Per Second), (a) 75 Hz "Signal; (b) 175 Hz Signal Yields Same Sample Values as 75 Hz; (c) 250 Hz Signal Yields Samples of Constant Value (0Hz) (R.E S h e r i f f and L. P. Gel d ar t , 1983)  Half the sampling frequency i s represented by the Nyquist frequency, f N , where fN=l/(:2At),- At being the sampling i n t e r v a l .  Frequencies which are  greater than fN by the amount Af w i l l be i n d i s t i n g u i s h a b l e from the lower frequency fN-Af.  In summary, one can state that wavelet  recovered exactly from the sampled values provided g(t) frequencies higher than the Nyquist frequency. transform of g ( t ) ,  G(f),  g(t)  can be  does not contain  When taking the Fourier  i t w i l l only contain information about g(t) 23  for  frequencies within the i n t e r v a l -fN to fN. G(f)  The Fourier s e r i e s w i l l  repeat  in each i n t e r v a l of width 2fN. Having discussed the importance of the sampling r a t e , i t i s useful to  further expand on f i l t e r i n g concepts. from a source wavelet.  F i l t e r s pass the desired information  Usually f i l t e r s attempt to remove s p e c i f i e d frequency  components of a sampled wavelet, which sometimes r e s u l t s in the f i l t e r  alter-  ing the amplitude and/or phase spectra of signals which pass through i t .  This  a l t e r i n g r e s u l t s because a f i l t e r i s defined by i t s own amplitude and phase spectra j u s t l i k e a wavelet.  When these two i n t e r a c t , they produce an output  signal which i s a product of the i n t e r a c t i o n of the wavelet and f i l t e r . output i s then defined by i t s own amplitude and phase spectra.  This process  i s i l l u s t r a t e d on page 2-8 and repeated in Figure 3.3 for convenience.  A(uv) = |A( )|e Input w  j^2 Y(ai)=|YU) |e• '  -  X(cu) Output  where X(o>)  Phase Spectra  A( ) w  Output  |AH||Y(<o)|e'  i=l,2  Amplitude Spectra Angular Frequency  Figure  3 . 3 . Changes in Phase and Amplitude Spectra Resulting from Interaction of F i l t e r and Input  24  This  Frequency f i l t e r s are generally c l a s s i f i e d as low-pass, high-pass band-pass, r e s p e c t i v e l y , as they d i s c r i m i n a t e against frequencies above or below a c e r t a i n l i m i t i n g frequency or outside of a given band of frequencies. An ideal f i l t e r i s one which has a gain of unity (|g("0| = 1) f ° within the passband and zero f o r those within the stopbands.  r  frequencies  Moreover,  it  i s advantageous that the f i l t e r produce zero phase s h i f t , so that the phase s p e c t r a , <f> (f), i s zero for a l l frequencies.  " I d e a l " f i l t e r s of these types  2  are the f o l l o w i n g Low-pass G ( ) L  W  = +1, |u)|<|a>| 0  0, | |<| | u  High-pass G U ) H  = 0, + 1,  Band-pass G U ) g  Uo  |co|<|a> | |u|>|a> |0  = +1, | 0, 1 I >(i>  or o» | ^21  a) = Angular Frequency These f i l t e r s are discontinuous at U K , W , , and w . o l 2 The low-pass f i l t e r i s obtained by applying a box car function with a 0  c e r t al i n passband oi .  The box car function i s i l l u s t r a t e d in Figure 3 . 4 , with  i t s associated transform. _-«/:» \  sine {  I  -8»/o  - 6 tia  where  Figure 3 . 4 .  2* la  O  sinc(.5at) =  Hn/a  2siri(;5at) at  Example of a Box Car Function, G(ai), and Its Transform, Sine (.5at) (R. E. S h e r i f f & L. P. Geldart, 1983) 25  The low-pass f i l t e r i s defined by the f o l l o w i n g expression: G (u) L  = box w (w) 2  (u /Tr)sinc(ai t) 0  0  0  = g (t)  (3.1)  L  where «- -»• denotes the Fourier transform. For d i g i t a l f u n c t i o n s , provided | u | < o  = ir/At,  = Angular Nyquist  frequency  L §\ = (1/ir) £ n=-L where  w  o  sinc(nw At)  (3.2)  0  L* -• «  It can mathematically be proven that the expression of the high-pass f i l t e r i s e s s e n t i a l l y the same as that f o r the low-pass f i l t e r , where g ( t ) = - g ( t ) , H  vided both f i l t e r s have the same c u t o f f frequency <d .  L  pro-  In deriving the passband  f i l t e r one j u s t manipulates the high and low-pass f i l t e r s to i n t e r a c t where the passband l i e s within the desired l i m i t s . The r e s u l t of using a f i n i t e series f o r gj; i s to introduce r i p p l e s , both within and without the passband, the e f f e c t i s e s p e c i a l l y noticeable near the cutoff frequency.  This i s referred to as the Gibbs' Phenomenon.  Gibbs'  Phenomenon describes the d i s t o r t i o n s which occur at d i s c o n t i n u i t i e s w i t h i n filtered signals.  This phenomenon can be conceptualized i n the time domain.  A signal truncated due to a short recording time, i s analogous to m u l t i p l y i n g the signal with a box car function which r e s u l t s i n d i s t o r t i o n s or " r i n g i n g " in the frequency domain.  As the s i g n a l ' s recording time increases " r i n g i n g "  becomes less prevalent because more information of the signal i s transformed and the f i n i t e series expressed i n Equation (3.2) i s c a r r i e d out f u r t h e r . These concepts are i l l u s t r a t e d i n Figure 3.5(b) where a 100 hz cosine signal i s being multipled with a box car f u n c t i o n .  Figure 3.5(a) i l l u s t r a t e s the  " r i n g i n g " e f f e c t i n the frequency domain. Figure 3.6 gives an i l l u s t r a t i o n of the harmonic makeup of s p e c i f i c waveforms.  In these figures a rectangular wave i s constructed by adding a suc-  cessive number of harmonics to the fundamental frequency.  As i s i l l u s t r a t e d ,  the more harmonics added, the better i s the approximation to a rectangular wave and/or a box car f u n c t i o n . 26  (a)  100  ZOO FREQUENCY,  300 Hertz  400  500  1.0  (b)  0.5  to g  0.0  -0.5  1.0 /SO TIME (N = 200,DELT=.00l)  Figure 3 . 5 .  100 , 150 NOTE:FREQUENCY = 100 HZ.  200  I l l u s t r a t i n g Gibbs' Phenomenon (a) Frequency Spectrum and (b) Truncated Signal  27  (a)  Rectangular Wave Produced by Adding Four Harmonics to the Fundamental. Shape of Composite Wave is Determined by Relative Amplitudes and Phases of Harmonics.  ripple  (b)  Rectangular Waves Produced by Combining 100 Harmonics. Note Reduction in Ripple.  Figure 3.5.  Composition of Rectangular Wave from Harmonics (after 1983). 28  0. W. Stimson,  Having discussed the problems associated with designing f i l t e r s , one should now explore e x i s t i n g f i l t e r s which lessen or remove the s e v e r i t y of undesired f i l t e r i n g e f f e c t s . i s the Butterworth Butterworth  A f i l t e r which has many advantageous q u a l i t i e s  filter.  Filter  The Butterworth f i l t e r i s a common form of low-pass f i l t e r , and i t can be defined by |G( )| W  2  = 1/ {.l+(o,/ ) Uo  (3.3)  2N}  where u> i s the " c u t o f f " frequency and N determines the sharpness of the c u t o f f , o p Curves of |G(w)| f o r various values of N are shown i n Figure 3 . 7 .  0  1  2  3  W.'oJ  0  normalized frequency Figure 3 . 7 . Amplitude Response of a Butterworth F i l t e r (R. E. S h e r i f f & L. P Gel d ar t , 1983)  29  The advantages associated with the Butterworth f i l t e r are as f o l l o w s : o  Their t r a n s f e r functions are smooth and maximally f l a t both inside and outside the passband.  o  The squared f i l t e r ( i . e . , the input i s f i l t e r e d twice so that the 2 amplitude response i s |G(w)| ) produced  zero phase s h i f t and i t s  power i s down by 3 dB ( f a c t o r of 1/2) at the c u t o f f frequency. (The c u t o f f frequency determines the half-power point of the f i l t e r ) . The value of N s p e c i f i e s the rate of attenuation where a larger value of N gives a greater rate of attenuation and " r i n g i n g " e f f e c t ( i . e . , Gibbs' Phenomenon).  An example of a c a l c u l a t i o n to determine the value of N necessary  i s given below. If the required a t t e n u a t i o n , measured in units of power, at 2 w  o  i s 48  db down from attenuation at oj =0, what value of N . i s necessary for a Butterworth o  filter?  From a d e f i n i t i o n of the decibel scale (decibel level=10 l o g ^ (power/  reference power) and using Equation 3 . 3 , we have w/ui  = W  o  10 l o g  (1/(1+W )) 2N  1 Q  W=0  /(1/(1+W ) 2N  N=2  = 48  so that (1/(1+2 )) = 1 0 " 2N  4  8  Approximately then (1/2 ) = l O " 2N  4  8  and N = 4.8/2 log 2 = 4.8/2(.301) = 8 Therefore N=8 i s required to obtain the required attenuation.  In general the  power changes by 6 db for a unit change in N. An important aspect that must be taken into account when applying the Butterworth f i l t e r i s the a l i a s i n g problem.  When one approaches the angular  Nyquist frequency of the sampled s i g n a l , a l i a s i n g errors may r e s u l t when the 30  t r a n s f e r function of Equation 3.3 i s a p p l i e d .  This i s conceptually i l l u s t r a t e d  in Figure 3.8 where a Fourier transform.of a signal with Nyquist of ui^ i s shown. It i s noted from Figure 3.8 that the Fourier Transform of t h i s function i s p e r i o d i c with period  2TT/W^.  A  Figure 3.8.  Fourier Transform of a Signal with Nyquist to  If we t r y to f i l t e r frequencies higher than the Nyquist ( i . e . , apply a cutoff frequency higher than the Nyquist), higher frequencies w i l l f o l d into the f i l t e r e d s i g n a l .  This i s i l l u s t r a t e d in Figure 3 . 9 .  - - Butterworth Window u i - C u t o f f Frequency  CD 3  -a  43 1.0  to  N  (0  C  to  Angular Frequency  Figure  3.9.  I l l u s t r a t i o n of the A l i a s i n g Problem When F i l t e r i n g ~~ 31  In Figure 3 . 9 , frequencies w i t h i n w* are incorporated into the s i g n a l . In order to avoid t h i s problem, a b i l i n e a r Z transform i s a p p l i e d .  The b i l i n e a r  Z transform converts a continuous t r a n s f e r function i n t o one which can be used on sampled data.  Furthermore i t eliminates a l i a s i n g errors by applying the  f o l l o w i n g frequency transformation:  w  d = | r tan ^  -Tr/At < (o < Ti/At  (3.4)  where bid = transformed or deformed angular frequency used to c a l c u l a t e the b i l i n e a r Z transform to = angular frequency i n the o r i g i n a l t r a n s f e r function Equation (3.3) This frequency transformation, wd, causes warping of the input frequencies when they approach the Nyquist frequency.  This warping i s i l l u s t r a t e d in  Figure 3.10.  1  •  tut  1  /  1  /  i 1  /  1 1 1  *  1 1  deformed frequencies  1  actual frequencies a f t e r a p p l i c a t i o n of b i l i n e a r Z transform  /  )  /  ,  Nyquist frequency  I  1  jd - d i s c r e t i z e d frequency  1  -ws  |  1  ' /  1 1  '  /  r  /  i i  / /  Figure 3.10.  I  1 1  1 1 1  Relationship Between the Deformed Angular Frequencies md, f o r the B i l i n e a r Z Transform and the Actual Frequencies, m, of the Digital F i l t e r (E. R. Kanasevich, 1981).  32  The Butterworth f i l t e r , when constructed in the time domain, i s a p h y s i c a l l y r e a l i z a b l e , recursive f i l t e r .  P h y s i c a l l y r e a l i z a b l e (causal) f i l t e r s are  those which do not demand future input values to compute output values.  A  recursive f i l t e r i s one which feeds back part of the output to compute new output. 3.3  The Butterworth f i l t e r having the t r a n s f e r function given by Equation  has the f o l l o w i n g time-domain r e l a t i o n s h i p N  N  Y = y " c X iTo " n  n  n  + m  BY  t l  m  n  (3.5) "  m  where the index N determines the order of the Butterworth  filter.  The Butterworth t r a n s f e r function G(f) can be constructed as a product (cascade) of second order (N=2) G2(f) and f i r s t order (N=l) G l ( f ) filters.  Butterworth  For instance, an eighth order Butterworth f i l t e r would be defined by  the f o l l o w i n g expression G(f;8) = G2(f)xG2(f)xG2(f)xG2(f)  (3.6)  Use of cascade f i l t e r i n g s i m p l i f i e s the algebra when designing the f i l t e r and reduces round-off e r r o r r e l a t i v e to a brute force c a l c u l a t i o n of G(f) f o r each order (Reference 14).  Thomson and Chow (1980) discuss how the Z.-domain t r a n s -  f e r functions of the f i r s t and second order Butterworth f i l t e r s are obtained, and determine the corresponding time-domain c o e f f i c i e n t s .  The time-domain  algorithms f o r the f i r s t order (N=l) and second order (N=2) Butterworth f i l t e r s are,  respectively.  Y  n  Y  n = o n  =  d  c  o n X  X  +  +  l n-l  d  X  c  l n-l X  +  +  e  c  l n-l  < ' ) 3  Y  l n-2  where X  n  = input  Y  n  = f i l t e r e d output  33  X  +  b  l n-l Y  +  b  2 n-2 Y  7  The Butterworth f i l t e r has zero phase lag and then back fed into the f i l t e r .  because a waveform i s front fed  C l e a r l y , i f there i s a phase s h i f t at a  given frequency in the f i r s t pass through the f i l t e r , there w i l l be an equal phase s h i f t of opposite sign at the frequency in the second pass.  By process-  ing the data in reversed chronological order through the f i l t e r , the two phase s h i f t s cancel e x a c t l y . Another consideration to take into account when processing the d a t a , i s the . f i n i t e length of the time s e r i e s .  This f i n i t e length can be viewed as  a t r u n c a t i o n , . where we m u l t i p l y an i n f i n i t e length s e r i e s by a box car f u n c t i o n . This m u l t i p l i c a t i o n r e s u l t s in Gibbs was previously discussed. 3.11 and 3 . 1 2 .  1  Phenomenon in the frequency domain, as  These concepts are again i l l u s t r a t e d in Figures  In Figure 3.11 we have an ideal i n f i n i t e length time series  with a 60 hz frequency content.  The Fourier transform of t h i s signal gives  ideal Dirac deltas at -60 and +60 hz in the frequency domain.  Figure 3 . 1 1 .  An I n f i n i t e l y Transform  Long Sinusoidal and I t s '  34  Fourier  Figure 3.12 i l l u s t r a t e s the frequency spectrum obtained from a time series truncated by a box car function in the time domain. g(t) A  G(t)  -60  Figure 3.12.  -Jr-  time, t  60 frequency, f  A Truncated 60 Cycle Signal and i t s Fourier Transform. The Signal was Switched on f o r " IT/2 Second  In order to remove the "ringing" e f f e c t , the data i s tapered.  The time  series can be tapered by a pair of cosine b e l l s , where the f i r s t and l a s t points of the trace approach the mean value of the s e r i e s , which i s u s u a l l y zero. Equation 3.9 gives us the cosine b e l l s , and Figure 3.13 i l l u s t r a t e s the tapering window. (1 + cos w.  ir(n+L)  )/2  1  (L+M)<n<-L Lin l L  (1 + cos Tr(n-L) )/2  (3.9.)  L < n < L+M  where M determines the period of the cosine b e l l ( i . e . , the r o l l - o f f of the window), and 2L i s 80% of the length of the time s e r i e s .  Hamming (1977)  recommends that M be about 10% of the e x i s t i n g data with 80% i n the f l a t part of the window.  The data can then be padded with zeros.  In t h i s way, there  are no d i s c o n t i n u i t i e s to i n i t i a t e transients during the Fourier transform.  35  w n  max of 80% | | <- A1 , sampl i ng  i  -(L+M)  Figure 3.13.  i  I  AX  j  a  k  A i  A"  L+M L d i s c r e t e time  -L  A Data Window Used to Taper a Discrete Function  Crosscorrelation The f i r s t step in determining the i n t e r v a l times, once seismic wavelets are adequately f i l t e r e d , i s to apply c r o s s c o r r e l a t i o n s .  The c r o s s c o r r e l a t i o n  of the successive wavelets i s defined as (3.10)  where x i s the time s h i f t between the two time s e r i e s . provides information between the two wavelets.  The c r o s s c o r r e l a t i o n  In t h i s case we are looking  f o r the a r r i v a l time differences between the dominant responses of the recorded wavelets. when the value of  With the above c o n s i d e r a t i o n , one can postulate that $  X  y ( . ) i n Equation 3.10 i s maximum, we have time s h i f t T  which i s representative of the time i n t e r v a l AT as applied from Equation 2,2a on page 2 - 2 . As stated previously, *  (x)  i s - t h e value which, represents the d i f f e r e n c e  in the a r r i v a l times of the dominant responses of the wavelets.  Therefore  t h e i r c r o s s c o r r e l a t i o n i s e s s e n t i a l l y an a u t o c o r r e l a t i o n i f the wavelets are the same with maximum value at x.  This concept i s i l l u s t r a t e d i n Figure 3.14.  36  Process  Autocorrelation  Figure 3.14.  Autocorrelation of a Process  An important aspect to consider when applying c r o s s c o r r e l a t i o n s to seismic wavelets i s D.C. s h i f t s . i s not centred at zero mean. two data sets X and t  D.C. s h i f t s occur when the recorded signal  Equation 3.10 defined the c r o s s c o r r e l a t i o n of  where t i s the displacement of Y r e l a t i v e to X^. t  This equation i s generally applied to the s i g n a l s when they have zero mean (which corresponds to zero D.C. o f f s e t ) .  A simple example i s i l l u s t r a t e d by  c o r r e l a t i n g the two f u n c t i o n s , X = { l , - 1 , - 1 / 2 } and Y = { 1 , 1 / 2 1 / 2 } , shown i n t  t  Figure 3.15.  (a) S h i f t of - 2  (b) S h i f t of - 1  (c) S h i f t of 0  i  +1/2 '  X! Z3 Q. E  QJ  4->  3  -fa-  0  > +1  QJ  [•1/2  0  Time  Time  -It  Time -3 •1/2 ^  CL  i  " +1/2  QJ  -a  4->  V  Time - -1  * +1 + 1/2  Time -1/2  0  E  E  <f> (-2)=0+0+l/2 = 1/2 xy  Figure 3.15.  P  x y  +1/2  0  • -1 +1  -a  CL  E  k  >+l iL  QJ  T3 3 +->  CL  E  ( - l ) = 0 - l + l / 4 =-3/4  +1/2 'Time  0 1 4,  , +1/2  (0) = l - l / 2 - l / 4 =  C a l c u l a t i n g the Crosscorrelation of Two Functions  37  +1/4  At the end of the signals the functions are padded with zeros in order that they w i l l provide no c o n t r i b u t i o n to c o r r e l a t i o n values. Let us now consider the two s i g n a l s , shown in Figure 3 . 1 6 , with negat i v e D.C. o f f s e t d, and d , r e s p e c t i v e l y . 0  CD  -a 3  -o  Time  Time  Q. E  CL  -AAMA-  E  I — W W  1  (b)  (a)  Figure 3.16.  T  Signals with D.C. S h i f t s d^ and d,.  Before taking the c r o s s c o r r e l a t i o n of these s i g n a l s , they are f i r s t padded with zeros, which i s i l l u s t r a t e d in Figure 3^17.  4 Padding  CD  -a 3  Padding  CD  /  -a  ime  CL  (b)  (a)  Figure  E  3.17. Padding Signals with Zeros at the End of Their Time Series 38  Time  When taking the c r o s s c o r r e l a t i o n of the f u n c t i o n s , shown in Figure 3 . 1 7 , we e s s e n t i a l l y have the autocorrelation of a box car f u n c t i o n , which i s i l l u s t r a t e d in Figure 3.18.  time s h i f t , x Figure  3,18. Autocorrelation of Box Car Function  Therefore, in order f o r the signals to give representative c r o s s c o r r e l a t ion  values, they should f i r s t be adjusted so that they are centred at zero  offset. Another way of looking at the D.C. s h i f t problem i s to consider Equat i o n .3.10. An o f f s e t can be represented by s h i f t i n g one signal r e l a t i v e to the other, that i s X^=X +C where C i s a constant. k  In t h i s case Equation 3.10  becomes  V > £ kkx T  =  X  Y  = £ (  V  T )  s  V  x  T )  ( 3  +  + c k  +  )  Y  l l a  ( 3  k + x  ? i<^ CY  (3  -  llb  llc)  Equation (3.11c) c l e a r l y shows that the D.C. s h i f t would r e s u l t in a m i s representative c r o s s c o r r e l a t i o n value. For a p r a c t i c a l example, consider the two wavelets shown i n Figures 3.19 and 3.20 with D.C. s h i f t s of - 0 . 0 6 and - 0 . 2 , r e s p e c t i v e l y .  39  )  )  -I  0  I  I  I  I  SO  100  150  200  Tl»F. UNIT..  Figure  3.19. Signal with D.C. Offset of -0.20  -0.08  Figure  3.20. Signal with D.C. Offset of -0.06  40  Taking the c r o s s c o r r e l a t i o n of these two s i g n a l s , we obtain the t r i a n g l e function i n Figure 3.21 (with p o s i t i v e time s h i f t shown).  40  60  POSTIVF:  80 TIME  SHIFT  100 UNITS...  Figure 3.21. Crosscorrelation of Offset Signals If we correct f o r D.C. o f f s e t s , the signals take the form shown i n Figure 3.22 and 3 . 2 3 , and t h e i r c r o s s c o r r e l a t i o n function i s shown i n Figure 3.24.  60  80 TIME  Figure 3.22.  UNIT  Signal 1 Corrected f o r D.C. Offset with the Number of . Points Changed from 200 to 150 Because Dominant Reponse Occurs Between 40 to 60 Time Units 41  0.05 0.04  60  80 TIUF UNIT  Figure 3.23. Signal 2 Corrected f o r D.C. Offset  0.06  0.04V  o.oz  0.0  -0.02Y  -0.04 V  -0.06 40  60  80  100  120  TIME SHIFT UNITS  Figure 3.24. Crosscorrelation Value f o r D.C. S h i f t Corrected S i g n a l ?  42  The maximum of t h i s c r o s s c o r r e l a t i o n occurs at 60 units of time s h i f t s , which is expected by v i s u a l l y picking time s h i f t or lag in Figure 3.25.  The procedure in obtaining the desired v e l o c i t y p r o f i l e can thus be summarized as f o l l o w s : (1).  Record data  (2)  Remove D.C.  (3)  Taper time-domain data with cosine b e l l s  (4)  Apply FFT to obtain frequency spectrum  (5)  Pick frequencies to be f i l t e r e d (low & high in the bandpass)  (6)  Apply Butterworth  (7)  Take inverse FFT  (8)  Take c r o s s c o r r e l a t i o n of f i l t e r e d time-domain s i g n a l s  (9)  Determine t , time o f f s e t  (10)  shift  filter  Obtain v e l o c i t i e s  The program. CROSSCOR was developed to apply the above steps.  CROSSCOR  i s a graphics i n t e r a c t i v e program which displays the frequency s p e c t r a , unf i l t e r e d and f i l t e r e d t r a c e s , and c r o s s c o r r e l a t i o n plots on a mainframe graphics 43  terminal and has been implemented on the IBM - P.C. perform Steps 2 through 9. (2)  SUBROUTINE REMAV  (4) & (7)  SUBROUTINE FASTF  (3)  SUBROUTINE SMOOTH  (5)  SUBROUTINE PL AND PLA - MAIN FRAME  (6)  SUBROUTINE BNDPAS  (8)  SUBROUTINE CROSS  (9)  SUBROUTINE MAXSN  Additional subroutines implemented in CROSSCOR are: (10)  SUBROUTINE NORME NORME normalizes a data series by i t s RMS energy, in order to bring the c o n t r i b u t i o n of each trace equal in magnitude when crosscorrelating.  (11)  SUBROUTINE P0W.2 P0W2 determine M in N=2**M f o r a s p e c i f i c trace length ( i . e . ,  L)  and then pads the rest of the t r a c e , L to N, with zeros.  The value  of N i s c r i t i c a l when implementing the FFT, because i t requires the length of the trace to be a power of two.  44  The f o l l o w i n g subroutines  SUBROUTINE FXFFT This subroutine saves computing time,.because i t allows for only one c a l l to FASTF.  In FASTF  one trace i s input as a real component and second trace i s input as an imaginary components. Reference 5 f u r t h e r expands on t h i s (Chapter 3, Section 3.6)  45  CROSSCOR was tested f o r performance by superpositioning several sinusoidal waves and f i l t e r i n g out the desired frequencies.  In Figure 3.26 there i s an  i d e n t i c a l traces containing fequencies of 30, 50, and 100 Hertz. are o f f s e t from each other by 20 time  units.  The s i g n a l s  The signals are then bandpassed  between 40 and 80 Hertz r e s u l t i n g in only the 50 Hertz signal remaining as i l l u s t r a t e d in Figure 3 . 2 6 ( c ) .  The c r o s s c o r r e l a t i o n function i s i l l u s t r a t e d  in Figure 3.26(d) where negative and p o s i t i v e time s h i f t s are i l l u s t r a t e d . determining the a r r i v a l time differences between the successive wavelets,  In it  is only necessary to obtain the maximum p o s i t i v e c o r r e l a t i o n value at a p o s i t i v e time s h i f t .  This i s due to the fact that the c r o s s c o r r e l a t i o n function i s com-  paring a r r i v a l time differences of the dominant responses ( i . e . , maximum p o s i t i v e c o r r e l a t i o n value) and the wavelet a r r i v a l time ( p o s i t i v e time s h i f t ) .  from the deeper layer w i l l always have a longer The c r o s s c o r r e l a t i o n function from the  f i l t e r e d signals in Figure 3.26(c) gave a correct maximum time s h i f t of 20 u n i t s .  46  Time Units  Time Units  g  (b) H  •J T  O.Ql  1  v 100 200 300 400 500 Frequency, Hz  100 200 300 400 500 Frequency, Hz  Time Units  Time Units CROSSCORRLA TION  (d)  -tOO  -150  -SO  0 SO TIUE SHIFT  100  Figure 3.26 Processing Synthetic Data with CROSSCOR. (a) Input Trace ~[b) Frequency Spectrums ( c l Bandpassed I races (4U to 80 Hz) (d) Crosscorrelation hunction of the F i l t e r e d Traces. 47  B.  F i e l d Program and Discussion of Results In t h i s section CROSSCOR i s applied to real data to assess i t s performance  and r e l i a b i l i t y .  Results from f i v e research s i t e s are given, and the v e l o c i t y  p r o f i l e s from d i f f e r e n t data and methods of data reduction are presented and compared.  The research s i t e s are i l l u s t r a t e d in Figure B . l and are as f o l l o w s :  1)  McDonald Farm S i t e  2)  Grant McConachie Way S i t e at Laing Bridge  3)  T i l b u r y Island Natural Gas Plant S i t e  4)  Annacis Bridge P i l e Research S i t e  5)  Lower Langley 232nd S t . Site on Trans Canada Highway  The McDonald Farm Site data was acquired on June 25 and July 2, 1987.  In  t h i s i n v e s t i g a t i o n f i v e sets of data were compiled and two d i f f e r e n t sensors were used.  The f i r s t sensor was the Hogentogler Super Cone with a Geospace*  v e l o c i t y transducer and the second sensor was the UBC Seismic Cone Pressuremeter, which had two ( P i e z o e l e c t r i c bender element) accelerometers i n s t a l l e d . The seismic data ( i . e . , hammer shear source) was reduced by both the reverse p o l a r i t y method and CROSSCOR.  In t h i s research s i t e , CROSSCOR was found to be  very helpful in respect to reducing noisy accelerometer data where picking c r o s s overs was very d i f f i c u l t or impossible. The Grant McConachie Way S i t e data was acquired on November 17, 1987, where the UBC #7 Seismic Cone was used with a P i e z o e l e c t r i c bender element a c c e l e r o meter sensor i n s t a l l e d .  Two d i f f e r e n t types of sources were used in t h i s  i n v e s t i g a t i o n , i . e . , the hammer shear source and the P-plate source.  From the  frequency spectrums obtained there were three d i s t i n c t responses of the a c c e l e r o meter observed. at 100 Hz,  The f i r s t response was due to the shear wave (S-wave) occurring  second response was a compression wave, (P-wave), at 900 Hz , the  t h i r d response resulted from the resonanting accelerometer at 3000 Hz. The T i l b u r y Island Natural Gas Plant Site data was acquired on July 16, 1987, where the UBC Seismic Cone Pressuremeter was used with a shear source.  In t h i s  research s i t e i t was not possible to use the reverse p o l a r i t y method due to the very noisy responses.  However, CROSSCOR proved to be extremely helpful i n reduc-  ing and analyzing the time s e r i e s . * See Laing (1985). 48  U.S.A.  Figure B . l .  General Site Location Map  49  The Annacis P i l e Research S i t e data was acquired on January 8 , 1988 and February 10, 1988.  On January 8, 1988 the UBC Seismic Cone Pressuremeter was  used f o r data a c q u i s i t i o n , and on February 10, 1988 the UBC #8 Seismic Cone was a p p l i e d .  In these investigations shear v e l o c i t i e s were obtained using  shear sources.  Since both sides of the shear source were used, i t was p o s s i -  ble to obtain two v e l o c i t y estimates for one depth increment using CROSSCOR as opposed to the reverse p o l a r i t y and visual pick method which allows f o r only one v e l o c i t y estimate. The Lower Langley 232nd St. data was acquired on November 19, 1987, where the UBC #7 Seismic Cone was used f o r data a c q u i s i t i o n .  In this i n v e s t i g a t i o n  f i l t e r e d (analogue) and u n f i l t e r e d reverse polarized hammer shear source data was obtained.  In a d d i t i o n , point source Buffalo Gun data was reduced with  CROSSCOR which resulted in separating and measuring both shear and compression wave v e l o c i t i e s . The analysis to follow in Section B i s as f o l l o w s : Section B . l : In t h i s section the research s i t e s investigated are presented along with t y p i c a l subsurface stratigraphy and cone penetration p r o f i l e s which give some i n d i c a t i o n of the geotechnical properties of s o i l s i n v e s t i gated. To obtain a more thorough understanding of the s i g n i f i c a n c e of p r o f i l e s presented, one should refer to the numerous UBC i n - s i t u Group reports. Section B.2: In t h i s section CROSSCOR i s compared to the established reverse p o l a r i t y method in obtaining v e l o c i t y estimates with low^-pass analogue f i l t e r e d data and s e l f - f i l t e r i n g geophone data. The purpose of t h i s i n v e s t i g a t i o n i s to prove that CROSSCOR works as well as the reverse p o l a r i t y method when reducing clean s i g n a l s , and CROSSCOR works better than the reverse p o l a r i t y method with signals containing many dominant low frequencies Section B.3: In t h i s section CROSSCOR i s investigated for i t s performance in removing P-waves and noise from u n f i l t e r e d accelerometer data. The v e l o c i t i e s obtained from these noisy traces are than compared to v e l o c i t i e s obtained from corresponding low-pass analogue f i l t e r e d data. Section B.4: In t h i s section CROSSCOR i s tested for i t s performance in obtaining v e l o c i t y estimates from nonpolarized sources ( e . g . , the Buffalo Gun Source and comparing one side of the Hammer Shear source to the other). This section i l l u s t r a t e s the f l e x i b i l i t y of CROSSCOR in working with d i f f e r e n t types of sources and giving more v e l o c i t y estimates f o r a downhole seismic p r o f i l e as opposed to the reverse p o l a r i t y method. Section B.5: In t h i s section CROSSCOR i s tested for i t s a b i l i t y to e x t r a c t compression wave events from seismic traces and give corresponding compression wave v e l o c i t i e s . Section B.6: In t h i s section specifying the Bandwidths ( i . e . , cutoff f r e quencies) when applying CROSSCOR i s addressed. 50  B.l  Research Sites  McDonald Farm Site McDonald Farm is situated upon Sea Island which is part of the Canadian DOT land at Vancouver International Airport.  As was shown in Figure B . l , Sea  Island is located between the North Arm and Middile Arm of the Fraser River. Therefore, the McDonald Farm Site is located on a river delta complex with a corresponding subsurface typical of this type of environment. The general geology consists of deltaic distributary channel f i l l and marine sediments (Armstrong 1978).  A typical CPT* profile from the site is  illustrated in Figure B . l . l and has the corresponding stratigraphy: 0-2m  soft organic silty clay  2-13m  loose to dense coarse sand; some layers of fine sand  13-15m  fine sand, some s i l t ; transition zone  15 > 80m soft normally consolidated clayey s i l t The groundwater table is generally found to l i e lm below ground surface. Additional information can be obtained from Campanella et a l , 1981.  * Cone Penetration Test  51  UBC S l t a Location! On  Site  Loci  FRICTION RATIO Rf (X) 0 5  I  M  U  SIT  T• E:  s ~r I  McOONALD FARM  CPT Data i  577-86-2  Cons Unodi  HOC SUPER STD PP  SLEEVE FRICTION (ba-) 0 t  CONE BEARING Ot (bar)  MC3  Pago Noi  88/09/25  CoB«ant«i  1 / 1 PIEZO CONE SNO  PORE PRESSURE U Ob of » a t « r ) 0 100  DIFFERENTIAL P.P.  RATIO -.2 0 0  AU/Ot  .8  «d soft orrparnc ay loose-d nse coarse and  10 (A L  fine with  ai  •P  Q) E  hd t  0_  LU  a  20  soft no mal l y consoli lated clayey ; i l t 30 Depth Incramont •  Figure B . l . l .  . 025 m  Max U*pth i  24. OS •  Cone Penetration P r o f i l e f o r McDonald Farm S i t e  Laing Bridge S i t e The Laing Bridge S i t e i s located at the eastern end of Sea Island adjacent to Laing Bridge approach f i l l at Grant McConachie Way as shown in Figure B . l . The research s i t e i s located next to over-pass embankments where extensive research has been conducted by Donna Le C l a i r (1988) into p r e d i c t i n g settlements. A t y p i c a l CPT p r o f i l e from the s i t e i s i l l u s t r a t e d in Figure B.1.2 and has the corresponding  stratigraphy:  0.2m  extraneous s o i l , top s o i l , and sandy s i l t y clay  2-20m  medium dense to very dense sands  20-60 to 65m  normally consolidated clayey  silt  The groundwater table i s generally found to l i e between lm and 1.5m below ground surface,  with f l u c t u a t i o n due to t i d a l i n f l u e n c e .  53  Si to On  Location. Slta  DIFFERENTIAL P.P. RATIO AU/Qt -.2 0 .8  T E S "r i  I M sin U -  UBC Loci  Laing  Br-idga  McConachla  S.  0/P  CPT  Cons  Data  i  Usadi  10/22/B7  UBC*7  Std/BFS  PP  Commantei  SLEEVE FRICTION (bar) 0 ' 1.5  CONE BEARING Qt (bar)  PORE PRESSURE U (a. of aalar) 0 175  NG  P a g o Noi  09i 3 0  1 / 1 CPTU-6 INTERPRETED PROFILE  FRICTION RATIO Rf (X) 0 5  sandy, silty clal medium dense to very dens sands 20  u.  normally consolida ed clayey s i t  0)  U1  E  a  Depth  40  40  4f>  Incrament  •  Figure B . 1 . 2 .  .025 m  Max  Depth  i  54. 775 m  Cone Penetration P r o f i l e f o r Laing Bridge Site  T i l b u r y Island Natural Gas Plant S i t e The T i l b u r y Island Natural Gas Plant i s located on T i l b u r y I s l a n d , as shown in Figure B . l .  The research s i t e i s located in the South Arm of the Fraser  River which has subsurface stratigraphy t y p i c a l of a r i v e r delta complex.  John  Howie (1987-88) conducted extensive work with the UBC Seismic Cone Pressuremeter at t h i s s i t e . A t y p i c a l CPT p r o f i l e from the s i t e i s i l l u s t r a t e d in Figure B . l . 3 and has the following  stratigraphy:  0-2m  Fill  2-25m  Sand  25-45m  Interbedded s a n d - s i l t layers  The groundwater table i s generally found to l i e lm below ground s u r f a c e , with f l u c t u a t i o n due to t i d a l  influence.  55  UBC Slto On  Location! Slto  Loci  PORE PRESSURE U ( » . of wotsr) 0 50  TILBURY  IM  SITU CPT  IS - L N G  Cono  T I L B U R Y if 1 SLEEVE FRICTION (bar) 0 1.5  Data  i  Usodi  CONE BEARING Qc (bar) .0  T E SB T I M G  31/03/B7  P a g o Noi  UBC SCP1  Commontoi  200  FRICTION RATIO Rf U> 0 5  I  PPi  / 1 BT&BS INTERPRETED PROFILE  DIFFERENTIAL P.P. RATIO AU/Oc -.2 0 .6  n  15 I/) L 01 01 E  0_ Ld  a  30  interbedded sand-sil 1ayers  45 Depth  Incromant  Figure B . 1 . 3 .  i  .025 m  Max  Dopth  i  39. 9 m  Cone Penetration P r o f i l e f o r Tilbury Island Site  Annacis P i l e Research S i t e The Annacis P i l e Research S i t e i s located on the extreme east end of Lulu Island which i s part of the Fraser River delta deposits.  This research s i t e ,  shown in Figure B . l , has been extensively investigated for geotechnical and geological parameters by Michael Davies (1987) and the B r i t i s h Columbia M i n i s t r y of Transportation and Highways (1984). The s u r f i c i a l geology of the Lulu Island region i s s i m i l a r to that of Sea Island and T i l b u r y Island ( i . e . , d e l t a i c d i s t r i b u t a r y channel f i l l and marine sediments).  A t y p i c a l CPT p r o f i l e from the s i t e i s i l l u s t r a t e d in Figure B . l . 4  and has the corresponding  stratigraphy:  0-2m  sand f i l l  2-15m  soft organic s i l t y clay  15-28m  medium dense sand with minor s i l t y sand lenses  28-60m  normally consolidated clayey s i l t with t h i n sand layers  57  UBC  I M  Slta Lacatlom  ANNA PLT  On  CPT PR1  3lt«  Loci  SITU CPT Data i Con. U » a d .  PORE PRESSURf SLEEVE FRICTION U U. of «<rtjr) Hw) 2.3 0 100 .0  CONE BEARING Ot (bar) 0 0  TESTING  850813 MO AS DV U0C0 STO TIP  200  Poga Noi Co«««nt.i  1 / 2 NEAfl CASING INTERPRETED PROFILE .  FRICTION RATIO DIFFERENTIAL P.P. Rf CD RATIO lU/Ot 0 3 . _0 . '  SAND f i l l  • 10"  soft oraanic silty CLAY  10-  19 L a CO  a E  a. LU a  -  20  med. dense SAND  20-  minor silty SAHp lenses  30 Dopth  Figure B.1.4.  Incr«»«nt  i  .025 n  Max Dap th i  35.075  spe  -  Cone Penetration P r o f i l e f o r Annacis Bridge P i l e Research S i t e  over  Lower Langley 232nd S t . The Lower Langley 232nd S t . S i t e i s located adjacent to the northern e x i t road on the westbound lane of the main Trans Canada Highway out of Vancouver which follows the Fraser V a l l e y .  The e x i t road i s the o f f ramp to 232nd S t r e e t .  This s i t e has been investigated f o r geotechnical and geological parameters by John S u l l y (1987) and James Greig (1985). The s u r f i c i a l geology of t h i s s i t e i s t y p i c a l of the Fraser Lowland area. A CPT p r o f i l e from the Lower Langley 232nd S t . S i t e i s i l l u s t r a t e d in Figure B . l . 5 and c o n s i s t predominantly of normally consolidated s e n s i t i v e s i l t y clay with occasional sand lense.  The near surface material i s overconsolidated due  to d e s s i c a t i o n .  59  UBC Site On  Location. Site  Loci  FRICTION RATIO Rf (Jt) 0 5  IN  SITU CPT  RGC-E8-AB  Cong  L0WER232 LANCLEY SLEEVE FRICflON (bar) 0 2  Data  i  Usedi  CONE BEARING Oc (bar)  TEST  11-19-B7  16iJ5  HOC SUPER  STD  u  PORE PRESSURE U (a. of later)  Page  IMG  Not  Comment«i  1 /  I  C77-8713  DIFFERENTIAL P.P. RATIO lU/Qc  5HHF1LT  INIERPREtEQ PROFILE  fill o.c. clay N.C. clay  N.C. clay with f i n e sand interbeddn N.C. clay  Interbedqed clay and snad  Depth  Increment  i  .025 m  Max O e p t h  •  2 9 . 655 m  * Overconsolidated **Normany consolidated Figure B . l . 5 .  Cone Penetration P r o f i l e f o r Lower Langley 232nd S t . S i t e  B.2  Cross-Over Method vs. CROSSCOR with Analogue F i l t e r e d Applied and D i f f e r e n t Seismic Probes Used In t h i s section the well established reverse p o l a r i t y or Cross-Over method  i s compared to the newly developed CROSSCOR in obtaining shear wave v e l o c i t i e s . The seismic traces have been recorded with the p i e z o e l e c t r i c bender element accelerometers (with a 300 Hz low-pass analogue f i l t e r applied) and the s e l f f i l t e r i n g Geospace geophone (Natural  Frequency = 28 Hz).  As was previously discussed in Section II, that in  picking a Cross-Over,  seismic p r o f i l e .  the Cross-Over method demands  a consistent method must be used throughout the  But i t was found that in some cases seismic traces may be  masked by several dominant low frequencies making Cross-Over picks d i f f i c u l t or impossible ( e . g . ,  Figure 2 . 8 ) .  The above problem i s futher i l l u s t r a t e d in Figure  B.2.1 with data recorded at Annacis I s l a n d .  In the discussion to f o l l o w , i t w i l l  be i l l u s t r a t e d that the Cross-Over method w i l l give i d e n t i c a l v e l o c i t i e s to that of CROSSCOR with clean s i g n a l s ( i . e . ,  containing one dominant frequency as was  i l l u s t r a t e d in Figure 2.9) and that CROSSCOR w i l l give more r e l i a b l e estimates with traces l i k e those of Figure 2.8 and B . 2 . 1 . McDonald Farm The f i r s t set of data presented in Table B . 2 . 1 i s from McDonald Farm, where the Hogentogler Super Cone v e l o c i t y transducer was used f o r data a c q u i s i t i o n on July 2, 1987.  Table B . 2 . 1 i l l u s t r a t e s the v e l o c i t i e s obtained from the  method and CROSSCOR.  Cross-Over  As i s shown, the v e l o c i t i e s are very comparable due to the  clean signals the slow responding v e l o c i t y transducer g i v e s .  One must keep in  mind that since CROSSCOR uses a l l the information contained in the seismic t r a c e s , the v e l o c i t i e s obtained from i t would be more accurate.  The v e l o c i t i e s from t h i s  set of data i s plotted along with the McDonald Farm cone bearing p r o f i l e in Figure B . 2 . 2 .  Figure B.2.2 i l l u s t r a t e s that the Cross-Over method and CROSSCOR  give v e l o c i t i e s that are very comparable and vary with the cone bearing p r o f i l e (i.e.,  increase i n v e l o c i t y with increase in bearing).  The second set of data presented i s also from the McDonald Farm S i t e , where the UBC Seismic Cone Pressuremeter (SCPM) was used for data a c q u i s i t i o n on June 25, 1987.  In t h i s i n v e s t i g a t i o n the signals recorded had several dominant f r e -  quencies present even with the 300 Hz low-pass f i l t e r a p p l i e d . 61  This character-  Time(msec) Figure B . 2 . 1 .  Seismic Data Acquired from Annacis Island Vibro-Compaction S i t e on May 18, 1988. Data Recorded with Accelerometer having a 300 Hz LowPass F i l t e r Applied. Diagram I l l u s t r a t e s D i f f i c u l t y in Obtaining Cross-Overs Due to Signal Being Masked by Many Dominant Low Frequencies.  Table B.2.1 Calculated Velocity Profiles from McDonald Farm Comparing the Cross-Over Method with CROSSCOR. The Data was Obtained with a Self-Filtering (Natural Frequency = 28 Hz) Geophone.  Depth (m) Interval  Cross-Over Hammer Source July 2, 1987 Hog-Super Vs (m-sec)  CROSSCOR Hammer Source July 2, 1987 Hog-Super Vs (m/Sec)  3.2  2.7-3.7  113  115  4.2  3.7-4.7  146  150  5.2  4.7-5.7  137  141  6.2  5.7-6.7  153  155  7.2  6.7-7.7  149  143  8.2  7.7-8.7  189  189  9.2  8.7-9.7  180  174  10.2  9.7-10.7  225  220  11.2  10.7-11.7  200  193  Average Depth (m)  63  Cone Bearing  Velocity ( m/s ) 100 •  t  i  ()  200  I  b  0  100  Qt (bar) 200  i—i—i—I—h  «~~~~ CROSSOVER CROSSCOR  CL  a  •  i  i  Figure B . 2 . 2 .  i—i—i—l  (a) V e l o c i t y P r o f i l e from McDonald Farm Comparing the r J , c c - ( W Method and CROSSCOR. The Data was Obtained with a S e l f F i l t e r i n g (Natural Frequency = 28 Hz) Geophone (b) Cone Bearing P r o f i l e from McDonald Farm.  64  i s t i c of these seismic traces i s t y p i c a l of the accelerometer transducer due to i t s high s e n s i t i v i t y .  \  The v e l o c i t i e s determined from the Cross-Over method  and CROSSCOR are given i n Table B . 2 . 2 .  The SCPM has accelerometers in both the  cone portion and the pressuremeter portion of the instrument.  The cone a c c e l e r o -  meter (SCPM CN) gives noticeable discrepancies between the two methods of v e l o city calculation.  But when these v e l o c i t i e s are compared to the SCPM pressure-  meter accelerometer (SCPM PM), the v e l o c i t i e s obtain from CROSSCOR are more similar. B.2.3(a).  The v e l o c i t i e s obtained from t h i s research s i t e are plotted in Figure Figure B.2.3(b) i l l u s t r a t e s the cone bearing p r o f i l e for McDonald  Farm and as i s shown, the v e l o c i t i e s in Figure B.2.3(a) vary with the cone bearing p r o f i l e . Laing Bridge S i t e The t h i r d set of data presented i s from the Laing Bridge S i t e , where the data was acquired on November 17, 1987 with the UBC #7 Seismic Cone. data i t was also not possible to accurately determine consistent  For t h i s  Cross-Overs,  due to the s e n s i t i v e nature of the accelerometer (as discussed in Section  I).  The shear wave v e l o c i t i e s obtained from the Cross-Over method and CROSSCOR are i l l u s t r a t e d in Table B . 2 . 3 .  Again we see noticeable discrepancies between the  two methods of v e l o c i t y determination (note depths 6 . 3 , 1 0 . 3 , 11.3 and 13.3m), but when compared to the v e l o c i t i e s obtained from the P-Plate source at 10.3 and 11.3m, CROSSCOR gives the most repeatable values.  The v e l o c i t i e s obtained from  t h i s s i t e are plotted in Figure B.2.4 along with the cone bearing p r o f i l e . The l a s t set of data presented in t h i s section was obtained from the Lower Lang ley 232nd S t . on November 19, 1987, where the UBC #8 Seismic Cone was used f o r data acquisition.  In t h i s i n v e s t i g a t i o n very clean signals were recorded with pre-  dominantly one S-wave frequency of 80 Hz present.  The calculated v e l o c i t i e s  from the Cross-Over method and CROSSCOR are given in Table B . 2 . 4 . t r a t e d , the v e l o c i t i e s obtained are in very close agreement. also plotted in Figure B.2.5 along  The v e l o c i t i e s are  with the cone bearing p r o f i l e .  65  As is i l l u s -  Table B.2.2  Calculated V e l o c i t y P r o f i l e s from McDonald Farm Comparing the Cross-Over Method and CROSSCOR. The Seismic Traces were Acquired with a 300 Hz Low-Pass F i l t e r Applied and Using SCPM with P i e z o e l e c t r i c Bender Accelerometers.  Cross-Over June 25, 1987 SCPM - CN Vs (m/sec)  CROSSCOR June 25, 1987 SCPM - CN Vs (m/sec)  Average Depth (m)  Depth (m) Interval  3.58  3.08-4.08  150  164  4.58  4.08-5.08  168  158  5.45  4.95-5.95  6.58  6.08-7.08  7.45  6.95-7.95  7.58  7.08-8.08  8.45  7.95-8.95  8.58  8.08-9.08  208  216  9.58  9.08-10.08  189  190  10.58  10.08-11.08  225  235  11.58  11.08-12.08  185  191  CROSSCOR June 25, 1987 SCMP - PM Vs (m/sec)  154 95  185 214  153  188 189  66  Velocity ( m/s ) 0  100  Figure B . 2 . 3 .  (t>)  200  Cone Bearing 0  100  Qt (bar) 200  (a) Ve.locity P r o f i l e from McDonald Farm Comparing the Cross-Over Method and CROSSCOR. The Seismic Traces were Acquired from an Accelerometer with a 300 Hz Analogue F i l t e r Applied, (b) Cone Bearing P r o f i l e from McDonald Farm.  67  Table B.2.3  Calculated V e l o c i t y P r o f i l e s from Laing Bridge, Comparing the Cross-Over Method with CROSSCOR. The Data was Obtained with an Accelerometer with a 300 Hz LowPass F i l t e r A p p l i e d . In t h i s Set of Data i t was not Possible to Obtain Accurate Cross-Over P i c k s .  Average Depth (m)  Seismic Cone Shear Wave Source Cross-Over Vs (m/sec)  Seismic Cone Shear Wave Source CROSSCOR Vs (m/sec)  Seismic Cone P-Plate Source CROSSCOR Vs (m/Sec)  4.3  94  140  5.3  144  120  6.3  124  171  7.3  154  167  8.3  143  138  9.3  188  187  10.3  152  197  189  11.3  210  150  152  12.3  145  144  13.3  120  197  14.3  193  183  15.3  197  224  16.3  246  68  Figure B . 2 . 4 .  V e l o c i t y P r o f i l e from Grant McConachie Way Comparing the Cross-Over Method with CROSSCOR. The Data was Obtained with an Accelerometer with a 300 Hz Low-Pass F i l t e r A p p l i e d . In t h i s Set of Data i t was not Possible to Obtain Accurate Cross-Over P i c k s , (b) Cone Bearing " P r o f i l e from the Laing Bridge"sTte.  69  Table B.2.4 Calculated V e l o c i t y P r o f i l e s from Lower Langley (232) Comparing the Cross-Over and CROSSCOR. The Seismic Traces were Acquired with an Accelerometer with a 300 Hz Analogue F i l t e r Applied. Cross-Over Hammer Source F i l t e r On Vs (m/sec)  CROSSCOR Hammer Source F i l t e r On Vs (m/sec)  Average Depth (m)  Depth (m) Interval  2.2  1.7-2.7  82  83  3.2  2.7-3.7  92  101  4.2  3.7-4.7  103  104  5.2  4.7-5.7  104  105  6.2  5.7-6.7  96  96  7.2  6.7-7.7  98  99  8.2  7.7-8.7  108  110  9.2  8.7-9.7  104  104  10.2  9.7-10.7  109  103  11.2  10.7-11.7  117  113  12.2  11.7-12.7  119  123  13.2  12.7-13.7  118  123  14.2  13.7-14.7  128  124  15.2  14.7-15.7  128  130  16.2  15.7-16.7  137  134  17.2  16.7-17.7  127  130  18.2  17.7-18.7  134  131  19.2  18.7-19.7 •  149  138  20.2  19.7-20.7  140  150  21.2  20.7-21.7  139  135  22.2  21.7-22.7  155  160  23.2  22.7-23.7  162  151  24.2  23.7-24.7  • 177  172  25.2  24.7-25.7  174  184  26.2  25.7-26.7  184  171  27.2  26.7-27.7  213  192  28.2  27.7-28.7  189  191  70  J  1  1  Figure B.2.5.  1  1  1  I  30  1—1  V e l o c i t y P r o f i l e from Lower Langley Comparing the CrossOver Method and CROSSCOR. The Seismic Traces were- Acquired with an Accelerometer with a 300 Hz Analogue F i l t e r A p p l i e d . (b) Cone Bearing P r o f i l e from the Lower Langley S i t e . "  71  B.3  Comparing CROSSCOR with F i l t e r e d (Analogue) and U n f i l t e r e d ( i . e . , Removing P-Waves and Noise With CROSSCOR) Seismic Traces In t h i s section CROSSCOR i s applied to u n f i l t e r e d seismic traces and com-  pared to traces f i l t e r e d with a 300 Hz low-pass analogue f i l t e r .  Figure B . 3 . 1  i l l u s t r a t e s a trace recorded at Laing Bridge in which the frequency spectrum c l e a r l y i l l u s t r a t e s three responses recorded by the accelerometer.  The  first  response i s the shear wave located at 100 Hz, the second response i s the compression wave at 800 Hz and the t h i r d corresponds to the accelerometer resonating at 3000 Hz.  From t h i s t r a c e , i t was possible to i s o l a t e the shear wave by apply-  ing a bandpass of 80 to 120 Hz, as i s i l l u s t r a t e d in Figure B . 3 . 2 .  The above  e x t r a c t i o n of the shear wave was possible with most of the u n f i l t e r e d seismic traces recorded and the s e l e c t i o n of the bandwidth was not a c r i t i c a l  criterion.  The important consideration was to make sure that the bandwidth was consistent throughout the p r o f i l i n g and contained a l l the frequencies w i t h i n the shear wave. In most of the data reduced, determination of the v e l o c i t y was r o u t i n e , because one only had to specify the same bandwidth throughout the p r o f i l e . Lower Langley 232nd S t . The f i r s t set of data i s from the Lower Langley S i t e .  In t h i s  investigation  the shear wave was located at 85 Hz, and the applied bandwidth was 70 to 100 Hz. The e f f e c t of time s h i f t s r e s u l t i n g from the 300 Hz low-pass analogue f i l t e r was determined by c r o s s c o r r e l a t i n g the f i l t e r e d and u n f i l t e r e d t r a c e s . procedure, i t was found that time s h i f t s were n e g l i g i b l e . c i t i e s are i l l u s t r a t e d in Table B . 3 . 1 .  Figure B.3.3(a)  From t h i s  The c a l c u l a t e d v e l o -  i s a plot of the v e l o c i t i e s  c a l c u l a t e d from CROSSCOR and the reverse p o l a r i t y method, and the cone bearing p r o f i l e i s given in Figure B . 3 . 3 ( b ) .  The Figure i l l u s t r a t e s that the u n f i l t e r e d  and f i l t e r e d r e s u l t s are nearly i d e n t i c a l and vary along with the cone bearing profile.  Thus, in t h i s set of data CROSSCOR separated out the shear waves from  the P-waves and n o i s e , and gave v e l o c i t i e s which were repeatable and varied as the cone bearing p r o f i l e . T i l b u r y Island Natural Gas Plant S i t e The second set of data i s from the T i l b u r y Natural Gas Plant S i t e , where i t was not possible to pick Cross-Overs.  Therefore,  invaluable tool in obtaining the shear v e l o c i t i e s .  72  CROSSCOR was found to be a In t h i s case the shear wave  3.8 Metre Seismic Trace  1.0  (3_  -1.0  0.0 (a)  1.0  TIME  (msec)  Seismic Trace Recorded at Laing Bridge with the P-Plate Source ( U n f i l t e r e d ) .  S (100 Hz)  j  M =3  i—• !3 •CC  P (800 Hz)  \L J ,  |  3000  FREQUENCY HERTZ  (b)  5000  Frequency Spectrum of the Above Tr.ace I l l u s t r a t i n g Three Responses Recorded by the Accelerometer: S-Wave at 100 Hz, P-Wave at 800 Hz, and Resonanting Accelerometer at 3000 Hz. .  Figure B . 3 . 1 .  Seismic Trace from Laing Bridge  73  3.8  metre  •1.0 3.8  TIME  metre  6.8  msecj (msec) FILTERED  '  (80-120 Hz) ' 1.01  metre  TIME  6.8 metre  (msec) FILTERED  (80-120 Hz) 1  -1.0 TIME  (msec)  TIME  (msec)  (a) Seismic Traces from Laing Bridge from Depths of 3.8 and 6.8 Metres. These Wavelets are then F i l t e r e d with a Bandpass of i 80 to 120 Hz Inorder to Extract the Shear Wave.  CROSSCORRELATION  0.0  TIME SHIFT  200 (msec)  (b) Corresponding Crosscorrelation Function  Figure B . 3 . 2 .  Seismic Traces from Laing Bridge  74  Table B . 3 . 1 Calculated Velocity P r o f i l e s from Lower Langley (232) Comparing F i l t e r e d (300 Hz Low-Pass) and U n f i l t e r e d Seismic Data. The Data was Acquired with an Accelerometer.  Average Depth (m)  Depth (m) Interval  Cross-Over Hammer Source F i l t e r On Vs (m/sec)  2.2  1.7-2.7  82  3.2  2.7-3.7  4.2  3.7-4.7  Cross-Over Hammer Source F i l t e r Off Vs (m/sec)  CROSSCOR Hammer Source F i l t e r On Vs (m/sec)  CROSSCOR Hammer Source F i l t e r Off Vs (m/sec)  92  * *  83 101  101  103  *  104  109  105  99  5.2  4.7-5.7  104  *  6.2  5.7-6.7  96  *  96  96  99  100  7.2  6.7-7.7  98  *  8.2  7.7-8.7  108  *  110  110  104  111  9.2  8.7-9.7  104  *  10.2  9.7-10.7  109  *  103  104  113  120  11.2  10.7-11.7  117  *  12.2  11.7-12.7  119  *  123  120  13.2  12.7-13.7  118  *  123  120  124  127  14.2  13.7-14.7  128  *  15.2  14.7-15.7  128  *  130  133  134  138  16.2  15.7-16.7  137  *  17.2  16.7-17.7  127  *  130  124  131  141  18.2  17.7-18.7  134  *  19.2  18.7-19.7  149  *  138  134  20.2  19.7-20.7  140  *  150  160  135  131  21.2  20.7-21.7  139  *  22.2  21.7-22.7  155  *  160  160  151  155  23.2  22.7-23.7  162  *  24.2  23.7-24.7  177  *  172  172  25.2  24.7-25.7  174  *  184  185  26.2  25.7-26.7  184  *  171  168  192  199  191  198  27.2  26.7-27.7  213  *  28.2  27.7-28.7  189  *  * => not possible to pick Cross-Overs because of varying s i g n a l s . 75  Figure B . 3 . 3 .  (a) Velocity P r o f i l e from Lower Langley Comparing F i l t e r e d . (300 Hz Low-Pass) and Unfiltered Seismic Data. The Data was Acquired w i t h i n Accelerometer. (b) Cone Bearing P r o f i l e from Lower;Langley.  76  maximum energies were centered at 100 Hz and the r e s u l t i n g v e l o c i t i e s are i l l u s trated in Table B.3.2 and compared in Figure B.3.4 with the cone bearing p r o f i l e . The calculated v e l o c i t i e s vary with the cone bearing p r o f i l e except f o r the v e l o c i t i e s at 3.66 and 6.64 metres.  This maybe due to poor responses being  recorded at these depths.  77  Table B.3.2  Calculated V e l o c i t y P r o f i l e s from Tilbury S i t e Comparing Shear Wave V e l o c i t i e s Obtained with CROSSCOR from Seismic Data Acquired by Two Different Accelerometers. In t h i s Set of Data i t was not Possible to Obtain Accurate CrossOver P i c k s .  Average Depth (m)  SCPM-CN Seismic Cone Accelerometer Vs (m/sec)  SCPM-PM Pressuremeter Accelerometer Vs (m/sec)  2.37  57  2.91  126  3.22  113  3.66  91  4.72  61  4.94  131  5.39  93  5.45  134  6.11  121  6.12  123  6.64  88  7.11  146  8.11  172  78  (a)  Velocity ( m/s ) 0 i  i  i  80 I  i  i  i  Cone Bearing Qc ( bar )  160 I  ***** SCPM CN ooooo SCPM PM  o *  o  CL  J  Figure B . 3 . 4 .  L  Calculated Velocity P r o f i l e s from T i l b u r y Site Comparing Shear Wave V e l o c i t i e s Obtained with CROSSCOR from Seismic Data Acquired by Two Different Accelerometers. The V e l o c i t y P r o f i l e s are also I l l u s t r a t e d with the Cone Bearing P r o f i l e . In t h i s Seismic Data i t was not Possible to Obtain Accurate Cross-Over P i c k s , (b) Cone Bearing P r o f i l e from the T i l b u r y Island S i t e .  79  B.4  CROSSCOR with Nonpolarized Sources In this section the a p p l i c a t i o n of CROSSCOR to nonpolarized shear sources  is addressed.  The two sources discussed are the Buffalo Gun and comparing one  side of the Hammer Shear Source to the other side ( i . e . , right and l e f t ) .  The  Buffalo Gun source i s thoroughly discussed in Reference 2 by Laing (1985) with regard to design and use.  The Buffalo Gun i s desirable because i t can be used  both on and o f f shore to obtain both shear and compression wave v e l o c i t i e s w i t h out the problems associated with a seismic cap or explosive source.  CROSSCOR  i s applied to data obtained from s t r i k i n g each Hammer Shear Source in one d i r e c t i o n , therefore, i t allows one to obtain two v e l o c i t y estimates for one depth increment (one for h i t t i n g the l e f t side and one for h i t t i n g the r i g h t s i d e ) . This i s opposed to the reverse p o l a r i t y method which allows for only one v e l o c i t y estimate for each depth increment. Lower Langley 232nd St. The f i r s t set of data discussed i s from the Lower Langley S i t e , where an extensive Buffalo Gun p r o f i l e was obtained.  This seismic section i s i l l u s t r a t e d  in Figure B . 4 . 1 .  From t h i s f i g u r e , one can c l e a r l y d i s t i n g u i s h the compression  and shear waves.  The shear waves.were obtained by bandpassing the data from 90  to 150 Hz.  The r e s u l t i n g shear wave v e l o c i t i e s are given in Table B.4.1 where  they are compared to the Hammer Shear Source v e l o c i t i e s . are comparable, while others deviate considerably.  Some of these values  This deviation was found to  r e s u l t from t r i g g e r i n g problems with the Buffalo Gun.  The Buffalo Gun source i s  triggered with a geophone which may r e s u l t in inconsistency.  This i s i l l u s t r a t e d  in Figure B . 4 . 2 , where two seismic traces from the same depth (13 metres) were crosscorrelated and a r e s u l t i n g 2.71 msec s h i f t was recorded.  For a Vs of 100  m/sec the i n t e r v a l time i s 10 msec over a 1 metre distance and 2.71 msec error represents +27% error on measured Vs.  If the t r i g g e r error were a constant 2.72  msec then the potential maximum e r r o r increases as Vs increases.  The measured  t r i g g e r error of 2.72 msec appears to be a maximum when compared to the percent differences between Hammer and Buffalo Gun data in Table B . 4 . 1 .  This t r i g g e r i n g  problem w i l l be resolved by using an accelerometer (having a fast response time) instead of the geophone as well as gounding the Buffalo Gun f i r i n g pin to create a contact closure ( s i m i l a r to the shear p l a t e ) .  80  0  20  40  60  80  TIME  Figure B . 4 . 1 .  100  120  140  {m sec)  Seismic Section from Lower Langley Where a Buffalo Gun was Used as the Source and An U n f i l t e r e d Accelerometer Recorded the Data. 81  cd C  C/3 i-l O  >  00  X) V •1-1  1-1  0.  low-pass analogue f i l t e r applied  0 Figure B . 4 . 2 .  20  40  60' TIME  80  100  120'  140  (msec).  Two Buffalo Gun Seismic Traces from Lower Langley I l l u s t r a t i n g t r i g g e r i n g Prop ems with the Buffalo Gun Source. The Data was A c o u i r e T ^ T r i r t r , Accelerometer at a Depth of 13 Metres. The Resulting Time W f f e r i H F e was 1  Table B . 4 . 1  Calculated V e l o c i t y P r o f i l e s of Lower Langley (232) Comparing Shear Wave V e l o c i t i e s Obtained from the Hammer Shear Source and the Buffalo Gun Source. The Data was Acquired with an Accelerometer ( U n f i l t e r e d ) .  Depth (m) Interval  CROSSCOR Hammer Source F i l t e r Off Vs (m/sec)  CROSSCOR Buffalo Gun F i l t e r Off Vs (m/sec)  % Difference (Hammer Reference)  2.7-3.7  101  106  +5  3.7-4.7  109  94  -14  4.7-5.7  99  90  -9  5.7-6.7  96  96  0  10.7-11.7  120  114  -5  11.7-12.7  120  120  -0  12.7-13.7  120  101-low  -16  13.7-14.7  127  108-low  -15  14.7-15.7  133  155-high  + 17  15.7-16.7  138  133-low  -18  16.7-17.7  124  122  - 2  17.7-18.7  141  165-high  +17  18.7-19.7  134  114-1ow  + 15  Note:  One t r i g g e r error creates two incorrect successive measurements of Vs.  83  Annacis Bridge P i l e Research S i t e The second set of data is from the Annacis P i l e Research s i t e , where v e l o c i t i e s obtained from s t r i k i n g each side of the pads for the Hammer Shear Source are compared.  The v e l o c i t i e s from the SPCM-CN are given in Table B.4.2 and  plotted along with the cone bearing p r o f i l e in Figure B . 4 . 3 .  The low v e l o c i t i e s  are expected for t h i s very soft organic s i l t and tend to l i n e a r increase with depth as the i n - s i t u stresses increase ( i . e . , e s s e n t i a l l y constant cone bearing p r o f i l e which has been corrected for overburden s t r e s s e s ) .  The higher  velocity  values at 12.0 metres are i n d i c a t i v e of the variable (in depth) sand unit at approximately 15.0 metres.  The v e l o c i t i e s obtained from the UBC #8 cone are  given in Table B.4.3 and plotted in Figure B.4.4 along with the cone bearing profile.  A special note must be made with respect to the cone bearing p r o f i l e s  shown in Figures B . 4 . 3 and B . 4 . 4 .  These two p r o f i l e s are the same and are from  a d i f f e r e n t i n v e s t i g a t i o n ( i . e . , d i f f e r e n t CPT locations) than the v e l o c i t i e s shown in Figures B . 4 . 3 and B . 4 . 4 .  Thus, the p r o f i l e s are only a rough guide for  the change in s o i l p r o f i l e with depth ( i . e . , one should expect some v a r i a b i l i t y with d i f f e r e n t CPT l o c a t i o n s ) .  The v e l o c i t i e s shown for the Annacis P i l e Re-  search S i t e are very comparable in each set of data and give greater confidence in the r e s u l t s .  84  Table B.4.2  Calculated V e l o c i t y P r o f i l e s for Annacis P i l e S i t e , where V e l o c i t i e s Obtained from S t r i k i n g each Side of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer ( U n f i l t e r e d ) .  Depth (m)  January 8 , 1988 Pressuremeter Right Side Vs (m/sec)  January 8 , 1988 Pressuremeter Left Side Vs (m/sec)  1.97-2.97  55  58  +5  2.97-3.97  41  39  -5  3.97-4.97  50  52  +4  4.97-5.97  63  64  +2  5.97-6.97  84  105  +25  7.97-8.97  90  100  +11  8.97-9.97  64  54  -16  9.97-10.97  58  56  -3  10.97-11.97  152  138  -14  85  % Difference (Reference-Riqht Side)  Velocity ( m/s )  Cone Bearing Qt (bar)  100  * * * * * CROSSCOR S C P V - P U  (b)  Rtciht S M «  acsaB CROSSCOR S C P M - P M un  Figure B . 4 . 3 .  0  SJ6«  (a) Calculated V e l o c i t y P r o f i l e s from the Annacis P i l e Research S i t e , where V e l o c i t i e s Obtained from both Sides of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer ( U n f i l t e r e d ) . (b) Cone Bearing P r o f i l e of the Annacis Bridge P i l e Research S i t e  86  Table B.4.3  Calculated V e l o c i t y P r o f i l e s for Annacis P i l e S i t e , where V e l o c i t i e s Obtained from S t r i k i n g each Side of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer ( U n f i l t e r e d ) .  Depth (m)  CROSSCOR January 8 , 1988 Cone #8 Right Scale Vs (m/sec)  CROSSCOR January 8 , 1988 Cone #8 Left Scale Vs (m/sec)  % Difference (Reference-Right Side)  1.84-2.84  19  24  +25  2.84-3.84  41  41  0  3.84-4.84  22  21  -4  4.84-5.84  75  85  +13  5.84-6.84  82  82  0  6.84-7.84  91  87  -4  7.84-8.84  114  122  +7  8.84-9.84  94  97  +3  9.84-10.84  123  128  +4  10.84-11.84  137  152  +11  87  Figure B . 4 . 4 .  (a) Calculated Velocity P r o f i l e s from the Annacis P i l e Research S i t e , where V e l o c i t i e s Obtained from both Sides of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer ( U n f i l t e r e d ) . (b) Cone Bearing P r o f i l e of the Annacis Bridge P i l e Research S i t e  88  B.5  CROSSCOR in Obtaining Compression Wave V e l o c i t i e s In t h i s section the a p p l i c a t i o n of CROSSCOR to obtaining compression wave  v e l o c i t i e s i s addressed.  As was discussed previously,  the seismic traces  obtained from u n f i l t e r e d accelerometer data contain information of both shear and compression body waves. on sources used.  The extent of t h i s information l a r g e l y depends  A point source such as a Buffalo Gun would give more i d e n t i -  f i a b l e shear and compression waves, as opposed to a shear source ( i . e . , Hammer Shear Source) which gives predominantly shear waves.  In order to i s o l a t e the  i n d i v i d u a l body wave responses, i t i s necessary to obtain the frequency spectrum of the seismic traces recorded and bandpass the desired frequencies.  In  Figure B . 5 . 1 a recorded wavelet i s i l l u s t r a t e d along with i t s frequency spectrum. This data was obtained with a P - P l a t e source at Laing Bridge S i t e at a depth of 6.8 metres. responses.  From the frequency spectrum obtained, one can i d e n t i f y three dominant The f i r s t response i s the shear wave occurring at 150 Hz, the second  response i s the compression wave at 900 Hz, and the t h i r d response resulted from the resonanting accelerometer at 3000 Hz.  In order to i s o l a t e the compres-  sion wave, a Bandpass ( e . g . , 800 to 1200 Hz) must be applied which focuses on 900 Hz. The body wave v e l o c i t i e s are l a r g e l y a function of the i n - s i t u stresses and the composition of the sediments.  Figure B.5.2 i l l u s t r a t e s a plot of shear and  compression wave v e l o c i t i e s versus density for a l l types of sediments and rocks. As i s i l l u s t r a t e d in t h i s f i g u r e , compression waves can vary from 1100 m/sec to 4000 m/sec in sediments.  In saturated materials (no gas present) the compression  wave v e l o c i t y must be at l e a s t greater than that for water or about 1500 m/sec. Table B . 5 . 1 shows that t y p i c a l compression wave v e l o c i t i e s for clay range from 1100 to 2500 m/sec, and those for loose sand are predominantly 1800 m/sec i f saturated. Lower Langley 232nd S t . The f i r s t set of data discussed i s from Lower Langley where a Buffalo Gun source was used.  As was i l l u s t r a t e d in Figure B . 4 . 1 , the seismic traces obtained  from Buffalo Gun sources c l e a r l y contain i d e n t i f i c a b l e S and P waves.  The present  problem with using t h i s data (as was shown in Section B.4) was that t r i g g e r i n g was i n c o n s i s t e n t .  The compression wave v e l o c i t i e s obtained from the B u f f a l o Gun  data are shown i n Table B.5.2 where a bandpass of 200 to 1000 Hz was a p p l i e d . The nonsensical values obtained are most l i k e l y a r e s u l t from t r i g g e r i n g problems. 89  INPUT HAUELET  o.o (a)  TIME  (msec)  Seismic Trace from Laing Bridge where a P-Plate Source was Used and an Accelerometer ( U n f i l t e r e d ) Recorded the Data at 6.8 Metres.  FREQUENCY SPECTRUM OF TRACE  1.0  •-3  3000  FREQUENCY HERTZ (b)  5000  Corresponding Frequency Spectrum of Seismic Trace I l l u s t r a t i n g the Three Dominant Responses: Shear Wave at 150 Hz, Compression Wave at 900 Hz, and Resonanting Accelerometer at 3000 Hz. Figure B . 5 . 1 .  Seismic Trace from Laing Bridge  90  Figure B . 5 . 2 .  V e l o c i t y Versus Density f o r Compressional and Shear Waves in a l l Types of Sediments and Rocks. Poisson's Ratio Versus Density at Top of Figure. (Milton B. Dobrin, 1976).  91'  Table B . 5 . 1 Compressional and Shear V e l o c i t i e s in Sediments and Rocks (Milton B. Dobrin, 1976)  Compressional velocity  Shear velocity  Material and Source  m/s  ft/s  Granite: Barrieficld, Ontario Quincy, Mass. Bear Mt., Tex.  5640 5880 5520  18,600 19,400 17,200  2870 2940 3040  9470 9700 10,000  Granodioritc, Weston, Mass.  4780  15,800  3100  10,200  Diorite, Salem, Mass.  5780  19,100  3060  10,100  Gabbro, Duluth, Minn.  6450  21,300  3420  11,200  Basalt, Germany  6400  21,100  3200  10,500  Dunite: Jackson City, N.C. Twin Sisters, Wash.  7400 8000  24,400 28,400  3790 4370  12,500 14,400  2880 3030  9500 10,000  500  1,650  Sandstone Sandstone conglomerate, Australia  1400-1300 4620-14,200 2400 7920  Limestone: Soft Solenhofcn, Bavaria Argillaceous, Tex. Rundle, AJberta  1700-1200 5610-13,900 5970 19,700 6030 19,900 6060 20,000  Anhydrite, U.S. Midcontinent, Gulf Coast Clay Loose sand  Saturated Saturated  SOURCE: Sydney V. Mem. 97, 1966.  4100  13,530  1100-2500 1800  3630-8250 5940  m/s  ft/s  Clark, Jr. (ed.), "Handbook of Physical Constants," rev. ed., Geol.  92  Soc. Am.  Table B.5.2  Calculated V e l o c i t y P r o f i l e from Lower Langley, where an Accelerometer was Used for Data A c q u i s i t i o n . The Compression Wave V e l o c i t i e s were Obtained from the Buffalo Gun Source and by Applying a Bandpass of 200 to 1000 Hz.  CROSSCOR Buffalo Gun F i l t e r Off Vp (m/sec)  Depth (m)  * =>'  1.7-2.7  243  2.7-3.7  1200  3.7-4.7  1443  4.7-5.7  330  5.7-6.7  353  10.7-11.7  330  11.7-12.7  332  12.7-13.7  *  13.7-14.7  *  14.7-15.7  760  15.7-16.7  3051  16.7-17.7  3060  17.7-18.7  383  nonsensical compression wave v e l o c i t i e s due to poor and/or uncorrelated responses.  93  I l l u s t r a t e d in Figure B . 5 . 3 . are traces obtained from the Buffalo Gun data along with the i s o l a t e d compression wave responses.  This figure i l l u s t r a t e s that the  compression wave responses are c l e a r l y i d e n t i f i a b l e and correlated ( i . e . , lar).  simi-  Thus, one can conclude that i f t r i g g e r i n g was repeatable, accurate com-  pression wave v e l o c i t i e s could be obtained. The next set of data was from Lower Langley where the Hammer Shear was used. The compression waves v e l o c i t i e s obtained from t h i s set of data  (unfiltered  accelerometer) are shown in Table B.5.3 and were obtained with a bandpass of 750 to 1000 Hz a p p l i e d .  The nonsensical v e l o c i t i e s are most l i k e l y due to the f a c t  that cl e a r compression wave responses were not obtained.  Figure B.5.4 i l l u s -  trates that the i s o l a t e d compression wave responses are very noisy and poorly correlated.  94  4.7  0.0  TIME 4.7  5.7  metres  metres  TIME  (msec)  metres (200-10Q0 Hz)  5.7  1.0  (msec)  metres (200-1000 Hz)  CD 3  CL  E  -1.0  TIME 11.7  metres  TIME 11.7  1.0  TIME  (msec)  12.7  0,0  ( msec)  metres (200-1000 Hz)  1 0 0  metres  TIME 12.7  (msec)  (msec)  1 0 0  metres (200-1000 Hz)  j P-Wave  CD  CL  E <:  •1.0 0.0  Figure B . 5 . 3 .  T T UV  1 IME(msec)  TIME  1 0 0  (msec)  1 0 0  Seismic Traces Obtained at Lower Langley where a Buffalo Gun Source was Used. The Traces were F i l t e r e d (200 to 1000 H z ) i n order to Extract the Compression Wave Responses The r e s u l t i n g Compression Wave Responses are Highly Correlated. g 5  Table B.5.3 Calculated Velocity P r o f i l e s from Lower Langley, where an Accelerometer was Used for Data A c q u i s i t i o n . The Compression Wave V e l o c i t i e s were Obtained by Applying a Bandpass of 750 to 1500 Hz.  CROSSCOR Hammer Source F i l t e r Off Vp (m/sec)  Depth (m) 2.7-3.7  1042  3.7-4.7  558  4.7-5.7  579  5.7-6.7  676  6.7-7.7  1600  7.7-8.7  692  8.7-9.7  542  9.7-1Q.7 10.7-11.7  2448 *  11.7-12.7  *  12.7-13.7  4937 *  13.7-14.7 14.7-15.7 15.7-16.7  1651 *  16.7-17.7  *  17.7-18.7  *  18.7-19.7  *  19.7-20.7  *  20.7-21.7  *  21.7-22.7 22.7-23.7  1244 *  23.7-24.7  *  24.7-25.7  *  25.7-26.7 26.7-27.7  1246 *  27.7-28.7  *  * => nonsensical compression wave v e l o c i t i e s are due to poor and/or uncorrelated responses . Q(  9.7  metres  1.0 cu -o 3 +->  -1.0 Tint  8.7  1.0  0.0  (msec)  metres (750-1000 Hz)  TIKE  g  7  m  e  t  r  e  (msec)  s  ( 7 5 0  _ i 0 0 Hz) 0  cu -o 3  •o  -t->  3  -t->  CL  QE  E  -1.0  °-°  TIKE  14.7  Figure B . 5 . 4 .  (msec)  metres  200  -1.0  °-°  TIKE  (msec)  200  15.7 metres  Seismic Traces Obtained at Lower Langley Where a Hammer Shear Source was Used. The Traces were F i l t e r e d (1750 to 2000 Hz) in Order to Extract the Compression Wave Responses. The R e s u l t ing Compression Waves are Poorly Correlated. 97  B.6  Bandwidth Considerations when Applying CROSSCOR In Section IIIA, i t was stated that the crosscorrelation function provides  information about the s i m i l a r i t y and time difference between two successive wavelets.  In t h i s case we are looking f o r the a r r i v a l time difference between  the dominant responses of the recorded wavelets.  The main idea behind CROSSCOR  is to i s o l a t e s p e c i f i c body waves in the frequency domain and then crosscorrel a t e the f i l t e r e d wavelets in the time domain.  The shear and compression body  waves are i d e n t i f i e d by t h e i r own band of frequencies where a dominant frequency characterizes the respective wavelet.  From the extensive amount of seismic data  processed in this research, i t was found that in specifying the bandwidth in f i l t e r i n g the most important considerations are t o : (1)  Include dominant frequency characterizing wavelet, because the other contributing frequencies do not e f f e c t the dominant responses s i g n i f i cantly.  (2)  The bandwidths s p e c i f i e d must be consistent throughout the seismic p r o f i l e in order that the crosscorrelated dominant responses are s i m i l a r in transient response s p e c i f i c a t i o n s ( e . g . , page 1.4)  In Figure B.6.1 there i s a seismic trace (at 2.7 metres) from Lower Langley i l l u s t r a t e d along with i t s ' frequency spectrum.  From t h i s frequency spectrum  the dominant shear wave response i s i d e n t i f i e d at 50 Hz.  Figure B.6.2 shows a  seismic trace recorded at a depth of 3.7 from the same seismic cone p r o f i l e . From the frequency spectrum of t h i s seismic wavelet, the dominant shear wave r e s ponse i s also i d e n t i f i e d at 50 Hz.  In Figure B.6.3 the two above wavelets are  i l l u s t r a t e d along with t h e i r f i l t e r e d responses.  The f i r s t f i l t e r e d set of data  was bandpassed from 20 to 120 Hz and the second set of f i l t e r e d data was bandpassed from 30 to 70 Hz.  These bandwidths were chosen in order that the f i l -  tered wavelets contain the dominant shear wave response of 50 Hz.  When cross-  c o r r e l a t e d , the two d i f f e r e n t sets of f i l t e r e d data gave the same time s h i f t and corresponding v e l o c i t y of 102 m/sec. Figure B.6.4 i l l u s t r a t e s a seismic trace recorded at a depth of 3.7 metres at the Lower Langley S i t e .  The frequency spectrum i l l u s t r a t e d here i s s l i g h t l y  more complex than those shown in Figures B.6.1 and B . 6 . 2 , because a Buffalo Gun source was used as opposed to the Hammer Shear Source. dominant shear wave responses are from 50 to 200 Hz. 98  From Figure B.6.4 the Figure B.6.5 i l l u s t r a t e s  2.7  metres  FREQUENCY HERTZ (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response at 50 Hz. Figure B . 6 . 1 .  Seismic Trace from Lower Langley Where a Hammer Shear Source Was Used. 99  3.7  Metres  0.0  200  TIME (a)  ( msec)  Seismic Trace at 3.7 Metres Acquired with An U n f i l t e r e d Accelerometer.  50  500  FREQUENCY HERTZ (b)  Figure B . 6 . 2 .  Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response at 50 Hz.  Seismic Trace Recorded at Lower Langley Where a Hammer Shear Source was Used.  100  2.7 metres  •0  3.7 metres  TTUP i TIME (msec)\  200  0.0  2.7 Metres (20-200 Hz)  TTiitv \ TIME(msec)  200  3.7 Metres (20-200 Hz)  -1.0 °-°  Figure B . 6 . 3 .  Tiyp  /  \  l i n t (msec)  200  0.0  TTUP/'  ^  TIME (msec)  200  (a) Seismic Traces Recorded at Depths of 2.7 and 3.7 Metres, (b) Trace in (a) F i l t e r e d Between 20 to 200 Hz Resulting in a V e l o c i t y of 102 F i l t e r e d Between 30 to 70 Hz Resulting in a Velocity ' o f 102 m/sec.  101  3.7  Metres  1.0  TIME (a)  (  m s e  c)  Seismic Trace at 3.7 Metres Acquired with an U n f i l t e r e d Accelerometer.  S-Wave  50  200  300  400  500  FREQUENCE HERTZ (b)  Figure B . 6 . 4 .  Frequency Spectrum of ( a ) , I l l u s t r a t i n g Dominant Shear Wave Responses Situated Between 50 to 200 Hz. Seismic Trace Recorded at Lower Langley Where a B u f f a l o Gun Source was Used 102  4.7  Metres  0.0  TIME (a)  (msec)  Seismic Trace at 4.7 Metres Acquired with an U n f i l t e r e d Accelerometer.  S-Wave  200  500  FREQUENCY HERTZ (b)  Figure B . 6 . 5 .  Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response Situated Between 50 to 200 Hz.  Seismic Trace Recorded at Lower Langley where a Buffalo Gun Source was Used  103  another seismic trace at a depth of 4.7 metres from the same seismic cone p r o f i l e as the trace in Figure B . 6 . 4 .  From the frequency spectrum i l l u s t r a t e d in Figure  B . 6 . 5 , the dominant shear wave responses are situated from 50 to 200 Hz.  Figure  B.6.6 i l l u s t r a t e s the above two Buffalo Gun traces along with t h e i r f i l t e r e d responses.  The f i r s t set of f i l t e r e d data was bandpassed from 20 to 200 Hz,  and the second set of f i l t e r e d data was bandpassed from 60 to 180 Hz.  Cross-  c o r r e l a t i n g these two d i f f e r e n t sets of f i l t e r e d data resulted in the same time s h i f t and v e l o c i t y of 94 m/sec. bandpassed between 100 to 200 Hz.  Figure B.6.7 shows the same Buffalo Gun Traces The r e s u l t i n g v e l o c i t y (93 m/sec) obtained  from t h i s set of data was n e g l i g i b l y d i f f e r e n t from those previously obtained.  104  3.7  1.0  Metres  4.7  Metres  1.0  cu -a 13  m •a  •4->  4->  1  CL  CL  E  E  -1.0  -1.0 0.0  TIME  3.7  1.0  160  (msec)  0.0  Metres (60-180 Hz)  TIME  4.7  (msec)  160  Metres (60-180 Hz)  cu -a CL  E  •1.0  °-° Figure B . 6 . 6 .  TIME  (msec)  160  0.0  TIME  (msec)  (a) Seismic Traces Recorded at Depths of 3.7 and 4.7 Metres, (b) Trace in (a) F i l t e r e d Between 20 to 200 Hz Resulting in a V e l o c i t y of 94 m/sec. (c) Traces in (a) F i l t e r e d Between 60 to 180 Hz Resulting in a V e l o c i t y of 94 m/sec.  105  4.7 metres  3.7 metres (100-200 Hz)  °'°  TIME  Figure B . 6 . 7 .  (msec)  0-0'  1 6 0  T  lint T  M  r  (100-200 Hz)  , ^ 160 (msec)  Trace in Figure B.6.6(a) F i l t e r e d Between 100 to 200 Hz Resulting in a V e l o c i t y of 93 m/sec.  106  B.7  CONCLUSIONS  The r e s u l t s from the previous discussions suggest that CROSSCOR i s a very b e n e f i c i a l tool to obtain seismic wave v e l o c i t i e s .  In these i n v e s t i -  g a t i o n s , CROSSCOR was also found to be very useful to reduce the noisy accelerometer data where picking cross-overs was d i f f i c u l t or impossible. The program allows one to separate but d i f f e r e n t seismic events and use a l l the information contained in the wavelets in obtaining shear and compression wave v e l o c i t i e s .  When conducting the data r e d u c t i o n s , i t was found that the  d i f f e r e n t wavelet signatures were consistent throughout a cone p r o f i l e .  Thus,  separating out the shear and compression wave instrument responses was routine and CROSSCOR proved to be user f r i e n d l y .  In a d d i t i o n , since CROSSCOR i s a  d i g i t a l f i l t e r i t removes most human bias in determining a r r i v a l times, the program allows one to compare v e l o c i t i e s from s i g n a l s obtained on opposite sides of the source, and one may also have the f l e x i b i l i t y in using d i f f e r e n t sources in obtaining seismic wave v e l o c i t i e s . In applying CROSSCOR to Buffalo Gun Source data i t was found that compression and shear waves could be i d e n t i f i e d and i s o l a t e d .  The problem with  processing Buffalo Gun data was that t r i g g e r i n g was not consistent out a seismic p r o f i l e r e s u l t i n g in erroneous v e l o c i t y c a l c u l a t i o n s .  thorough It was  suggested that the t r i g g e r i n g problem could be resolved by using an a c c e l e r o meter (having a f a s t response time) instead of the geophone as well as gounding the Buffalo Gun f i r i n g pin to create a contact closure ( s i m i l a r to the shear p l a t e ) . The performance of CROSSCOR i s summarized in the f o l l o w i n g l i s t of advantages/disadvantages with comparisons made to the reverse p o l a r i t y method in the l i s t of advantages and comparisons made to the following Kalman F i l t e r formulation in the l i s t of disadvantages. Advantages (1)  CROSSCOR i s able to work with noisy accelerometer data e x t r a c t i n g shear and compression waves.  (2)  The c r o s s c o r r e l a t i o n function in CROSSCOR makes use of a l l the information in the s i g n a l s (averaging out i r r e g u l a r i t i e s and putting s i g n i f i c a n c e or dominant responses) as opposed to the reverse p o l a r i t y method which r e l i e s only on one cross-over point to determine i n t e r v a l times. 107  (3)  CROSSCOR works as well as the reverse p o l a r i t y method when reducing clean s i g n a l s , and CROSSCOR works better than the reverse  polarity  method with s i g n a l s containing many dominant low frequencies  (e.g.,  Figures 2.8 and B . 2 . 1 ) . (4)  CROSSCOR can obtain v e l o c i t y estimates from nonpolarized sources  (e.g.,  the Buffalo Gun Source and comparing one side of the Hammer Shear Source to the other).  Thus, CROSSCOR has greater f l e x i b i l i t y than the reverse  p o l a r i t y method in working with d i f f e r e n t types of sources and g i v i n g more v e l o c i t y estimates f o r a downhole seismic p r o f i l e as opposed to the reverse p o l a r i t y method. (5)  CROSSCOR can extract compression wave events from seismic traces and give corresponding compression wave v e l o c i t i e s .  (6)  CROSSCOR i s a user f r i e n d l y program which has f l e x i b i l i t y in s p e c i f y i n g  the cutoff frequencies i f there i s a repeatable t r i g g e r of f a s t response and s u f f i c i e n t accuracy and r e s o l u t i o n . Pi sadvantages (1) CROSSCOR cannot give s o i l damping p r o p e r t i e s . (2)  CROSSCOR cannot give any i n d i c a t i o n of the shape of the wavelets being recorded (which i s a very important parameter in earthquake engineering).  108  IV.  KALMAN FILTER FOR ESTIMATING AMPLITUDES & ARRIVAL TIMES IN P AND S WAVE ANALYSIS - A PRELIMINARY ASSESSMENT This section gives a b r i e f description of the Kalman F i l t e r , KF, and  presents the equations for using the KF to estimate P-wave and S-wave amplitudes and a r r i v a l times.  More d e t a i l s of the KF equations, with  examples and applications to geophysical problems, can be found in Reference 24. The KF should be investigated for i t s possible a p p l i c a t i o n to obtaining accurate estimates of the P-wave and S-wave amplitudes and a r r i v a l times, for the p a r t i c u l a r problem we are considering, due to both the f l e x i b i l i t y of the KF approach and to l i m i t a t i o n s of the frequency domain approaches that were discussed in Section  III.  The Kalman F i l t e r i s an optimal (in a least squares sense) f i l t e r which is based on state-space, time-domain formulations of physical problems. Application of t h i s f i l t e r requires that the physical problem be modelled by a set of f i r s t order d i f f e r e n t i a l or difference equations which, with i n i t i a l conditions, uniquely define the system behaviour.  The f i l t e r u t i l i z e s know-  ledge of system and measurement dynamics, assumed s t a t i s t i c s of system noises and measurement e r r o r s , and s t a t i s t i c a l information about the i n i t i a l c o n d i tions.  Figure 4.1 i l l u s t r a t e s the basic r e l a t i o n between the system, the  measurements, and the Kalman F i l t e r , for a p p l i c a t i o n to physical problems. The Kalman F i l t e r takes into account the s t a t i s t i c s of the measurement and state e r r o r s , and the a p r i o r i model information provides for optimal use of any number, combination, and sequence of external measurements.  The Kalman  F i l t e r can be applied to problems with l i n e a r time-varying systems, and with non-stationary system and measurement s t a t i s t i c s .  Problems with n o n l i n e a r i -  t i e s can sometimes be handled by l i n e a r i z i n g the system and measurement equations.  Furthermore, the Kalman F i l t e r is r e a d i l y applied to estimation,  smoothing, and p r e d i c t i o n . One of the most important aspects that i s r e a l i z e d when applying the KF to solving a s p e c i f i c problem, is that a considerable e f f o r t is required to model the physical problem to f i t into the framework of the KF.  The physical  problem must be modelled/approximated by a set of f i r s t order d i f f e r e n t i a l difference equations.  or  In a d d i t i o n , the s t a t i s t i c a l properties of the system  109  and measurements need to be modelled by f i r s t and second order statistics (expected values and covariances). results are  not always obtained by a specific formulation of a problem.  may be necessary to obtained.  It should also be noted that satisfactory It  try different approaches before satisfactory results are  For example, different smoothing techniques may need to be considered  to improve the estimates. System Error Sources  System  Measurement Error Sources System State ' x(t)  , MaaCM romant  Apriori Statistics  Observations z(t)  Estimate of System State  1  Kalman Filter i  x(t)  x(t) P(t) Covariance of System State . Estimate  P(t)  Figure 4 . 1 .  Block Diagram of System, Measurements, and Kalman Filter  For the sake of completeness, a summary of the continuous and discrete KF are given in Sections B and C, respectively.  Since the modelling of the  physical system is an important aspect of the KF approach, Section D gives a fairly detailed discussion of the derivation of differential for a typical acceleration measuring device.  equations  Sections E and F discuss the  specific KF equations for the estimation of P-wave and S-wave amplitudes and arrival times.  110  A.  Kalman F i l t e r Results As stated e a r l i e r , when applying the KF i t i s important to model the  physical problem to f i t into the framework of the KF.  The KF f i l t e r proposed  would attempt to model a seismic receiver as a second order system with i t ' s c h a r a c t e r i s t i c natural frequency and damping.  The modelling of a seismic  receiver r e s u l t s i n the necessity of thorough c a l i b r a t i o n and performance indexing.  Guigne'(1986) also i l l u s t r a t e s the importance of instrument  c a l i b r a t i o n when applying his  Acoustic Sub-Seabed i n t e r r o g a t o r , he s t a t e s ,  "The spectral response of the transmitter and receiver i s a l s o being pursued to refine their overall s e n s i t i v i t y .  Providing these parameters are c a l i b r a t e d ,  i t i s believed that dynamic modelling of the truncated parametric array can be used as a r e l i a b l e approach f o r rapid measurements of attentuation i n marine sediments." The seismic receiver acts as a transform function in respect to the input of noise and body waves propagating through the e a r t h .  This i n t e r a c t i o n i s  i l l u s t r a t e d i n Figure 4.2 when we have n o i s e , state x ( 5 ) , and signal (amplitude ^ x ( l ) and a r r i v a l time ^x(2)) i n t e r a c t i n g with the instrument which i s represented by the second order system F(s) with natural frequency w , and damping ? . The output i s a superposition of the n o i s e , x ( 5 ) , passing through the instrument and giving us the second order responses x(6) and x ( 7 ) .  The body waves (S and  P-waves) represented by states x ( l ) to x(4) r e s u l t in the second order responses x(8) and x ( 9 ) .  Therefore, knowing the output, transform function  i t i s desired to determine the c h a r a c t e r i s t i c s of the input signal  (instrument), (i.e.,  maximum amplitude x ( l ) and a r r i v a l time x ( 2 ) ) . The modelling of a geophone was accomplished by applying the program KALDAT which i s given i n Appendix C,  The n o i s e , state x ( 5 ) , was modelled  by passing white noise through the earth, which acts as a low pass f i l t e r . The f i l t e r i n g e f f e c t of the earth is. due to i t s damping c h a r a c t e r i s t i c s , which, i s t y p i c a l l y modelled as a f i r s t order system ( i . e . , exponential decay).  The extent of the decay i s represented by the time constant Tc,  where an increasing time constant applies greater attenuation.  Figure 4.3  i l l u s t r a t e s the e f f e c t an increasing time constant has on a white noise input.  As shown, an increase i n Tc tends to smooth the signal i n d i c a t i n g  a removal of the higher frequencies.  Figure 4.4 i l l u s t r a t e s how the  second order system responds to varying input noise.  Figure 4.4 indicates  that as the time constant decreases, the instrument gives a higher frequency response.. Ill  Acceleration  Input Noise  %P-wave*\^  x^=arrival time  X2=arrival time  Geophone G(s) = ~2 5  Figure. 4.2.  + 2  *V  + u  o  Schematic I l l u s t r a t i o n of the States, x^ for the KF Formulation  112  SYNTHETIC  I  0.0  DATA:INPUT  I  I  0.5  1.0  NOISE  1  1.5  Time (sec)  Figure 4 . 3 .  The E f f e c t An Increasing Time Constant Has On White Noise (Note: T = .003, 005, . 0 1 , . 0 4 , . 0 7 , . 0 9 , . 1 , .2) c  113  i _  2.0  Time(.sec)  Input Noise T  f*°  = .01  Geophone Response  Time (sec)  Figure 4 . 4 .  Time (sec)  Time (sec)  Low-Pass Noise with Corresponding Geophone Response as the Constant Increases  114  Figure 4.4.  Low-Pass Noise with Corresponding Geophone Response as the Time Constant Increases (Cont d)  Figure 4.5 i l l u s t r a t e s the geophones response from the input of the P and S-wave impulses and noise with increasing variance.  As i s i l l u s t r a t e d ,  the impulses r e s u l t in an i n i t i a l displacement than decay which i s a function of the damping and natural frequency of the instrument.  Also shown i s that  as the variance of the noise increases, we have a greater geophone response due to noise ( i . e . , Figure 4.5(a) -  (c)).  The performance of the Kalman F i l t e r was analyzed by generating a source wavelet and passing i t through the second order instrumentation.  The input  wavelet and corresponding geophone response i s i l l u s t r a t e d in Figure 4 . 6 . impulse was generated by the following coding: REAL X1(4000),X3(4000),F1,Q,PIE PIE=3. 1459265 DO a 1=1,200 X1(I)=0.0 1 CONTINUE DO 3 1=66, 76 Q=.01*(I-66) Xl(i)=2500.*C0S(Fl*Q)*EXP(-Q) 3 CONTINUE DO 2 1=1,200 WRITE(1,*)  X1(I)  2 CONTINUE STOP END  116  The  Geophone Response  Input Noise  Figure 4 . 5 .  Geophone Response Due to P and S-Wave Impulses and Noise with Increasing Variance 117  Input Noise  Time (sec)  Figure 4 . 5 .  Geophone Response  Time (sec)  Geophone Response Due to P and S-Wave Impulses and Noise with increasing Variance (Cont'd!"  INPUT DATA  2000Y  1000 V  to § Or  -1000Y  -2000  TIME  SO (N-200.DELT-.01,TC=.  100  SYNTHETIC  50 TIME  Figure 4 . 6 .  150 10.VARlANCE-.09.R0-.ia) DATA:  100 (N=200.DELT=.0i.TC=.  Z  ISO 10,WO=8)  I l l u s t r a t i n g the Synthetic (a) Input and (b) Geophone Response 119  200  200  The signal i s defined by a 50 Hertz cosine wavelet exponentially decaying at .01 and having a maximum amplitude of 2500.  The signal i s i n i t i a t e d  at 66 time units and ends at 76 time units with the t o t a l trace l a s t i n g for 200 time u n i t s . The second order response i s then fed into the Kalman F i l t e r with the output shown in Figures 4.6 and 4 . 7 .  The desired parameters in this run are  x(l),  the signal a r r i v a l time.  the signal amplitude, and x(2),  In the KF i t  i s necessary to give a p r i o r i estimates of the above two s t a t e s , these i n i t i a l estimates were x(l)=2500 and x(2)=75 time u n i t s .  Therefore, we put the exact  amplitude of the input signal and deviated from the a r r i v a l time by 9 time , units ( i . e . , the actual a r r i v a l time was 66 time u n i t s ) .  In Figure 4.7 we  have eight i t e r a t i o n s or passes made through the f i l t e r , with each subsequent i t e r a t i o n having i t s i n i t i a l state value given by the l a s t i t e r a t i o n where the mean square e r r o r , Figure 4 . 8 , was minimum. time, state x ( 2 ) ,  The r e s u l t s show that the a r r i v a l  approached the actual value (66 time units) and reached i t  a f t e r seven i t e r a t i o n s .  The amplitude doesn't change s i g n i f i c a n t l y because  i t was already defined to be the correct value. The above r e s u l t s are s i g n i f i c a n t , because if. t h i s f i l t e r can be f i n e l y tuned , i t would be possible to obtain a r r i v a l times and amplitudes on l i n e r e s u l t i n g in v e l o c i t i e s and damping c h a r a c t e r i s t i c s , r e s p e c t i v e l y .  In a d d i t i o n ,  i t would also be possible to model and c a l i b r a t e the sensors used in data a c q u i s i t i o n , which i s necessary in any type of d i g i t a l processing. Additional e f f o r t needs to be made to v e r i f y the numerical r e s u l t s and models.  In a d d i t i o n , i t w i l l be necessary to investigate smoothing techniques  to obtain reasonable estimates of the P-wave and S-wave parametres x ( l )  to  x(4).  Since the e f f e c t s of the P-wave and S-wave wavelets pass through measurement i n t e r v a l i s a r e l a t i v e l y short time, recognition of t h i s fact to derive a suitable smoothing version of the KF needs to be made.  In a d d i t i o n , further  i n v e s t i g a t i o n needs to be made to optimize the method of working multiple passes of the data through the KF.  120  UPDATED STATES XZ (X0(2)-75  400 900 DISCRETE UNIT TIME (sampling  . Iterations - 3)  1200 interval".01)  1600  UPDATED STATES XI (XO(I)-ZSOO. Iterations - B)  400 BOO DISCRETE UNIT TIME (sampling  Figure 4 . 7 .  1200 interval-.01)  Updating the States (a) x(2) ( A r r i v a l and (b) x ( l ) (Amplitude)  121  1600  Time)  GENERATED UEAN SQUARE ERROR P2 (P0(2.2)-.990.Itrrations-3)  400 300 1200 DISCRETE UNIT TIUE (sampling intrrval-.OI)  1800  GENERATED UEAN SQUARE ERROR P1(P0(1.1 )-.0O0t.ItmtionM-3)  9.9999981S-0S  § 9.99994S2E-0S •«  I.  | * 9.99aa9S7E-0S  9.99934S2E-0S •  400 300 DISCRETE UNIT TIUE (sampling  1200 xntmal-01)  ISOO  Figure 4 . 8 . Mean Square Error f o r (a) State x(2) and (b) State x ( l )  122  B.  Description of the Continuous Version of the Kalman F i l t e r The Kalman F i l t e r , KF, i s a method for estimating a state vector x  from measured data z_.  The state vector may be corrupted by a noise  vector w and the measurement vector i s corrupted by a noise vector v_. KF f i l t e r i s applicable f o r systems that can be described by a f i r s t d i f f e r e n t i a l equation i n x. and a l i n e a r (matrix) equation in z.  The order  The state  and measurement equations assumed to be in the following canonical form: x ( t ) = F ( t ) x ( t ) + G(t)w(t)  (4.1)  z ( t ) = H ( t ) x ( t ) + v(t)  (4.2)  where x^ i s an n - v e c t o r , w i s a p - v e c t o r , z and v_ are m-vectors, and F, G and H are known matrices. The random (vector) processes w and v_ are assumed to be zero mean, white noise processes, so that E(w(t))=0 and E(v(t))=0  (4.3a)  E(w(t)w(t) ) = Q ( t ) 6 ( t - r )  (4.3b)  E(v(t)v(x) ) = R ( t ) « ( t - r )  (4.3c)  T  T  It i s further assumed that w and v_ are s t a t i s t i c a l l y independent of each other,  so that E(w(t)v(r) ) T  = 0  (4.3d)  In the above equations the superscript T denotes the vector transpose, the E denotes the expected value operation, and <5(t-r) denotes the Dirac delta function. Note:  The KF can be modified to treat the case when equations (4.3c) and (4.3d) are not t r u e ; however, the notation i s s i m p l i f i e d considerably when these assumptions are made.  123  It i s further assumed that the i n i t i a l s t a t e , x(t ) H X ^ , i s a zero mean, Gaussian random process and that the covariance of p o s i t i v e s e m i - d e f i n i t e nxn matrix Po.  Thus we have  E(x(t ))=0 a n d E ( x ( t ) x ( t ) ) = P Q  i s s p e c i f i e d as a  o  Q  T  (4.4)  q  We use the notation x_(t'|t) to denote an estimate of the state at a time t ' , t'>t , based on measurements Z ( T ) on the i n t e r v a l t <x<t. —  O  to make the following d e f i n i t i o n s : Definition:  It i s convenient  Q— —  —  1.  If t ' > t , x_(t'|t) i s a predicted estimate.  2.  If t ' = t , xjt'|t) i s a f i l t e r e d estimate.  3.  If t ' < t , x.(t'|t) i  s  a  smoothed estimate.  The estimation error i s defined as x(t'|t) = x(t')-x(t'|t)  .  (4.5a)  and the performance index of the estimation i s given by I(x(t'|t)) = E(e(x(t'|t)))  (4.5b)  where E ( . ) i s an error function of the estimation e r r o r . Although more general error functions can be assumed, we s h a l l assume the e i s the square of 7 ( t ' | t ) ; i . e . , I(x(t'|t)) = E(xU'|t)x(t'|t) ) T  (4.5c)  The estimation problem can now be stated as f o l l o w s : Estimation Problem:  Given the system defined by Equations 4.1 and 4.2 with  s t a t i s t i c s defined by Equation 4 . 3 , and measurements z{x) over the i n t e r v a l t ^ r < t , determine an estimate of £ ( t ' | t ) such that I ( 7 ( t ' | t ) ) i s minimized. The s o l u t i o n to t h i s problem was given by R. E. Kalman (Reference 11) in 1960.  This s o l u t i o n i s s p e c i f i e d by the following equations for the  f i l t e r e d estimate: 124  State Estimation equation x(t)  F(t) x(t) + K ( t ) ( z ( t ) - H ( t ) x ( t ) ) .  with  (4.6) given  Kalman Gain matrix K(t)  P(t)H(t)  R(t)  T  -1  (4.7)  Estimation Error Covariance (Matrix R i c a t t i ) equation P(t)  = F(t)P(t) + P ( t ) F ( t )  T  + G(t)Q(t)G(t)  T  - K(t)R(t)K(t)  T  (4.8)  with P(t ) Q  Note:  = P , given Q  The equations f o r the predicted and smoothed estimates are not given at t h i s time; however, i t appears that smoothed estimates should be used to obtain s a t i s f a c t o r y r e s u l t s the our p a r t i c u l a r problem. A block diagram of the system, measurement, and f i l t e r e d estimates i s  shown in Figure 4.9..  125  input matrix  state noi se  x(t)  measurement matrix  x(t)  OS(t)  >  W(t)  state vector x(t) initial  measurements  measurement  H(t)  z(t)  state matrix F(t)  Figure 4 . 9 a .  measurement noi se v(t)  state  Block Diagram of State and Measurement Equations (4.1 ) and (4.2 ). The Sol id Lines Indicate Vector Quantities  Kalman Gain  state estimate x(t)  measurement estimates  x(t)  z(t) •2(t)  F(t)  apnori estimate  P(t)  x —o apnori covariance P  covariance calculations  Figure 4.9b.  Block Diagram of Continuous Kalman F i l t e r , (4.6 ) and (4.7 )  126  Equations  C.  Description of Discrete Version of the Kalman F i l t e r In most cases i t i s more p r a c t i c a l to formulate the KF in d i s c r e t e  form; t h i s i s e s p e c i a l l y true for computational purposes.  Sometimes the  system and measurement equations are n a t u r a l l y in d i s c r e t e form; at other times i t may be necessary to transform the continuous equations into d i s crete form.  It i s assumed that the d i s c r e t e state and measurement equa-  t i o n s , can be written as  where x ^  +1  x  k + 1  = •(k+l,k)x  z  k + 1  - H(k l)x +  i s an n-vector, w  k  k + 1  + r(k+l,k)w v  +  k + 1  (4.11)  k  "  (4.12)  i s a p - v e c t o r , and z ^  k  +1  and v_^ ^ are m-vectors. +  The nxn matrix $(k+l,k) i s the state t r a n s i t i o n m a t r i x ; t h i s matrix can be obtained from the continuous case (Equation 4.1) as the s o l u t i o n to the homogeneous (zero input) s o l u t i o n to x = F x. the input t r a n s i t i o n matrix.  The nxp matrix r ( k + l , k ) i s  The subscript k refers to the d i s c r e t e time  t = t , k=0, 1 k  The random (vector) processes w^ and  are assumed to be zero mean,  white noise sequences, so that E^)  = 0 and E ( v  E^w/) E(v  j + 1  = Q  k  vj ) +1  *  k + 1  ) = 0  (4.13b)  j k  - R  (4.13a)  k + 1  *  (4.13c)  j k  It also assumed that w and v_ are s t a t i s t i c a l l y independent of each other, k  k  and that the i n i t i a l state x.(t )=x^ and that the covariance of x matrix.  i s a zero mean, Gaussian random process  i s s p e c i f i e d as a p o s i t i v e semidefinite nxn  Thus we have E(ljwJ) = 0 E(x.)  = 0 and E ^ )  (4.13d) =P 127  Q  (4.14)  As i n the continuous c a s e , we use the notation x.(k|j) to denote an estimate of the s t a t e at time t ^ , t ^ t , based on measurements 2^. on the interval t ^ t ^ t j .  As i n the continuous c a s e , i t i s convenient to make  the f o l l o w i n g d e f i n i t i o n s :  f i l t e r e d estimate. smoothed e s t i m a t e . The e s t i m a t i o n e r r o r i s defined as x(k|j) * X|< - x(k|j)  (4.15a)  and the performance index of the e s t i m a t i o n i s given by I(x(k|j)  = E(e(x(k|j)))  (4.15b)  where E ( - ) i s an e r r o r f u n c t i o n of the e s t i m a t i o n e r r o r s .  We s h a l l choose  e to be the square of 7(k | j ) , i . e . , I(£(k|j)) = E(x(k|j)x(k|j) )  (4.15c)  T  The d i s c r e t e e s t i m a t i o n problem can be stated as f o l l o w s : Estimation Problem:  Given the d i s c r e t e system defined by Equations 4.11  and 4.12 with s t a t i s t i c s defined by Equations 4 . 1 3 , and measurements z_. over the i n t e r v a l t . f t ^ _ f t j , determine an estimate x_(k|j) such that I(x_(k|j)) i s 0  minimized. Estimation Equations f o r the F i l t e r e d E s t i m a t e , k=j State Estimation equation  £  k + 1  = •(k+l,k)x  k  + K  k +  with  id  k +  r k l H  +  *(  k + 1  >k)i ) k  (4.15)  128  Kalman Gain Matrix k+1  K  =  p  (k+l|k)H(k+l) (H(k+l)P(k+l|k)H(k+l) T  Estimation Error Covariance (Matrix Ricatti) P(k+l|k) = *(k+l,k)P(k|k)*(k+l,k)  T  T  + R(k+l))  (4.17)  _1  equation  + r(k+l,k)Q r(k+l,k)  T  (4.18) P(k+l|k+l) = (l-K  H  k+1  k+1  )P(k+i|k)  with P(0|0) = P  Q  given.  It is noted that P(k+l|k+l) is the covariance of ^ =_x^ i~^<+i• +1  +  Note: The equations for the prediction and smoothed estimates are not given at this time.  Application of smoothed estimates to our problem is discussed  in Section VI. The computation sequence for the discrete KF is as follows: •  At t  •  At t j ,  Q  specify^, P  and Q , and compute *(1,0), H^, and R^. Q  Q  compute the projected estimate of the covariance matrix  P(l|0)=*(l,0)P *(l,0) +Q . T  o  o  0  Compute the gain matrix K =P(l| 0 ) H ^ ( H P ( 1 , 0 ) H | + R ) "  •  Using the measurement  1  at t  1  at t=t^,  x  = •d.0)x  (J  + K^-H^l.O)^)  The estimation covariance matrix at t^ is given by 129  1  1  the best estimate of the state  is given by l  t  1  P ( l | l ) = P(1|0) •  K H P(1|0) 1  At t = t , a new measurement 2  c y c l e i s repeated.  130  1  i s o b t a i n e d , and the computational  D.  Derivation of Transfer Function for Geophone An electromagnatic device whose output i s proportional to the input  a c c e l e r a t i o n of i t s  case i s shown schematically in Figure 4 . 1 0 .  When the  case of the device moves, in the d i r e c t i o n of i t s s e n s i t i v e a x i s , the r e l a t i v e motion of the c o i l and the magnetic f i e l d induces a voltage across the c o i l .  The recorder across the c o i l measures the c o i l v o l t a g e , e , as a Q  f u n c t i o n of time.  Figure 4 . 1 0 . Schematic of an A c c e l e r a t i o n Measuring Device (Reference 12 ) In Figure 4 . 1 , the f o l l o w i n g parametres are d e f i n e d : M  -  mass of suspended magnet  B  -  flux  K  -  s t i f f n e s s of spring which supports the magnet  D  -  viscous damping between magnet and the case  N  -  number of turns in c o i l  d e n s i t y between magnets  131  d  diameter of coil  I  length of c o i l , ir-d-N  R  resistance of coil inductance of coil resistance of recorder acceleration of case relative to a fixed reference  a (t) c  equilibrium position of M relative to the case dynamic displacement of M relative to the case  The dynamic equations for the device can be derived by summing the voltages induced in the electrical measuring c i r c u i t .  These two situations  are shown schematically in Figure 4.11.  I U/(</<t) + x'l c  Ac 4-  Jilx  r.T i  - mi  Kx  Eft  (b) measuring circuit  (c) forces on M  Figure 4.11. Schematic of Forces on Magnet and Voltages in Measurement Ci rcuit 132  In Figure 4.4, the following variables are used: k  - velocity of M  x  - acceleration of M  i  - current in measurement circuit  BJii  - electrically induced force on the coil  The direction of i is such that a negative voltage is induced in the coil when  M moves downward (x>0).  Summing the forces acting on M gives  Mx + Dx + Kx = B£i - Ma  ^  c  2 1 a  ^  and, summing the voltages in the circuit gives  dT c  L  + R  1 + e  o "** o V =  B  e  (- )  =  4 21b  It is noted that i f L-K) in (4.21b), then iR -B«.x, i . e . , i is proportional to x and / i dt is proportional to x.  Also (4.21a) becomes x+x(D+B £ /R)/M + 2  2  K x/M = - a . Taking the Laplace Transform of these equations gives c  (Ms +Ds+K)X(s) = -Bt I(s)- MA (s) 2  (4.22a)  c  (Ls+R +R )I(s) = - B i sX(s) c  (4.22b)  r  E (s) = R I(s) 0  (4.22c)  r  Since our goal is to obtain the transfer function of the output, E (s), Q  over the input, A (s), that is T(s) = E (s)/A (s) 0  c  it is necessary to eliminate X(s) and I(s)  in Equations (4.22).  accomplished by substituting X(s) from (4.22b) into (4.22a). algebraic manipulations we obtain the desired result 133  This is  After some  T(s)  = E (s)/A (s) = R BAMs/P(s) 0  c  (4.24a)  r  where P(s) = MLs + (DL+MR)s + (DR+KL+B * )s + KR  (4.24b)  R = R + R c r  (4.24c)  3  2  2  2  This derivation of the transfer function for a typical accelerometer helps one obtain a better feel for the relationship between the geophone dynamics and the second order model speicifications given in Section  II.  However, one needs to consider the actual geophone that is used to collect data when one eventually tries to process some of this data with a Kalman Filter. It is noticed that the current in the coil circuit (or voltage across the load resister e =i*R ) is used as the output variable. Q  r  However, in this  case the transfer function that results has an s (Laplace operator) in the numerator - - this is not convenient for applying a second order model. The integral of this output voltage is more convenient to use (this is analogous to the charge across a capacitor). Let q (t)  = J  e (t)  = q (t)  0  e (t')dt'  (4.25a)  0  Then  Q  Q  and E (s) - sQ (s) Q  Q  (4.25b)  i . e . , i t is assumed that the measurement is the integral of the output vol tage.  In this case the transfer function is sQ (s)  Ks  o  Q (s)  r  KJJT rur> =  K;  n  o r  orr ^T  134  =  (4  -  26)  where  s R B£M f  and P(s) is as defined above.  Next i t is assumed that the coil inductance, L, can be neglected. This approximation is valid for low frequency input accelerations, but not necessarily for high frequency, inputs.  It might be informative,to deter-  mine the frequency limits at which this approximation holds when considering P and S waves inputs. In this case P(s) = MRs + (DR+B * )s + KR 2  2  2  (4.27a)  so that the transfer functions can be written as Q (s)  K  r  • °, . = _ E c^ ' s +2C<" s+«  (4.27b)  s  0  0  where K  c  = BAR/R r  w = K/M o  (4.27c)  2  c = (D+B ^ /R)/2Mw 2  In Section II,  2  Q  Table 2.1, the geophone specifications are given as  and 5=0.18.  135  u ~2?*28.0 0  E.  KF Equations for Estimating Amplitudes and Arrival Times In this Section the KF equation for estimating P-wave and S-wave ampli-  tudes, Ap and A , and arrival times, Tp and T , at the geophone location are s  s  derived. In Appendix D the details of the derivation of a discrete model of a second order system with white noise inputs is given.  Comparing the transfer  function obtained in the" previous section (Equation 4.27) with the transfer function in Reference 24 i t is seen that they are in the same form, with the following relationship  K  c  a  c  ( t )  =  w  o  w  {  t  )  (4.28)  or w(t)  = K a (t)/^ c  c  In other words, with the proper scaling the solution equation for the discrete problem (Equation (D.6.f), Reference 24 can be applied here.  A block diagram  of the transfer function is shown below.  K, c  u> and ; are constants Q  magnet position  scaled acceleration w(t)  y(t) s  input  + 2;u» s + a) o o  9-  2  output w(t)  Figure 4.12.  2  = K  c  a (t)/" c  Block Diagram of Second Order Geophone Transfer Function  136  Remember now that P (t)  is an output of the instrument that is related  Q  to the integral of the voltage across a load resistor R in the coil c i r r  cuit.  This could be the voltage across some sort of a capacitance circuit  i t depends on how a particular instrument is designed. One could assume that the measurements (for the KF formulation) are cfU )+v(t ), where v(t ) k  k  k  is white noise which could be due to computer  roundoff, quantization, circuit noise, etc. to assume that q(t ) k  Actually i t might be better  times a scale factor (say h) is observed to be  more general in the KF formulation - - h can always set to 1.0 later on. NOTE:  The quantity q(t) will be proportional to what would probably called geophone "position". If you have a geophone that gives "acceleration", then this will be q(t),  or the second component (y^U)) in the state vector formula-  tion which follows. With the above comments in mind, we define y(t)  = q (t)  (4.29a)  y(t)  + 2 « y(t)  Q  then we can write  C  y(t)  0  = K. a (t)  (4.29b)  c  Next define j  y (t) :  = y(t),  y (t) 2  = y(t),  i  ={y ,y ) l  2  > and  (4.29c) u(t)  = a (t) c  so that the following vector differential equation (d.e.) can be written y = F y + b u  (4.29d)  137  Note that y(t)  corresponds to the instrument output, which is referred to  as the "position". NOTE:  The geophone case acceleration a  c  is assumed to consist of back-  ground noise and P-wave and S-wave accelerations.  Since the geo-  phone transfer function is linear, the noise and signal (P-wave and S-wave) can be treated separately. For the KF formulation define the following state variables  Xg = geophone input due to background noise ( f i r s t order correlation will be assumed) Xg = geophone output due to background noise  x  7  =  (x,-)  *6  Xg = geophone output due to P&S waves  x  NOTE:  x, & x  7  9  =  x  8  and x „ & x  Q  are represented by the same second order d.e.  If we assume that the input noise and P-wave and S-wave acceleration "impulses" are constant over a sampling interval A, then Equation (4.6f) (Reference 4) can be used to obtain the discrete equation. P and S waveforms can be more general than "impulses".  138  Actually, the  "x (k+D" 6  ll  a  =  ^  6  _x  _ 21 2 2 .  7  a n c  "x  12  a  ^  s  (k)"  V  x (k) 5  _x (k+l)_  where ^ t ^ j " ^  a  t  h  ?  a  9 P o  e  e o  n  n e  (k)_  (4.30a)  - 21.  input noise.  b  The constants a^j  b.- are defined by Equation (4.6f) of Reference 4, and will be discussed ij later. b  a  n  The variables x^ and Xg are governed by the same discrete equation, except that i t is clearer to write the outputs as a function of tp and t < s  Thus (assuming tp<t ). $  "o"  "Xg (k+1)" It  (4.30b)  _Xg (k+l)_  "x  8  (k+1)'  _0_  " n" b  x  _x  g  NOTE:  (k+l)_  I f  (4.30c)  b  5-function for the P-wave.  _Xg (k+l)_  (k+l)/A  - 12-  X£(k+1) = tp and x^(k+l) = ap.  "Xg (k+1)"  1  a  We divide by A to approximate a  Generalizations are discussed later.  ll  a  12  _ 21 22a  a  Vl s s t  139  'xg(k)'  (4.30d)  "x  (k+1)"  8  a  "x (k)"  l l 12  _x  (k+l)_  g  NOTE:  x (k+l) = t 4  a n d x3 (k+1)  s  =  4  . 21.  0  b  s  Thus the "nonlinearity" occurs in the times t x (k+l)).  x,(k)'/4  Xg(k)_  a  a  b  +  - 21 2 2 . a  " u"  8  a  —  and t  $  ( i . e . , x (k+l) and 2  However, this nonlinearity is subtle, and will be discussed in  more detail in a later section. It will be assumed that the input noise to the instrument (geophone) is f i r s t order correlated, i . e . ,  x  5  where Y = l / T , and T  +  c  Y X  5  =  Y W  i t is assumed that (4.31a)  a  is the time constant of the input noise to the instru-  ment, and w is a white noise process with 3  Exp(w ) = 0 and Cov(w ) = q a a a a  a  This noise model is discussed in more detail in a separate section.  (4.31b) If  w is assumed constant over a sampling interval, A , then a discrete model a for the input noise can be written as 3  x (k+l) = a x (k) + b w (k) 5  w  5  w  a  (4,31c)  where a  w  = e"  - A, Y  Y A  and b„ = ( l - e " ) TU  140  (4.31d)  Note:  to x^ are assumed to be constants, i . e . , x.(k+l) = x.(k)  i = 1,2,3,4  One can always assume that they are random constants and associate expected values, variances, etc. , to them. Defining the state vector x_ as _x = ( x , x , . . . , x ) 1  2  9  (4.32a)  T  One can write the discrate state recurrence equation as - I  (JVV  +  (4.32b)  where \  +  l  = (0,0,0,0,b ,0,0,0)w (k+l) w  a  (4.32c)  the uncorrelated input noise to the geophone The vector jF needs to be written as a function of x^(k+l) and x^(k+l) t  (i.e.,  and t ) to see how i t is incorporated into the KF. $  The measurement equation for this problem is assumed to be (as discussed above)  i . e . , assume we are measuring velocity. z(k+l) = H x_(k+l) + v(k+l)  (4.33a)  where H = (0 0 0 0 0 0 1 0 l)m fa is a scalar) y  cov (v(k+l))= r  ra  v  and exp (v(k+l)) = 0  (4.33b) (4.33c)  i . e . , we have a scalar measurement m (x^+Xg) with uncorrelated measurement v  noise having variance r m  141  We w r i t e  1(*K>V  =  F  (4.34a)  *k  where  F =  (4.34b)  B  s  o  A;  0  0  where 1 0 I  c  =  0  1 0  0  0  1 0  0  0  0 1  b  ll  L 21 b  a  a  a  an i d e n t i t y matrix to  represent  the 4 constant state v a r i a b l e s  12  the a . , and b , .  s p e c i f i e d above  21 22.  0  '21  It  l l  0  11  0  q  b  b  i i o" 0  21  a  and A' = s  ll  _ 21 a  a  12~  a  22_  i s noted that the c o e f f i c i e n t s b'. . , b'.'. depend on the constants x ( k + l ) ?  and x ^ ( k + l ) .  In the KF implementation one makes estimates of x  t=0, and c a l l them x ( 0 ) 2  until Note:  t  k + 1  ix (0), 2  & x (0).  = x.(0)  One can see the d i s t i n c t i o n with  (i=2,4) f o r  t  k + 1  estimates x\(0)  at  <x (0). 2  between generating the geophone data  the models, and estimating the parameters with the  In one c a s e , x ^ O ) ,  & x,  The KF w i l l not change these estimates  4  and x . ( k + l )  2  models.  i=l to 4 are known, and in the other case, the  are used to obtain updated estimates x ^ k + 1 ) , i = l to 4. 142  If x.(T)  denotes the  estimates of x.,  i  s  1 to 4, at the end of  the geophone trace, then one could reprocess the data with the KF using x\(0)  = x\(T)  for the next pass.  Hopefully this  process will converge!! Note:  It is convenient to define F i n terms of matrices I . A , B , and c n s A because when you propagate ahead in the KF you need to make T s  operations likeFPF , and, by partitioning fications over the intervals is the end time of  (O^),  ^JX^),  the geophone trace.  143  F, one can make simplia n  d ( ^,T), where T X  The Kalman Filter Equations for Velocity Measurements We have the state equation, at time t  • KvV  +  \  k+1  =(k+l)  ,  n  (  4  -  3  5  3  5  b  )  a  )  where w  = (0,0,0,0,b ,0,0,0,0) w (k+l) T  w  | ( + 1  cov(w (k+l)) = q  a  a  fl  and the measurement equation z  k l " iik 1 H  +  +  +  v  m^  k+1  )  c o v  (  v  m  (  k + 1  )  =m r  (  4  '  where H = (0 0 0 0 0 0 m 0 m ) y  NOTE:  w (k+1) and  v  m  ( k + l ) are both zero mean, white noise sequences.  Let P denote the 9x9 state estimate covariance matrix Assume i n i t i a l  estimate of the state and estimation matrix, x^ and P , are  available. Step 1 Compute $>=3f/3_x i f x (k+l) <_ t 2  $ (k+l) < t 4  L + 1  <_x (k+l) + N A. 4  s  k + 1  <_ x ( k + l ) + N A, or i f 4  p  Otherwise, *=F.  Step 2 Propagate the estimate of the state ahead one step, i . e . ,  » V i • I<VV  <4  -  36a)  Step 3 Propagate Covariance matrix ahead  P  'k+1  =  *k * P  T+  (4.36b)  Q  2 where a l l the elements of Q are zeros except Q(5,5) = b q W  144  Q  Step 4 Compute the updated measurement noise covariance  (4.36c) -  »»  ( 77 P  +  P  79  t  97  P  The p^j refer to elements of matrix  t  P  99  )  * m r  a  s  c  a  1  a  r  P'^+j  Step 5 Compute the Kalman Gain matrix K,  P ' . •. H  T  T  +  •Pi?  Pi " 9  + P »26  k+l  +  LP 7  (4.36d)  29  P j 9 9  9  Step 6 Compute measurement residuals  A z  k+l  =  z  =  k+l "  z  i'k+l  H  k+r v m  ( x  '7  (4.36e)  ( k + l )  +  x  g'  ( k + l ) )  Step 7 Compute new State Estimate  *k+l  =  x  'k+l  +  K  k+l  A z  k+l  (4.36f)  Step 8 Update State Estimation Error Covariance matrix  P  k+l  =  P  'k+l  "  K  k+l  H  P  'k+l  145  (4.36g)  Note:  Since  H P'  = K  +  j  1  n v  r  K  k+1  H  P  (  p  i 7  l l P l 7  K  (  K  H  Pj^  + 1  is  a  1x9  +  p  l 9  +  P  19  }  K  P Q 9Q 9) ''  +  and  H  n1 1 Q ^ 977 (  p  +  PQQ^ 99 H  J  k+1 - 91  A considerable  P I 9 ' - - - ' P Q 79 7  +  number vector  of  17  ( p  the  and  K^  }  K  computations + 1  is  a  9x1  146  91  ( P  are  97  +  saved  vector.  P  99)  by  using  the  fact  that  F,  Comments on Computing the Estimation Error Covariance M a t r i x , P Since a considerable portion of the c a l c u l a t i o n s required to obtain the  KF estimates are r e l a t e d to updating the Estimation Error Covariance M a t r i x , P, i t i s useful to note possible computational e f f i c i e n c i e s which might be r e a d i l y implemented.  The e s s e n t i a l  computational e f f i c i e n c i e s are r e -  lated to the s t r u c t u r e of the state t r a n s i t i o n m a t r i x , <t>, and to the rather simple structure of the measurement m a t r i x , H.  That i s , the upper l e f t 4x4  portion of the * matrix i s an i d e n t i t y m a t r i x , the upper r i g h t 4x5 matrix i s a n u l l m a t r i x , and the H matrix i s a 9x1 vector with mostly zero elements. The P matrix i s f i r s t updated over a sampling i n t e r v a l , A , to the P' matrix by the f o l l o w i n g equation: * P  P' k+1  (4.37a)  *' + Q  k  K  where, as p r e v i o u s l y noted, a l l elements of the state noise covariance matrix Q  are z e r o , except for the Q(5,5) element.  The state t r a n s i t i o n matrix «  i s of the f o l l o w i n g form:  *44 * = B  54  °45 A  (4.37b)  55  where denotes a 4x4 i d e n t i t y matrix O^g  denotes a 4x5 n u l l matrix  a w b  55  l l  21 0  b  0  0  0 a  a  l l  21 0 0  a  a  12 22 0 0  a  a  and 147  0  0  0  0  0  0  l l  a  12  21  a  22  '54  0  0  0  0  0  0  0  0  0  0  0  0  11  c  ll  b  h  c  h  21  c  21  b  21  C  21  The elements of the  and B^ matrix have been p r e v i o u s l y d e f i n e d . 4  i s noted that i f the time t t  k + 1  k + 1  It  f a l l s into any of the f o l l o w i n g i n t e r v a l s  < x (k l) 2  +  or x (k+l) + N .A 2  p  <t  <x (k+l)  k + 1  4  (4.38)  or t  k + 1  > x ( k l ) + N -A 4  +  s  then B g = 0 ^ , i . e . , a 5x4 n u l l matrix, 4  The  4  updated P matrix f o r time t  P  k 1 +  ••  P  M  "  K  M  H  P  k +  ^ i s given by the f o l l o w i n g equation:  M  ^  3  9  >  where K ^ i s the Kalman Gain matrix given by k+  K  k+1  =  P  (4.40)  k+1 " W + l  where  r  k l = +  H  P  k 1 +  +  r  (4.41)  m 148  To simplify the notation, the time subscripts k and k+1 are dropped for the rest of this section since i t is clear that P and K are being updated from time t^ to time  The subscripts on the matrices  and  will  also be dropped. Since H = m (0 0 0 0 0 0 1 0  1), for the velocity measurement case,  v  it  is evident that  Pi7 Pi7 +  P'H  T  = m.  (4.42a)  L.P97 99 +  p  and  r' = ^ ( P ) P ) P 7 P 9 9 ) m 7  +  9  +  +  (4.42b)  + r  9  where the p'. • are elements of P'.  The gain matrix has the f a i r l y simple  form  Pi7 Pig' +  1 T irP'H r  m  =A r  1  (4.43)  P97 P99 +  It is next necessary to partition the P and P' matrices into the same structure as the partitioning of the * matrix. a  P ^ Lb  t c-  Thus we define  (4.44a) 1  where  149  P, a  denotes a 4x4 symmetric matrix  P. b  denotes a 4x5 matrix  P  denotes a 5x5 symmetric matrix  c  With these definitions i t is seen that the following matrix product is obtained P *  P *  P  a A  a  B +P.A T  b  T  (4.44b)  =  _  BP + A P J  a  b  c  BP B +BP. A +APIB +AP T  b  T  b  T  A  T  c-i  For convenience of notation, the following matrix products are defined: X =  BP A'  b K  Y = P B  a  T  + P. A  b  (4.44c)  T  and Z = BP^B'  c  150  Since the matrix B has three rows which are all zero (top three), i t is readily seen that the only non-zero elements of matrix Z are the bottom right 4 elements.  The details of these calculations are given in  the "Summary of Equations" section.  It is noted that i f  f a l l s into  one of the intervals defined above, then the computations are reduced considerably since X = 0 (4.44d)  Y = P. A b and  It is also noted that some of the calculation involving products with the B matrix are eliminated by the fact that the top three rows of B are all zeros, With the above definition of X, Y, and Z i t is seen that  +Q •Y  (4.45)  Z+X+X +AP A -  T  T  C  T  In the final step of updating the P matrix, i t is convenient make a different partitioning of P'. of the H matrix.  This partitioning is based on the structure  In particular, i t is readily seen that  P = P'  (4.46a)  P"  where (4.46b)  p« = i,. P' (H H)P' T  Since T  H'H  =  2 m v  ro  1  o  2  c  L  °2  '3 151  where denotes a 6x6 n u l l matrix 0  2  denotes a 6x3 n u l l matrix  and  l  3 -  1  0 1  0  0 0  1  0 1  i t i s convenient to p a r t i t i o n P' i n t o the same s t r u c t u r e ,  i  P  P'  where P^, P  2  P  l  P  2  P  3  i.e.,  (4.46c)  and P^ are dimensioned the same as 0 ^ , 0 , and I^, 2  respectively.  With t h i s p a r t i t i o n i n g , i t i s seen that °1  H HP' =  °2  (4.46d)  T  L i P  T  P  3J  and P P'(H'H)P  2  P  2  T  P  2 3 P  (4.46e)  1  L  P ^  P ^  where, i t i s noted that Pj^ and P^ are symmetric matrices. The d e t a i l s of these c a l c u l a t i o n s , and the computational e f f i c i e n c i e s achieved, are discussed in the "Summary of  152  Equations" Section VI.  G.  Evaluation of P a r t i a l s with Respect to  and  According to L e i b n i t z ' s r u l e , we have  1  f(t..)dt-y  f . ,  +  f  U  , . , _ 1 . f(  l  In the present case, ^ ( a ^ x ^ n A - A  dx  + 2 > c  ,) _ i  (4.47)  and <j> (a)=x +nA, n=l to N . 2  2  Thus  iand_j_i = i dx  2  (4.48a)  2  Also f ( t , a ) i s not a f u n c t i o n of a ( i . e . ,  3 f  |i(t..)dt  x ) 2  thus  ( ) = Q 3a  (4.48b)  t > o t  Thus x +nA 2  4-r--  ,/  fir,  J  d a  f(t)dt = f(x„+nA-A)-f(x +nA)  (4.48c)  0  x +riA-A  2  2  1  Since the present problem i s time i n v a r i a n t ( i . e . ,  the s o l u t i o n f o r the  second order system depends only on the d i f f e r e n c e  x +A-x =A, 2  2  we need only  consider the l i m i t s 0 and A . NOTE:  One needs to be c a r e f u l on what the values of a (t) over t h i s i n t e r v a l .  2  are  That i s , the l i m i t s 0 and A d o not change, but  the values of ap(t) <j> (a) = x  (and a ( t ) )  do.  2  Thus the l i m i t s of i n t e g r a t i o n are  + n A  (4.48d) ^(a)  = x  2  + (n-l)A  153  n=l,2,... ,N . p  However, since the system is time-invariant, one needs only to consider the sampling interval, A. Typically, as previously noted, the input wavelets will have the form shown below.(the same  'P  arguments are also true for the S-wave).  W  *  W'  =  w h e r e  p  t h e w  are  predefined scalars  x +NpA 2  T—l  »—f-ff^  time, t.  r-  A  U  p-w (k) p  According to Equation (4.6b) of Reference 4, we have 2  A F-(A-t) e  ~C A - A _<-'  t r  -  sin p'(A-t)  dt-a (t ) p  sinp'(A-t)+cosp*(A-t)  01  Define f(t)  r  =e  - I sinp' (A-t) 01  - c sin p'(A-t) + pcosp'(A-t)  p(V  a  (4.49a)  then -  r  f(0)  -  sin  p  0)  .-c sin p"+ cosp"_ p  yv  (4.49b)  and f(A) = e  0 P  a  P  < W  154  (4.49c)  k  Thus -  f(0)-f(A) =  sinp"  _COS(p"+a)-pe"  A -(w (t. )-w (t. .)) p p k' pVk+1'' v  (4.50a)  A  v  where p" =p'A = p-cu-A (4.50b)  C" =5 A = ?-oj-A k = x /A to x /A + Np 2  Note:  x  2  and x  4  2  should be selected/estimated so that they are multiples of  the sampling time A. Therefore aie  _d_ dx.  f  sinp" • ap(t")  f(t)dt cos( "+a) • ap(t')-paj e" "ap(t ) 2  P  - sin 01  a) e  P  "  2?  ap(t ) (4.51a) -r"  cos(p"+a)  +  +  • ap(t ) - p e ^ ap(t )  where t  = x  2  + n A -A  t  = x  2  + n  (4.51b)  A  n = 1 to Np Alternatively,  i f we use Equation (D.6.e) of Appendix D we have  155  / .2  r" 4  cosp"  I sinp"  sin(p"+a)  ii  - c ' t I sinp't  COS(p"+a) i  icosp't . I  4.  I  dt-w. K  where a = sin Applying similar notation, i . e . ,  c  letting  sinp't f'(t)  (4.52a)  =e COSp't  we obtain f'(0)  =e  f(A)  = e"  (4.52b)  " (J)-p(O 0  ,.  /sinp"\  (4.52c)  c  \cosp" /  p  Thus 0 i  aj e  sinp .a (t')-e' "(-cosp"sinp"+sinp"cosp")a (t ) M  c  p  p  +  2  cos(p"+a).a (t")-e" "(sin +a)sinp"+cos(p"+a)cosp") . ;  n  OJ  e  i s i n p » . a  p  P  a  ( t - )  (4.53a) cos(p"+a)-a (t')- e' "a (t ) p  which is the same result as above.  P  ?  p  +  Since this result holds for both the P-wave  and the S-wave, we can write, for i=2 and 4,  156  c f(t)dt = S_i  A.  J L  f  J  sin ".w.  2  P  OJ  x^+nA-A  ,(t )  (4.53b)  1-1  O  cos( "+a)-w._ (t")-pe" "w _ (t )_ P  ?  1  i  1  +  where w^(t") = x^(t")-Wp(t")  t" = x +n-A-A, n=l to N 2  p  (4.53b) w (t~) = x-,(t~)-w (t~)  Note:  t" = x.+nA-A,  n=l to N  When the w^ are both zero, these partials are zero and there is no need to evaluate them in the computational algorithms.  Further-  more, since these partials are required before updated estimates can be obtained, the estimates from the previous sampling interval, i . e . , x^t'-A),  x (t"-A), x (t -A), 2  and x ( t ' - A )  _  3  4  (4.53c)  are used. If the state vector is defined by  then F is used for updating x^, and $ = tion error covariance matrix P^. $ = F  for  t^<x  and *=F f o r x <_t^x +NA 2  82  2  t  c  > x  ^ N -A +  s  except for the following elements:  4  3F /3x ,  is used for updating the estima-  It is noted that  and |  2  3F , .  3 F / 3 x , 3 F / 3 x , and 3 F / 3 x g 2  2  8 4  4  g4  2  In summary, these elements can be written for i=2 and 4, as  8i = c. L/9iJ  c  2  •w.^t )  c  3  • w-.^t") - c  (4.54a)  157  4  •W ^ t  )  c  = — sin p"  0  p" = pcuA (4.5b)  C, = COS(p"+a)  X  i-1^0 _  w^t*)  p  2  5  at t", i =2 and 4  n = 1 to N  p  = w^t'+A)  w (t") = w (x +nA-A) 3  1  = best estimates of  w (t ) = w (x +nA-A) x  a - sin'  p  4  (4.54c) n = 1 to N  $  w (t ) = w (t +A) 3  +  3  _  The parametres c^, i=l to 4, are constants and can be precompiled. The weighting factors w and w are known functions of time. $  158  V.  SUMMARY & CONCLUSIONS  In t h i s research a new approach of d i g i t a l f i l t e r i n g in determine shear wave and compression wave v e l o c i t i e s i s introduced with great success.  The  technique currently being used at U.B.C. (Campanella and Robertson, 1986) r e l i e s upon the p o l a r i t y of shear waves to determine a r r i v a l times of shear waves at s p e c i f i c depth i n t e r v a l s to c a l c u l a t e shear wave v e l o c i t i e s over each depth i n t e r v a l . The reverse p o l a r i t y method requires the shear waves to be polarized 180 out of phase.  These waves are obtained by generating shear waves on opposite  sides of the seismic source.  Thus, i t i s ' necessary to use two seismic waves  to determine one reference a r r i v a l time as opposed to d i g i t a l f i l t e r i n g methods where only one seismic wave i s necessary.  Once the waves are recorded, i t i s  necessary to v i s u a l l y pick a cross-over or reference a r r i v a l time from the oppositely p o l a r i z e d waves.  The determination of the cross-over can be d i f f i -  c u l t when low or high frequency noise or signal i s present.  Therefore,  analogue f i l t e r i n g i s often applied in order to a s s i s t in the understanding and reduction of the time series recorded.  Stokoe and Hoar (1978) suggest that  f i l t e r i n g should be minimized because i t may s i g n i f i c a n t l y d i s t o r t the signal and erroneously a l t e r a r r i v a l times.  An additional consideration i s that the  reverse p o l a r i t y method demands that in picking a cross-over, a consistent method must be used.  For instance, one may decide that the second cross-over  w i l l be used throughout the seismic i n v e s t i g a t i o n .  But i t has been found that  a s i n g l e c r i t e r i o n may not be appropriate throughout the seismic p r o f i l e due to varying instrument responses and frequencies changing s l i g h t l y with s o i l l a y e r ing. The two f i l t e r i n g methods considered in t h i s thesis are the one based on frequency domain f i l t e r i n g , where Fast F o u r i e r , Butterworth, and c r o s s c o r r e l a t i o n algorithms are implemented, and one based on time domain techniques, where a Kalman F i l t e r i s designed to model the instrument and the physical environment. The c r o s s c o r r e l a t i o n method r e l i e s upon f i l t e r i n g response data such that a s p e c i f i c wave i s analyzed ( e . g . , SH-wave) when applying the crosscorrelation function.  The c r o s s c o r r e l a t i o n function i s applied to two successive waves  (corresponding to the successive depths) where the maximum value of the c r o s s c o r r e l a t i o n function i s d i r e c t l y related to the difference in travel times 159  between the waves. calculated.  With t h i s a r r i v a l time d i f f e r e n c e , a v e l o c i t y can be  The method i s very advantageous because i t allows one to focus  on a s p e c i f i c wavelet and use a l l the information of the wavelets present averaging out any noises or i r r e g u l a r i t i e s and r e l y i n g upon dominant responses as opposed to the reverse p o l a r i t y method where there i s dependence on one point ( i . e . , the cross-over)  to define a r r i v a l time d i f f e r e n c e s .  Some impor-  tant aspects of d i g i t a l f i l t e r i n g to consider when implementing the crossc o r r e l a t i o n function are:  sampling r a t e , Gibbs' phenomenon, a l a i s i n g , the  Butterworth f i l t e r , and d . c . The-crosscorrelation  shifts.  f i l t e r was coded in Fortran language, and the program  i s referred to as CROSSCOR. CROSSCOR i s a graphics i n t e r a c t i v e program which displays the frequency s p e c t r a , u n f i l t e r e d and f i l t e r e d times series and crossc o r r e l a t i o n s on a mainframe graphics terminal which has been adapted to run on the IBM P.C.  CROSSCOR was tested by superpositioning several sinusoidal  waves and f i l t e r i n g out the desired frequencies. data to assess i t s performance and r e l i a b i l i t y .  CROSSCOR was applied to real Results from f i v e research s i t e s  are given, and the v e l o c i t y p r o f i l e s from d i f f e r e n t methods of data a c q u i s i t i o n and reduction are presented and compared along with the respective cone p r o f i l e s . The f i r s t research s i t e analyzes data acquired at the McDonald Farm S i t e in Richmond, B.C.  In t h i s case study f i v e sets of data were compiled and two  d i f f e r e n t sensors were used.  The f i r s t sensor was the Hogentogler Super Cone  Velocity transducer and the second sensor was the UBC Seismic Cone Pressuremeter, which had two accelerometers i n s t a l l e d .  The seismic data was reduced  by the reverse p o l a r i t y method and CROSSCOR, and compared with cone bearing profiles.  The second s i t e was located at the Grant M.cConachie Way in Richmond,  B . C . , where UBC #7 seismic cone was used with an accelerometer sensor i n s t a l l e d . Two d i f f e r e n t types of sources were compared in t h i s i n v e s t i g a t i o n , i . e , the shear and compression waves sources.  The t h i r d research s i t e analyzes data  acquired at the T i l b u r y National Gas Plant S i t e on T i l b u r y Island where the UBC seismic cone pressuremeter was used for data a c q u i s i t i o n .  The fourth s i t e was  the Annacis Bridge P i l e Research S i t e where the UBC Seismic Cone Pressuremeter and UBC #8 Seismic Cone were used for data a c q u i s i t i o n . shear v e l o c i t i e s were obtained using shear sources.  In t h i s  investigation  Since both sides of the  source were used, i t was possible to obtain two v e l o c i t y estimates for one depth increment using CROSSCOR as opposed to the reverse p o l a r i t y method which allows for only one v e l o c i t y estimate. 160  The l a s t research s i t e was the Lower  was the Lower Langley 232nd S t . S i t e .  In t h i s investigation f i l t e r e d  (analogue)  and u n f i l t e r e d reverse polarized Hammer Shear Source data was obtained.  In  a d d i t i o n , point source Buffalo Gun data was reduced with CROSSCOR which resulted in separating and measuring both shear and compression wave v e l o c i t i e s . The results from the synthetic data and research s i t e s indicate that CROSSCOR i s a very b e n e f i c i a l tool to obtain seismic wave v e l o c i t i e s .  The pro-  gram allows one to separate out d i f f e r e n t seismic events and use a l l the i n f o r mation contained in the wavelets in obtaining shear and compression wave v e l o c i ties.  The procedure in extracting the desired instrument responses proved to  be routine and user f r i e n d l y due to the d i f f e r e n t responses being consistent throughout a cone p r o f i l e .  In a d d i t i o n , since CROSSCOR i s a d i g i t a l f i l t e r  it  removes most human bias ( e . g . , determining cross-overs in the reverse p o l a r i t y method) in determining a r r i v a l times, the program allows one to compare v e l o c i t i e s from signals obtained on opposite sides of the source, and one may also have the f l e x i b i l i t y in using d i f f e r e n t sources in obtaining seismic wave velocities.  The problem with t r i g g e r i n g when applying the Buffalo Gun Source  may be resolved by using an accelerometer (having a f a s t response time) as the t r i g g e r as well as grounding the Buffalo Gun f i r i n g pin to create a contact closure ( s i m i l a r to the Hammer Shear Source). The Kalman F i l t e r was applied in the manner in which i t modelled the sensors used and the physical environment of the body waves and noise generation.  The  KF was investigated for i t s possible a p p l i c a t i o n to obtaining accurate estimates on the P-wave and S-wave amplitudes and a r r i v a l times.  The KF i s a very f l e x i -  ble tool which allows one to model the problem considered accurately.  In a d d i -  t i o n , the KF works in the time domain which removes many of the l i m i t a t i o n s of the frequency domain techniques. The Kalman F i l t e r i s an optimal (in a l e a s t squares sense) f i l t e r which i s based on s t a t e - s p a c e , time-domain formulations of physical problems.  Application  of t h i s f i l t e r requires that the physical problem be modelled by a set of f i r s t order d i f f e r e n t i a l or difference equations which, with i n i t i a l c o n d i t i o n s , uniquely define the system behaviour.  The f i l t e r proposed attempts to model a  seismic receiver as a second order system with i t ' s c h a r a c t e r i s t i c natural f r e quency and damping.  The instrument acts as a transform function in respect to  the input of noise and body waves propagating through the earth. 161  The performance of the Kalman F i l t e r was analyzed by generating a source wavelet and passing i t through the second order instrumentation.  The second  order response i s than fed i n t o the Kalman F i l t e r with the a r r i v a l time and maximum amplitude being determined.  The f i l t e r was found to perform well  and i t has much promise in respect that i f i t i s f i n e l y tuned, i t would be possible to obtain a r r i v a l times and amplitudes on l i n e r e s u l t i n g in v e l o c i t i e s and damping c h a r a c t e r i s t i c s , r e s p e c t i v e l y .  There are many other forms of the  f i l t e r which could be implemented such as smoothing.  Smoothing would allow one  to use a l l the information in the seismic trace to determine the location and shape of the wavelets being recorded (which i s a very important parameter in earthquake engineering).  The damping obtained from the KF would be a d e f i n i t e  advantage over conventional techniques which r e l y on power spectrum measurements to determine the damping c h a r a c t e r i s t i c s of a s o i l p r o f i l e with depth.  Since  the power spectrum depends on more variables ( e . g . , i n i t i a l energy of source, frequency content of body waves, and coherency between successive t r a c e s ) , i s greater chance of error than the Kalman F i l t e r formulation.  162  there  VI.  SUGGESTIONS FOR FUTURE RESEARCH  The research in t h i s t h e s i s indicates that d i g i t a l f i l t e r s are extremely b e n e f i c i a l in i n - s i t u t e s t i n g .  As was o u t l i n e d , there are many d i f f e r e n t  types of f i l t e r s one can implement; therefore, i t i s important to become fluent in a l l the d i f f e r e n t f i l t e r i n g techniques and apply the ones which best s u i t the problem considered.  The frequency domain approach i s a well established  method of manipulating the data to obtain desired frequencies, while the state space approach i s the " s t a t e - o f - t h e - a r t " technique in time series a n a l y s i s . Future research at UBC should concentrate having the i n - s i t u group becoming more f l u e n t on the above techniques in d i g i t a l f i l t e r i n g in order to expand the possible goetechnical parametres that can be obtained in i n - s i t u t e s t i n g . of these parametres i s damping of the s o i l .  One  S o i l damping and the shape of the  seismic wavelets are very important parametres in earthquake engineering. Another area of i n t e r e s t in geotechnical engineering i s seismic tomography. Seismic Tomography combines the data from a large number of waves in a given s o i l p r o f i l e to construct a three dimentional image of the medium ( e . g . , veloc i t y and damping) that the waves have t r a v e l l e d through.  Therefore, i f one  establishes accurate methods f o r determining P and S-wave v e l o c i t i e s and attenuation, the next step would be to apply tomographic concepts. Further research at UBC should also concentrate on making newly developed programs user f r i e n d l y .  The programs should be tested extensively in the f i e l d  with industry and e f f o r t made i n t o developing the best, most r e l i a b l e , and e f f i c i e n t techniques and programs f o r seismic data reduction.  163  LIST OF REFERENCES 1.  R i c h a r t , F. E . , J r . ,  H a l l , J . R., J r . ,  and Woods, R. D., (1970),  "Vibrations of S o i l s and Foundations," P r e n t i c e - H a l l , I n c . ,  Englewood  C l i f f s , New J e r s e y , 414 pp. 2.  Mooney, H. M. (1974), "Seismic Shear Waves in Engineering," Journal of the Geotechnical Engineering D i v i s i o n , ASCE, V o l . 100, No. GT8, Aug. Proc, Paper 10745, pp. 905-923.  3.  Borm, G. W., (1977), "Methods from Exploration Seismology:  Reflection,  R e f r a c t i o n , and Borehole Prospecting/' Proceedings of DMSR 77, Karlsruhe, 5-16 Sept. 1977, V o l . 3 , pp. 87-114. 4.  R i c e , A . , 1984, "The Seismic Cone Penetrometer," M.A.Sc. Thesis, Dept. of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver, B.C.  5.  Laing, N. L., 1985, "Sources and Receivers with the Seismic Cone T e s t , " M.A.Sc. T h e s i s , Dept. of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver, B.C.  6.  Close & F r e d r i c k , Modelling and Analysis Dynamic Systems, Houghton Mif1 in C o . , 1978.  7.  Ogata, J . , Discrete Process C o n t r o l , Cambridge University P r e s s , 1987.  8.  Ewing, J a r d e t z , and Press, E l a s t i c Waves in Layered Media, McGraw-Hill Series in the Geological Sciences, 1957.  9.  K. H. Stokoe and R. S. Hoar, "Variables A f f e c t i n g I n - S i t u Seismic Measurements," Proceeding of Conference on Earthquake Engineering and S o i l Dynamics, ASCE Goetechnical Engineering D i v i s i o n Specialty Conference, Pasadena, CA, Vol. I I ,  pp. 919-939, 1978.  10.  Lancaster, Active F i l t e r Cookbook, Howard W. Sands and Company, 1982.  11.  R. E. S h e r i f f and L. P. Geldart, Exploration Seismology, Volume 2 , Data Processing and I n t e r p r e t a t i o n , Cambridge University Press, 1983.  12.  George W. S t i n s o n , Introduction to Airborne Radar, Hughes A i r c r a f t C o . , 1983.  13.  E. R. Kanasevich, Time Sequence Analysis in Geophysics, Third E d i t i o n , the University of Alberta Press, 1981.  164  14.  R. E. Thomson and K. Y. Chow, "Butterworth and Lanczos-Window Cosine Digital F i l t e r s :  with Application to Data Processing on the UNIVAC 1106  Computer," P a c i f i c Marine Science Report 8 0 - 9 , I n s t i t u t e of Ocean Sciences, Sidney, B.C. 15.  E. A. Robinson, " C o l l e c t i o n of Fortran II Programs for F i l t e r i n g and Spectral Analysis of Signal Channel Time S e r i e s , " Chapter 22 of textbook by Unknown.  16.  U.B.C. P l o t t i n g Programs, Unpublished Notes.  17.  Campanella, R. G . , Robertson, P. K., and G i l l e s p i e s , D., 1982, "Cone Penet r a t i o n Testing in D e l t a i c S o i l s , " Canadian Geotechnical J o u r n a l , V o l . 20, February.  18.  Howie, J . A . , 1989, Ph.D. Thesis (in preparation), University of B r i t i s h Columbia, Vancouver,  19.  B.C.  Davies, M . , 1987, " P r e d i c t i n g A x i a l l y and L a t e r a l l y and P i l e Behaviour Using I n - S i t u Testing Methods," M.A.Sc. Thesis, Department of C i v i l i n g , University of B r i t i s h Columbia, Vancouver,  20.  Engineer-  B.C.  G r e i g , J . , 1985, "Estimating Undrained Shear Strength of Clay from Cone Penetration Tests," M.A.Sc, Thesis, Department of C i v i l  Engineering,  University of B r i t i s h Columbia, Vancouver, B.C. 21.  Dobrin, M. B . , 1975, "Introduction to Geophysical Prospecting, McGraw-Hill Company.  22.  Guigne', J . Y . ,  "Acoustic Sub-Seabed Interrogator,"  P h . D . , University  of  Bath-England, 1986. 23.  Asten, M. W., "Theory and Practice of Geophone C a l i b r a t i o n I h - S i t u Using Modified Step Method," IEEE Trans. Geosc., V o l . GE-15, No. 4 , pp. 208-214, October, 1977.  24.  Baziw, E. J . , "Concepts of Kalman F i l t e r i n g and the Application to Seismic Time Series A n a l y s i s , " Bachelors of Applied Science Thesis,  University  of B r i t i s h Columbia, A p r i l 1986. 25.  Kalman, R. E . , "A New Approach to Linear F i l t e r i n g and Prediction Problems," J . Basic Energy, Series D, V o l . 45, pp. 25-36, March 1960.  165  APPENDIX  A  Program l i s t i n g of CROSSCOR  166  Listing 1 2 3 4 5 6 7 8 9 10 1 1 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 4 1 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58  o f CROSSCOR  REAL REAL REAL REAL INTEGER REAL  c  c c  c .  c  c c c c c  on JUL 2 7 , 1988 f o r C C i d = S E I S on G  C I 1 ( 2 0 4 8 ) ,CI 2 ( 2 0 4 8 ) , X ( 2 0 4 8 ) , Y ( 2 0 4 8 ) , W ( 2 0 4 8 ) XX(2048),YY(2048),FC1,FC2,FC3,FC4.D1(8),G1 AVRGX,AVRGY,X2(2048),YZ(2048),D2(8),G2,DELT DLT,G(5000),GA(10000),CMAX,EX(10000) LG,LA,LMAX,HMAX X 1 ( 2 0 4 8 ) ,Y 1 ( 2 0 4 8 ) ,FN  READ(3,*) N1,DELT DLT=1000.0*DELT FN=0.5/DELT DO 1 0 0 1=1,N1 READ)1,*) X ( I ) READ(2,*) Y ( I ) 100 CONTINUE  CALL NORME(N1,X) CALL NORME(N1,Y) CALL REMAV(N1,X,AVRGX) CALL REMAV(N1,Y,AVRGY) CALL P0W2IN1.X.Y.N)  c  c  a t 16:44:17  DO 1 0 2 1=1.N1 XX(I)=X(I) YY(I)=Y(I ) 102 C O N T I N U E  101  CALL F A S T F ( X , Y , N ) CALL F X F F T ( X 1 , Y 1 , X , Y , N ) CALL SMOOTH(XT.N) CALL SMOOTH(Y1,N) CALL RAD(W,N,DELT) CALL NORME(N,X1) CALL NORME(N,Y1 ) NN=N/2 DO 101 1=1,NN WRITE! 7 , * ) W ( I ) , X 1 ( I ) W R I T E 1 8 , * ) W ( I ) .Y 1 ( I ) CONTINUE CALL  PL(X1,W,NN)  WRITE(6,4) 4 FORMAK'What f r e q u e n c i e s READ(5,*)FC1,FC2  CALL  do y o u want  PL(Y1,W,NN)  WRITE(6,4) READ(5,")FC3,FC4  C A L L B N D P ' S ( F C 1 , F C 2 , D L T , D 1 ,G1 ) C A L L F I L !. , ( X X , N 1 , D 1 ,G1 , I G 1 ) CALL BNDPAS(FC3,FC4,DLT,D2,G2) CALL FILTER(YY,N1,D2,G2 , IG2) . CALL NORME(N1,XX)  167  filtered?')  Listing 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 1 10 111 112 11 3 114 1 15 1 16  o f CROSSCOR  a t 16:44:17  CALL  c  201  on JUL 27.  1988 f o r C C i d = S E I S on G  N0RME(N1,YY)  DO 2 0 1 I = 1 , N 1 WRITEO,*) X X ( I ) WRITEMC,') Y Y ( I ) CONTINUE  C  c  LG=2*N1-1 CALL  CROSS(N1,XX,N1,YY,LG,G)  C L A = 2 * LG DO 30 J = 1 . L G IJ =J- 1 GA(J+LG)=G(J) 30 CONTINUE C CALL CR0SS(N1,YY,N1.XX.LG.G) LH=LG+1 DO 4 0 1 1 = 1 , L G GA(II)=G(LH-II) 40 CONTINUE C  c  c c  c  c c  DO 5 0 1 = 1 . L G EX(I)=I-LG 50 CONTINUE LGG=LG+ 1 DO 5 3 I = L G G , L A EX(I)=I-LG 53 CONTINUE DO 1 0 9 1 = 1 , L A W R I T E ( 1 1 ,*) E X ( I ) , G A ( I ) 109 C O N T I N U E CALL  MAXSN(LA,GA,CMAX,LMAX)  HMAX=LMAX-LG- 1  c c c c c c c c c  CALL  PLA(GA,EX,LA)  W R I T E ( 6 , 3 0 0 ) HMAX 3 0 0 FORMAT ( ' T h e M a x i m u m  Time  W R I T E ( 6 , 3 0 3 ) FN I F NOT 303 FORMAK'CORRECT  Shift  is  :  ' ,17)  F I L T E R E D PAST NYQUIST  Fn =  STOP END SUBROUTINE F A S T F ( F R . F I , N ) N i s t h e number o f d a t a points=2 FR i s t h e r e a l d a t a s e t FI i s the imaginary p a r t o f d a t a set(=0.0  168  i f only  Listing 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 17 1 172 173 174  o f CROSSCOR a t 1 6 : 4 4 : 1 7 C C  First  compute  on JUL 2 7 , 1988 f o r C C i d = S E I S on G M  REAL F R ( N ) , F I ( N ) , G R , G I , E R , E I . E U . E Z M=0 KD=N 1 KD=KD/2 M=M+1 I F ( K D . G E . 2 ) GO TO 1 ND2 = N/2 NM1 = N-1 L=1 C C C  Shuffle  2 3  4 C C c  data  in binary  digit  reverse  order  DO 4 K=1,NM1 I F ( K . G E . L ) GO TO 2 GR = F R ( L ) GI = F I ( L ) FR(L) = FR(K) F K L ) = F K K ) F R ( K ) = GR F K K ) = GI NND2 = ND2 I F ( N N D 2 . G E . L ) GO TO 4 L = L-NND2 NND2 = N N D 2 / 2 GO TO 3 L = L + NND2 P I = 3 . 14 1 5 9 2 6 5 First DO 6 NJ = NJD2 EU = EZ = ER = EI =  C C C  input  arrange  accounting  of M  stage  J= 1 , M 2"* J = NJ/2 1 .0 0.0 COS(-PI/NJD2) SINK - P I / N J D 2 )  Compute  Fourier  transform  i n each  M  stage  DO 6 I T = 1 , N J D 2 DO 5 I W = I T . N . N J I W J = IW + N J D 2 GR = F R ( I W J ) * E U - F I ( I W J ) " E Z GI = F I ( I W J ) * E U + F R ( I W J ) * E Z F R ( I W J ) = F R ( I W ) - GR F I ( I W J ) = F I ( I W ) - GI F R ( I W ) = F R ( I W ) + GR 5 F I ( I W ) = F K I W ) + GI S E U = EU EU = S E U ' E R - E Z ' E I 6 E Z = E Z * ER + S E U ' E I RETURN END C C  FXFFT  corrects  the amplitudes  169  o f FASTF,  as o u t l i n e d  by K a n a s e w i c h  Listing 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232  o f CROSSCOR  at  16:44:17 on  JUL 27,  1988 f o r C C i d = S E I S  on  G  C SUBROUTINE F X F F T ( X 1 , Y 1 , X , Y , N ) REAL X(N),Y(N),X1(N),Y1(N) ,A(3000) ,B(3000) ,C(3000) REAL D ( 3 0 0 0 ) , E ( 3 0 0 0 ) ,F(3000) ,G(3000) ,H(3000) DO 10 1 = 1 ,N A(I) = .5*(X(I)+X(N+2-I ) ) B(I)=.5*(Y(I)-Y(N+2-I)) C(I)=A(I)'A(I) D ( I ) = B ( I ) *B( I ) X1(I)=SQRT(C(I)+D(I)) E(I)=.5*(X(I)-X(N+2-I)) F(I)=.5*(Y(I)+Y(N+2-I)) G(I)=E(I)*E(I) H(I)=F(I)"F(I) Y1(I)=SQRT(G(I)+H(I) ) 10 C O N T I N U E RETURN END C C C C  P0W2 d e t e r m i n e s M i n N=2**M f o r a s p e c i f i c t r a c e l e n g t h ( i . e , L ) , a n d t h e n p a d s t h e r e s t o f t h e t r a c e , L t o N, w i t h z e r o s . SUBROUTINE POW2(L,X,Y,N) DIMENSION X ( L ) , Y ( L ) REAL TT,ST,SS,S SS=0.30103 TT=FL0AT(L) ST=L0G10(TT) S=ST/SS+1.0 Q=IFIX(S) N=1 DO 2 0 1 = 1 ,Q N=2'N 20 CONTINUE DO 10 I = L , N X(I)=0.0 Y(I)=0.0 10 C O N T I N U E RETURN END  C C C  .  C  NORME  i s normi1ization  o f an a r r a y  SUBROUTINE NORME(LL,F) REAL F(4000),FS,2 INTEGER FF I F ( L L ) 30,30,10 CALL M A X S N t L L , F , F S , F F ) WRITE(6,*) FS Z=FS DO 2 0 1 = 1 , L L 20 F ( I ) = F ( I ) / Z RETURN END  C C DOT  I S DOT PRODUCT S U B R O U T I N E DOT ( L , X , Y , A N S ) I M P L I C I T REAL'8 (A-H.O-Z)  170  by  it's  RMS  energy.  Listing 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 26 1 262 263 264 265 266 267 268 269 270 27 1 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290  o f CROSSCOR  10 20 30 C C C  at  16:44:17  on  JUL  27,  1988  f o r C C i d = S E I S on  G  R E A L * 8 X ( 2 ) , Y ( 2 ) ,ANS ANS=0.0 IF(L) 30,30,10 DO 2 0 1 = 1 ,L ANS=ANSi-X(I)'Y(I) RETURN END REMAV  removes  the arithmetic  average  i n o r d e r DC  shifts  a r e removed  SUBROUTINE REMAV(LY,Y,AVERG) DIMENSION Y ( 2 ) S=0. DO 10 1 = 1 , L Y 10 S = S + Y ( I ) AVERG=S/FLOAT(LY) DO 2 0 1 = 1 ,LY 20 Y(I)=Y(I)-AVERG RETURN END C  c  c c c c c c c c c  SUBROUTINE RAD(W,N,DELT) D I M E N S I O N W(N) REAL SS.ST DO 10 1=1,N SS=FL0AT(I) ST=FLOAT(N) W(I)=SS/(ST'DELT) 10 C O N T I N U E RETURN END Smooth computes smoothed s p e c t r u m u s i n g Turkey+Hamming f o r m u l a . The SPECT=unsmoothed spectrum LS=its lenght The s u b r o u t i n e o u t p u t i s SPECT=smoothed spectrum SUBROUTINE  SMOOTH(SPECT,LS)  DIMENSION S P E C T ( L S ) MM=LS- 1 A=.54*SPECT(1)+.46*SPECT(2) B=.54*SPECT(LS)+.46*SPECT(MM) SJ=SPECT( 1) SK=SPECT(2) DO 10 J=2,MM SI=SJ SJ=SK SK=SPECT(J+1) SPECT(J)=.54*SJ+.23'(SI+SK) 10 C O N T I N U E SPECT(1)=A SPECT(LS)=B RETURN END  171  from the c o s i n e t r a n s f o r m subroutine inputs are  Listing 29 1 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 31 1 312 313 314 315 316 317 318 3 19 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 34 1 342 343 344 345 346 347 348  o f CROSSCOR  at  16:44:17 on  JUL 27,  1988 f o r C C i d = S E I S  on  C SUBROUTINE P L ( A M P , Z , L ) REAL AMP(3000),Z(3000),SS(3000),TT(3000) INTEGER L.FSMPL 1 N=C DO 1 0 0 I - 1 . L  s s m = z m TT(I)=AMP(I) I F ( I . E Q . L ) GO 1 0 0 N=I C C C C  IF  EQUAL  TO  ZERO  TO  200  TERMINATE  PROGRAM  200  I F ( N . E O . O ) GO TO 3 0 0 C A L L S U B R O U T I N E TO DO THE P L O T T I N G CALL SUMFUN(SS,TT,L) C AND R E T U R N FOR MORE DATA C P R O G R A M COMES H E R E WHEN F I N I S H E D . 300 CONTINUE C T E R M I N A T I N G P L O T T I N G THEN S T O P CALL PLOTND RETURN END C SUBROUTINE SUMFUN(X,Y,N) C S U B R O U T I N E DRAWS A L I N E THROUGH N P O I N T S REAL X(N),Y(N) C REAL A M P S { 3 0 0 0 ) , Z S ( 3 0 0 0 ) C  C C  AMPMIN=AMP(1)  C C C C C  DO 10 1=1,N AMPS(I)=SNGL(AMP(I)) ZS(I)=SNGL(2(1)) AMPMIN=AMIN1(AMPMIN,AMPS(I)) 10 C O N T I N U E c F I R S T S C A L E THE P O I N T S CALL SCALEtX,N,10.0,XMIN.DX, 1 ) CALL SCALE!Y,N,10.O.YMIN.DY,1) C A L L A X I S t O . , 0 . , ' F R E Q U E N C Y ' , - 9 , 1 0 . 0,0. .XMIN.DX) C A L L A X I S I O . ,0. . 'AMPLITUDE' ,9, 10.0 ,90. ,YMIN,DY) C F I N A L L C P L O T THE L I N E S C A L L L I N E U ,Y ,N,1) c MOVE THE O R I G I N AND R E T U R N RETURN END  c  SUBROUTINE P L A ( A M P , Z , L ) REAL A M P ( 3 0 0 0 ) , Z ( 3 0 0 0 ) , S S ( 3 0 0 0 ) ,T T ( 3 0 0 0 ) INTEGER L,FSMPL 1 N=0 DO 100 1 = 1 , L SS(I)=Z(I) TT(I)=AMP(I) I F ( I . E Q . L ) GO TO 2 0 0 1 0 0 N=I  172  G  Listing 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 3/8 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406  o f CROSSCOR C C C C  IF  at  EQUAL  16:44:17 on  TO  JUL  27,  ZERO T E R M I N A T E  1988 f o r C C i d = S E I S  on  G  PROGRAM  200  I F ( N . E Q . O ) GO TO 3 0 0 C A L L S U B R O U T I N E TO DO THE P L O T T I N G CALL SUMFN1(SS,TT,L) C AND R E T U R N FOR MORE DATA C P R O G R A M COMES H E R E WHEN F I N I S H E D . 300 CONTINUE C T E R M I N A T I N G P L O T T I N G THEN S T O P CALL PLOTND RETURN END C SUBROUTINE SUMFN1(X.Y.N) C S U B R O U T I N E DRAWS A L I N E THROUGH N P O I N T S REAL X(N),Y(N) C REAL A M P S ( 3 0 0 0 ) , Z S ( 3 0 0 0 ) C  C C  '  AMPMIN=AMP(1)  DO 10 1=1,N c C AMPS(I)=SNGL(AMP(I)) C ZS(I)=SNGL(Z(D) C AMPMIN=AMIN1(AMPMIN,AMPS(I)) C 10 C O N T I N U E C F I R S T S C A L E THE P O I N T S CALL SCALE(X,N,10.0,XMIN,DX,1) CALL SCALE(Y,N,10.0,YMIN,DY,1) C A L L AX I S ( 0 . , 0 . , 'TIME S H I F T ' , - 9, 1 0 . 0 , 0 . , X M I N , D X ) C A L L A X I S 1 0 . , 0 . , 'CROSS C O R R E L A T I O N V A L U E ' , 9, 10 . 0 , 9 0 . , Y M I N , D Y ) C F I N A L L C P L O T THE L I N E S CALL L I N E ( X , Y,N,1) C MOVE THE O R I G I N AND R E T U R N RETURN END  c c c c c c c c c c c c c c  SUBROUTINE  BNDPAS(F1,F2,DELT,D,G)  Subroutine  by Dave G a n l e y  on March  5,  1977.  The p u r p o s e o f t h i s s u b r o u t i n e i s t o d e s i g n a n d a p p l y a r e c u r s i v e B u t t e r w o r t h band pass f i l t e r . In order to d e s i g n t h e f i l t e r a c a l l m u s t b e made t o BNDPAS a n d t h e n t h e f i l t e r may b e a p p l i e d b y c a l l s t o f i l t e r . T h e f i l t e r w i l l have 8 p o l e s i n t h e s p l a n e and i s i n f o r w a r d and r e v e r s e d i r e c t i o n s so as t o have z e r o phase c u t o f f f r e q u e n c i e s . T h e a t t e n u a t i o n w i l l b e -6DB a n d t h e r o l l o f f w i l l b e a b o u t 9 6 DB p e r O c t a v e . A B i l l i n e a r Z t r a n s f o r m i s used i n d e s i g n i n g the f i l t e r t o prevent a l i a s i n g problems. COMPLEX P ( 4 ) , S ( 8 ) , Z 1 , Z 2 DIMENSION D ( 8 ) , X ( 2 0 4 8 ) , X C ( 3 ) , X D ( 3 ) , X E ( 3 ) DATA I S W / 0 / , T W O P I / 6 . 2 8 3 1 8 5 3 /  c c  This  section  calculates  173  the f i l t e r  and must  be  called  Listing 407 408 409 410 411 412 413 4 14 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464  o f CROSSCOR C C C C C C C C  C C  C C C C C C  at  16:44:17 on  before F1 = F2 = DELT D = W G = W  filter  JUL  27,  1988  f o r C C i d = S E I S on  G  is called.  L o w f r e q u e n c y c u t o f f ( 6 DB D o w n ) H i g h f r e q u e n c y c u t o f f ( 6 DB d o w n ) = Sample i n t e r v a l i n m i l l i s e c o n d s i l l c o n t a i n 8 Z domain c o e f f i c i e n t s of r e c u r s i v e i l l contain the gain of the f i l t e r  filter  WRITE(6,1) F1,F2,DELT 1 F O R M A T ! ' B A N D P A S S F I L T E R D E S I G N FOR A BAND F R O M ' , F 8 . 3 , ' T O ' , F 8 1.3,'HERTZ.',11' SAMPLE I N T E R V A L I S ' , F 5 . 2 , ' M I L L I S E C O N D S . ' ) DT = D E L T / 1 0 0 0 . 0 TDT = 2.0/DT FDT = 4.0/DT ISW = 1 P(1) = CMPLX!-.3826834,.9238795) P(2) = CMPLX!-.3826834,-.9238795) P(3) = CMPLX!-.92388795,.3826834) P(4) = CMPLX(-.99238795,-.382684) W1 = T W O P I ' F I W2 = T W O P I * F 2 W1 = TDT * T A N ( W 1 / T D T ) W2 = T D T * T A N ( W 2 / T D T ) HWID = ( W 2 - W U / 2 . 0 WW = W1'W2 DO 19 1=1,4 Z1 = P ( I ) ' H W I D Z 2 = Z1'Z1-WW Z2 = C S Q R T ( Z 2 ) S ( I ) = Z1 + Z 2 19 S ( I + 4 ) = Z 1 - Z 2 WRITE(6,2) S 2 F O R M A T ( ' - S P L A N E P O L E S ARE A T : ' ,/' ' , 8 ( / ' ' , E 1 2 . 6 , ' +I ' ,E 1 2 . 6 ) ) G = .5/HWID G = G"G G = G" G DO 29 1 = 1 , 7 . 2 B = -2.0*REAL(S(I) ) Z1 = S ( I ) * S ( I + 1 ) C = REAL(Z1) A = TDT+B+C/TDT G = G*A D ( I ) = (C'DT-FDTl/A 29 D( 1 + 1 ) = ( A - 2 . 0 ' B ) / A G = G*G WRITEC6.3) 3 FORMAT!'-FILTER IS (1-Z*'2)"*4/B1*B2*B3*B4') WRITE(6,4) D 4 F 0 R M A T ( 4 ( / ' B ( I ) = 1 + ' . E 1 2 . 6 , ' Z+' . E 1 2 . 6 , ' Z " 2 ' ) ) WRITE(6,5) G 5 FORMAT!'-FILTER GAIN I S ' . E 1 2 . 6 ) RETURN  C C ENTRY F I L T E R ( X , N , D , G , I G ) C C C  X = D a t a v e c t o r o f l e n g t h N c o n t a i n i n g d a t a t o be D = F i l t e r c o e f f i c i e n t s c a l c u l a t e d by BNDPAS  174  filtered  Listing 465 466 467 468 469 470 47 1 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 49 1 492 493 494 495 496 497 498 499 500 501 502 50 3 504 505 506 507 508 509 510 51 1 512 513 514 515 516 517 518 519 520 521 522  o f CROSSCOR C C C C  at  16:44:17  on  JUL  27,  G = F iIter Gain IG = 1 M e a n s t o r e m o v e is unity. I F ( I S W . E Q . 1 ) GO TO 31 WRITE(6,6) 6 F O R M A T f ' 1 B N D P A S MUST B E CALL EXIT  1988  f o r CCid=SEIS  the f i l t e r  CALLED  gain  BEFORE  so  W R I T E ( 6 , " ) X( 1 ) Apply f i l t e r in forward 31  direction  XM2=X(1) XM1=X(2) XM=X(3) XC(1)=XM2 XC(2)=XM1-D(1)*XC(1) XC(3)=XM-XM2-D(1)*XC(2)-D(2)*XC(1) XD(1)=XC(1) X D ( 2 ) = X C ( 2 ) - D ( 3 ) ' X D ( 1) XD(3)=XC(3)-XC(1)-D(3)'XD(2)-D(4)*XD(1) XE(1)=XD(1) XE(2)=XD(2)-D(5)*XE( 1) XE(3)=XD(3)-XD(1)-D(5)*XE(2)-D(6)'XE(1) X ( 1 ) = X E ( 1) X ( 2 ) = X E ( 2 ) - D ( 7 ) * X ( 1) X(3)=XE(3)-XE(1)-D(7)*X(2)-D(8) X(1) DO 39 1=4,N XM2=XM1 XM1=XM XM=X(I) K = I - ( ( I - 1 ) / 3 ) *3 GO TO ( 3 4 , 3 5 , 3 6 ) ,K M=1 M1 = 3 M2=2 GO TO 37 M=2 M1 = 1 M2 = 3 GO TO 37 M=3 M1 = 2 M2=1 XC(M)=XM-XM2-D(1)"XC(M1)-D(2)'XC(M2) XD(M)=XC(M)-XCIM2)-D(3)*XD(M1)-D(4)*XD(M2) XE(M)=X0(M)-XD(M2)-D(5)'XE(M1)-D(6)"XE(M2) X(I)=XE(M)-XE(M2)-D(7)'X(I-1)-D(8)'X(I-2) ,  34  35  36  37  39 C  c  Filter  in reverse  direction  C XM2=X(N) XM1=X(N-1) XM=X(N-2) X C ( 1 )=XM2 XC(2)=XN1-D(1)"XC(1) XC(3)=XM-XM2-D(1)'XC( XD(1)=XC(1)  175  2 ) - D ( 2 ) ' XC( 1)  that  FILTER')  C C C  on  G  Listing 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 54 1 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 56 1 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580  o f CROSSCOR  a t 16:44:17 on JUL 2 7 , 1988 f o r C C i d = S E I S  on G  XD(2)=XC(2)-D(3)*XD(1) XD(3)=XC(3)-XC(1)-D(3)'XD(2)-D(4)"XD(1) XE(1)=XD( 1 ) X E ( 2 ) = X D ( 2 ) - D ( 5 ) ' X E ( 1) XE(3)=XD(3)-XD(1)-D(5)'XE(2)-D(6)*XE(1) X ( N ) = X E ( 1) X ( N - 1 ) = X E ( 2 ) - D ( 7 ) 'X( 1 ) X(N-2)=XE(3)-XE(1)-D(7) X(2)-D(8)'X(1) DO 4 9 1 = 4 , N XM2=XM1 XM1=XM J=N-I+1 XM=X(J) K=I - ( ( I - 1)/3)* 3 GO TO ( 4 4 , 4 5 , 4 6 ) ,K 44 M=1 M1 = 3 M2=2 GO TO 47 4 5 M-2 M1 = 1 M2=3 GO TO 47 46 M=3 M1=2 M2=1 47 X C ( M ) = X M - X M 2 - D ( 1 ) ' X C ( M 1 ) - D ( 2 ) * X C ( M 2 ) XD(M)=XC(M)-XC(M2)-D(3)*XD(M1)-D(4)*XD(M2) XE(M)=XD(M)-XD(M2)-D(5)"XE(M1)-D(6)* XE(M2) 49 X ( J ) = X E ( M ) - X E ( M 2 ) - D ( 7 ) * X ( J + 1 ) - D ( 8 ) * X ( J + 2 ) I F ( I G . N E . I ) RETURN DO 5 9 1 = 1 ,N 59 X ( I ) = X ( I ) / G W R I T E ( 6 , " ) ( X ( I ) , 1=1,N) RETURN END ,  C C C  Subroutine  DOT  SUBROUTINE  finds  the dot product  o f two a r r a y s .  DOT(L,X,Y,ANS)  C I M P L I C I T REAL (A-H.O-Z) REAL X(2),Y(2),ANS A N S = 0 .0 IF(L) 30,30,10 10 DO 2 0 1 = 1 ,L 2 0 A N S = A N S + X( I ) * Y( I ) 30 R E T U R N END C C C  Subroutine  CROSS  SUBROUTINE  finds  the cross  product  o f two a r r ,  CROSS(LX,X,LY,Y,LD,D)  C I M P L I C I T REAL (A-H.O-Z) REAL X(2),Y(2),D(2) DO 10 1 = 1 , LD D0T(MIN0(LY+I-1,LX)-I+1,X(I),Y,D(I)) 10 C A L L  176  Listing 581 582 583 584 585 586 587 588 589 590 59 1 592 593 594 595 596 597 598  o f CROSSCOR  at  16:44:17 on  JUL 27,  1988 f o r C C i d = S E I S  on G  RETURN END C C C C  S u b r o u t i n e MAXSN f i n d s t h e m a x i m u m taking into account the a l g e r b r a i c SUBROUTINE  value signs  MAXSN(LX,X,XMAX,INDEX)  C REAL X(2),XMAX XMAX=X(1) DO 10 1 = 1,LX 10 X M A X = A M A X 1 ( X M A X , X ( I ) ) DO 2 0 J = 1 , L X INDEX=J IF(X(J)-XMAX) 20,30,20 20 CONTINUE 30 R E T U R N END  177  o f an a r r a y , of the elements.  APPENDIX  (3  Program l i s t i n g of KALDAT, the d i s c r e t e model of the second order system.  178  L i s t i n g of KALDAT at 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 59.5  C C C C C C  C C C C C C C C C C C C C C C c  c c c  c  18:25:29 on DEC 12. 1987 for  CCId-SEIS  This program generates s y n t h e t i c geophone data Inorder to determine the performance of the Kalman f i l t e r . The Input n o i s e , wk. is white, and has a gausslan d i s t r i b u t i o n with mean=0.0 and standard d e v l a t 1 o n « S I G M A REAL REAL REAL REAL REAL REAL INTEGER  A 11.A12.A21.A22.B11.821.WO.R0.PO.R01.P01 TC.GAMA.HO.HI,H2,KI,WSIGMA,WMEAN,WS.ALPH,DELT DM,012,VS,AW, BW, OA ,SSS,SA.SS.BB11,BB1,ST,SH,SP Xl(40OO).X2(40OO).X3(4OOO).X4(40OO).X5(400O) X6(4000),X7(4000).X8(4000).X9(4000).2(4000) X81(40O0),X82(40O0).X91(4000).XS2(4000).V(4000) T  PIE=3.1415927 length of d e s i r e d time s e r i e s T = 200 sampling Interva1.de1ta DELT=.01 standard d e v i a t i o n of noise present WSIGMA=.09 VSIGMA=.009 natural frequency of WO=8»2»PIE damping of R0=.IB  Insrument  Instrument  KI f a c t o r . KI=1.0 gain on measurement H HO=.29 time constant T O . 10 DO IOO I » 1 . T READ( 15.•) X1(I),X3(I ) 100 CONTINUE WMEAN-0.0 GAMA- 1/TC ALPH*ASIN(RO) R01=R0»W0*DELT P0«S0RT(1-R0«R0) POI-PO'WO'DELT AA = EXP(-R01 )/P0 A11=AA«(C0S(P01-ALPH)) A !2-(AA*SIN(P01))/WO A21»-AA»(W0»SIN(P01)) A22"AA*(COS(PO1+ALPH)) B 11-( 1 0 - A 1 1 ) » K I / ( W 0 * W 0 )  179  L i s t i n g of KALOAT at 18:25:29 on DEC 12. 1987 f o r CCid=SEIS 59.7 60 61 62 63 64 65 66 67 68 70 72 73 74 75 76 77 78 79 80 81 82 83 84 85 85 . 5 86 87 88 89 90 91 92 93 94 9596 97 98 99 10O 101 102 105 106 107 108 1 1 1 112 113 1 14 1 15 1 16 1 17 1 18 1 19 120 12 1  B12=-A21'KI/(WO*WO) H1«-W0*W0 H2=>-2»R0«WO QA=GAMA*DELT AW=£XP(-OA) BW=1.0 -AW X5(1)=0.0 X6(1)=0.0 X7(l)=0.0 X8(l)=0.0 X9(1)=0.0 2(1)=0.0 WRITE(6.') AW.BW C C  00  1 I- 1 .T  c  At t h i s p o i n t i t i s necessary to generate the G a u s s i a n n o i s e C as input input our f i r s t order transform funct1on(which r e p r e s e n t s c a low pass f i l t e r ) which has X5 as our o u t p u t . X5 i s then f e e d c i n t o our second order model.  c  c c c c c c c c c c  TK=I ST=1.0«1 SH=2.0'I SS=SQRT(D£LT) SP=GAMA/2.0 CALL RNGAUS(WS.WMEAN,WSIGMA.TK,ST) NN=I+1 X5(NN)=X5(I)»AW+BW»WS X6(NN ) = B 11 *X5(I) + A 11 *X6(I) + A12*X7(I) X7(NN)=B12*X5(I ) + A2 1*X6(I ) + A 2 2 » X 7 ( I ) BB 1t*Bl1*X1(I) BB1=B1 1'X3(I ) X8(NN)=A11*X8(I ) + A12*X9(I)+BB11+BB1  X9(NN) A21*X8(I)+A22*X9(I)+B12*Xl(I)+B12»X3(I) i  c c c c c c  CALL RNGAUS(VS.WMEAN.VSIGMA.TK.SH) Z(NN)=H0«(X7(NN)+X9(NN))+VS/SS WRITE(2. • ) XS(I ) 1 CONTINUE DO 2 J=1.T  180  Listing 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168  KALDAT at  C C C C C C C C  C C C C C C C C C  18:25:29 on DEC 12. 1987 for  WRITE(1.•) 2 CONTINUE  CCId-SEIS  Z(J)  DO 2 J - l . T W R I T E ( 2 . » ) X5(J)+1 .4 2 CONTINUE STOP END SUBROUT INE RNGAUS ( W. W.ME AN. WS IGMA . TK . SS ) REAL INTEGER DATA DATA  RNSEO,WMEAN,C1,C2.C3,WSIGMA,W,SS TK IRNSED/1/ TW0P1/6.283185307/  IF  (TK.LT.10.0) THEN SS-TK*1000.0+1.0 ELSE SS-TK»100.0+1 END IF GO TO (10.20.30)  IRNSED  10 RNSED-ABS(RNSED) CI-RAND(RNSED)  C C F i r s t c a l l for a random number. C C1=S0RT(-2.*AL0G(RAN0(O.))) 20 C2-TW0PI»RAND(0.) C3'C1»C0S(C2) IRNSED-3 GO TO 40 C C Second c a l l for random number. C C3-C1«SIN(C2) 30 IRNSED-2 C 40 W-WMEAN+WSIGMA-C3 C RETURN END  181  APPENDIX C  Program l i s t i n g of the Kalman F i l t e r , KALMAN 4  182  Listing 1 2  . KALMAN4 at  10:07:47 on DEC 13. 1987 for  CCId-SEIS  C C  3  £•»•*•*••*»•••*»»»•'****»***«*•*********•**••*•#««**•*••••**«••*••»#**  4 5 6 7 8 9 10 11 j2 13 14 15 16 16.5 16.7 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56  C* • C» MAIN ' PROGRAM * C* • C* • C* • C« • C* • C« . • £••••«••••»*•»•***••»•*•*****••*•*••»*•»•••*••••*•»**••»••»*•#«»•**•» C C C INTEGER SNP.SN.ST.SI.SP.NT.FT,KT.QT.T.IDEX4.IDEX2 INTEGER TE(4) REAL X2.X22.X4.X44 C REAL WO.RO.KI.BB.LL.RR.R.0ELT.A11.A21.A12.A22 REAL AA.ALPH.ESP.B1t,B12.0R.GAMA.C1.C2.C3.C4,P01.R01 REAL PIE.AN.RW.HH.TC,AW.QW.RA.PW.MV,HV(9) REAL XMIN1.XMIN2,XMIN3.XMIN4.XMAX4.XMAX2.SNN.SPP REAL PH1(500).PH2(500).PH3(500),PH4(500),XH1(500).XH2(SOO) REAL XH3(500).XH4(50O).KK.QD.N1.NPP.NSS.SD.N3 INTEGER INOEX1.INDEX2.INDEX3.IN0EX4,N2.N4,NP(6).NS(6) C REAL Z(2000).WP(200).WS(200).BB!1.BB21.BB1.BB2 REAL KS(9).KH(9.9).BETA(10).QN(2.2).PS0(9.9) REAL P S ( 9 , 9 ) . P S 1 ( 4 . 4 ) . P S 2 ( 4 . 5 ) . P S 2 T ( 5 . 4 ) . P S 3 ( 5 . 5 ) REAL P U ( 9 . 9 ) . P U 1 ( 4 . 4 ) . P U 2 ( 4 . 5 ) , P U 2 T ( 5 . 4 ) , P U 3 ( 5 , 5 ) REAL FS(9.9).FS1(5.4),FS2(5.5),XU(9),XE(9).X0(9) REAL F L ( 9 . 9 ) . F L l ( 5 . 4 ) . F L 2 ( 5 . 5 ) . H 1 1 ( 2 0 0 ) . H 1 2 ( 2 0 0 ) C C C D e f i n i n g parameters: C C WO=natural frequency of Instrument R0=damplng of Instrument C KI= c o e f f 1 c I e n t which takes into account instrument inductance. C r e s i s t a n c e and magnetic f i e l d . C 0ELT=samp1tng Interval AN=white noise Input into second order C system. C ER=covariance on measurement noise QR=covarlance on s t a t e noise C C i n i t i a l estimates of the s t a t e s C REAO(11.*)X0(1),X0(2),X0(3),X0(4).X0(5).X0(6).X0(7).X0(8) 1 .X0(9).FT C C DO 76 1=1.9 DO 77 J=1.9 PS0(I.J)=0.0 77 CONTINUE 76 CONTINUE C READ(11.•)PS0(1.1).PS0(2.2).PSO(3.3).PS0(4,4).PSO(5.5). 1PS0(6.6),PSO(7.7).PSO(8.8).PS0(9.9).PSO(6,7).PS0(8.9).  183  L i s t i n g of KALMAN4 at 10:07:47 on DEC 13. 1987 f o r CC1d*SEIS 57 58 59 eo 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 8t B2 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 too 101 102 103 104 105 106 107 108 109 1 10 1 1 1 1 12 113 1 14  1PS0( 1 . 8 ) . P S 0 ( 1 . 9 ) . P S 0 ( 2 . B ) . P S 0 ( 2 1PS0(3.8).PS0(3.9).PS0(4.8).PS0(4 C C  PS0(8.1)-PS0(1.8) PS0(9.1)-PS0(1.9) PS0(8.3)»PS0(3.8) PS0(9.3)-PS0(3.9) PS0(8.4)»PS0(4.8) PS0(9.4)=PS0(4.9) PS0(8.2)-PS0(2.8) PS0(9.2)-PS0(2.9) PS0(7.6)«PS0(6.7) PS0(9.8)=PS0(8.9)  C C C C  DELT-.01 KI-1.0 T = 200 PIE=3.1415927  C W0=8»2*PIE C R0 =.18 C C NP: C C C C NS:  length of NS( 1) = 1 NP(2)=2 NP(3)=3 NSS=1  Input P-wave  length of Input S-wava DO 232 1 = 1. 10 NP(I)=I C NS(3)=3 232 CONTINUE NPP=9 c c WP: c o e f f i c i e n t s d e f i n i n g P-wave WS( 1 )=0.0 c WP(2)=1.0 c WP(3)=0.0 c WP(4)=0.0 c WP(5)-.6 c WP(6)=.8 c c WS: c o e f f i c i e n t s d e f i n i n g S-wave 00 323 1 = 1. 10 F1=50.0 QAS= .01 *(I ) WP(I)=C0S(F1»QAS)*EXP(-QAS) 323 CONTINUE c WS(2)=1.0 c WS(3)«0.0 c WS(4)=0.0 . c WS(5)=.6  184  L i s t i n g d i KALMAN4 at 10:07:47 on DEC 13. 1987 f o r  I  I  ! |  115 116 1 17 1 18 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 163.5 163.7 163.75 163.77 163.8 163.9 164 165 166  C WS(6)=8 C C C MV: g a i n or a m p l i f i c a t i o n on recorder MV=.29 C C FT: number of I t e r a t i o n s d e s i r e d C FT = 4 C 00 78 1=1.T R E A D ( 8 . « ) Z(I) 78 CONTINUE C C TC: time constant TC=.10 C C ESP=1/TC H1=-W0'W0*1.0 KI=1.0 H2=-2»RO»WO*1.0 P0=S0RT(1-R0»R0) R01«R0»WO*DELT P01=P0'WO»DELT GAMA=ATAN(RO/PO) ALPH=ASIN(RO) AA=EXP(-R01)/P0 A12«(AA'SIN(P01))/W0 A21=-AA«(W0»SIN(P01)) A22=AA»(C0S(P01+ALPH)) A11=AA«(C0S(P01-ALPH)) B1 1 = ( 1 - A 1 1 ) » K I / ( W O * W O ) B12=-A21«KI/(WO*WO) C 1 = ((VO*WO)»EXP(-R01))/P0 C2=SIN(P01)/W0 C3=C0S(P0t+ALPH) C4=P0»EXP(-R01) AW=EXP(-ESP'DELT) BW=1.0-AW 0R=.O9 0W=BW*BW»0R RW=.009 C DO 9 1=1.9 XE(I)=X0(I) 9 CONTINUE C DO 1 111=1.FT C TEC 1 ) = IFiX(XE( 1 )) TE(3) = IFIX(XE(3) ) XE( 1 ) = FLOAT(TE( 1 ) ) XE(3)=FL0AT(TE(3) ) C C X2=XE(2)/OELT TE(2)=IFIX(X2) X22=TE(2)'DELT v  185  CCId-SEIS  L i s t i n g of KALMAN4 at 167 167.5 168 169 170 171 171.5 171.6 171.7 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220  10:07:47 on DEC 13. 1987 f o r C C i d - S E I S  XE(2)-X22 C X4-XE(4)/DELT TE(4)-IFIX(X4) X44-TE(4)'DELT XE(4)»X44 C C C  WRITE(6.»)  XE(1).XE(2 ) .XE(3).XE(4)  XE(5)=X0(5) XE(6)=X0(6) XE(7)=XO(7) XE(8)«X0(8) .Xc(9)»X0(9) C DO 70 1 - 1 . 9 DO 80 J - 1 . 9 PS(I.J)=PS0(I.J) 80 CONTINUE 70 CONTINUE C DO 30 1=1.4 DO 31 J - 1 . 4 PS1(I.J)=PS(I.J) 31 CONTINUE 30 CONTINUE C DO 32 1=5.9 DO 33 J - 5 . 9 D = d-4 D0=I-4 PS3(DD,D)=PS(I.J) 33 CONTINUE 32 CONTINUE C DO 34 1=1.4 DO 35 d - 5 . 9 D=J-4 PS2( I . D ) » P S ( I .'J) 35 CONTINUE 34 CONTINUE C C DO 17 1=1.9 00 18 0 - 1 . 9 FS(I.J)=0.0 FL(I.J)=0.0 18 CONTINUE 17 CONTINUE - C C FS(5.5)=AW FS(6.5)=811 FS(6.6)=A11 FS(6.7)-A12 FS(7.5)=B12 FS(7.6)=A21 FS(7.7)=A22  186  L i s t i n g u. KALMAN4 at 10:07.47 on DEC 13. 1987 f o r CC1d=SEIS 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 . 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278  FS(8.8)=A11 FS(8.9)=A12 FS(9.8)=A21 FS(9.9)=A22 C  C  C C C C C C C C C C C C C C C C  C C C C C  DO 21 1*1.5 DO 22 J - 1 . 5 FL(I.U)-FS(I.J) 22 CONTINUE 21 CONTINUE DO 23 1-1.4 FS<I.I )- 1 . 0 FL(I.I)=1.0 23 CONTINUE SI » 0 SP = 0  DO 2 11 = 1 .T KT=II Lines 328-356 d e f i n e the e r r o r t r a n s i t i o n matrix. BS, which has nonnutl components d e f i n e d by the d i s c r e t e second order system equat1ons(I.e. B1I.B12). T h i s matrix Is a l s o a component of the s t a t e t r a n s i t i o n m a t r i x , F, and is a f u n c t i o n of the P and S-wave time of a r r i v a 1 s ( t p - X 2 . t s - X 4 ) .  N1=NPP»DELT 0D=XE(2)+N1 KK=KT *DELT IF (KK.GE.XE(2) FS(8.a)»At1 FL(8.8)=A11 FS(8.9)=A12 FL(8.9)=A12 FS(9.8)=A21 FL(9.8)-A21 Fi(9.9)=A22 FL(9.9)=A22  OR.KK.GE.XE(4) ) THEN  ELSE GOTO 109 END IF 109 IF (KK.LT.XE(2).0R.KK.GT.0D) BB11=0.0 BB21=0.0  187  THEN  L i s t i n g of KALMAN4 at 10:07:47 on DEC 13, 1987 f o r 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 29B 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336  CCIOSEIS  GOTO 173 C ELSE C SI=SI+1 C IF (SI.GT.NPP+1) THEN GOTO 173 ELSE GOTO 174 END IF 174 N2=NP(SI) BB 1 1 *B 1 1*WP(N2) BB21=B12*WP(N2) C 173 END  IF  C C N3=NSS*DELT SD=XE(4)+N3 C IF  ( K K . L T . X E ( 4 ) . O R . K K . G T . S D ) THEN  C C  BB1=0.0 BB2=0.0 ELSE SP=SP+1 IF (SP.GT.NSS+1) THEN GOTO 171 ELSE GOTO 172 END IF 172 N4=NS(SP) BB1=B11*WS(N4) BB2»812«WS(N4)  C C  WRITE(6.»)  KK.XE(4).SD.WS(N4)  •  17 1 END  IF  C C C  WRITE(6.•)BB11,BB21.BB1,BB2  C C C  WRITE(6.*) F S ( 8 . 1 ) , F S ( 9 . 1 ) . F S ( 8 . 3 ) , F S ( 9 . 3 )  FS(8. 1 )°BB 1 1 FS(9.1)=BB21 FS(8.3)=BB1 FS(9.3)=BB2 FL(8. 1)«BB1 1 FL(9.1)=BB21 FL(8.3)=BB1 FL(9.3)=8B2  DO 303 1=1.9 HV(I)=0.0 303 CONTINUE  C  188  L i s t i n g o. KALMAN4 at 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 3B5 386 387 388 389 390 391 392 393 394  C C C C C C C C  10:07:47 on DEC 13. 1987 f o r  CCId-SEIS  HV(7)=MV HV(9)-MV CALL XUPDAT(XU.XE.FS) CALL EXTEND(FL.C1,C2.C3.C4.XU.KK.N1,N2.WS.WP.DELT.N3,N4.FT)  CALL EXTRAP(PU.PS.FL) PJ(5.5)»PU(5.5)+QW  C C C C C  C  C C C C C C C C C C C  CALL GAIN(HV.PU.RW.KS) CALL NWSTAT(XE.KS.XU.MV.KT.Z) CALL PUPDAT(HV.KS.PU.PS) P H 1 ( K T ) « P S ( 1 . 1) PH2(KT)*PS(2,2) PH3(KT)=PS(3.3) PH4(KT)=PS(4.4) XH1(KT)=XE(1) XH2(KT)=XE(2) XH3(KT)=XE(3) XH4(KT)=XE(4) WRITE(1.*) WRITE(2.«) W3ITEO.*) WRITE(4.«) WRITE(5,«) WRITE(7,«) WRITE(10.») WRITE(9.»)  PHI(KT) PH2(KT) PH3(KT) PH4(KT) XE(1) XE(2) XE(3) XE(4)  2 CONTINUE  CALL MINSN(T.PHI.XMINI.INOEXI) W R I T E ( 6 . « ) XMIN1.XH1(INDEX 1).INDEX 1 XE( 1 ) » X H 1 ( I N D E X 1 ) CALL MINSN(T.PH2.XMIN2.INDEX2) XE(2)=XH2(INDEX2) CALL MINSN(T.PH3.XMIN3.INDEX3) W R I T E ( 6 . « ) XMIN3,XH3(INDEX3).INDEX3 XE(3)* XH3(INDEX3) CALL MINSN(T.PH4.XMIN4.IN0EX4)  189  L i s t i n g of KALMAN4 at 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 41 1 412 4 13 4 14 4 15 4 16 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 44 1 442 443 444 445 446 447 448 449 450 451 452  C C C C C C C C C C C C C C C C C C C  c c c c c c c c c c c c c c c c c c c c c c c c c c c c c  10:07:47 on DEC 13. 1987 f o r  CCfd=SEIS  XE(4)»XH4(INDEX4) W R I T E ( 6 . » ) XMIN4.XE(4).INDEX* CALL MAXSN(T.PH2.XMAX2.I0EX2) WRITE(6.•) XMAX2.XH2(IDEX2+2).IDEX2 CALL MAXSN(T.PH4.XMAX4.I0EX4) W R I T E ( 6 . » ) XMAX4.XH4(IDEX4).IDEX4 IF  (INDEX2.E0.1)  THEN  XE(2)»XH2(IDEX2) ELSE XE(2)=(XH2(INDEX2+2)+XH2(INDEX2)+XH2(INDEX2+1))/3 END  IF  IF  (INDEX4.E0.1)  THEN  XE(4)=(XH4(IDEX4+2)+XH4(IDEX4)*XH4(IDEX4+1))/3 XE(4)=XH4(IDEX4) ELSE XE(4)»(XH4(INDEX4+2)+XH4(INDEX4)+XH4(INDEX4+1))/3 END  IF  SNN»XE(2)/DELT SN=IFIX(SNN) XE(2)=SN*0ELT W R I T E ( 6 . » ) XE(2) SPP=XE(4)/DELT SNP=IFIX(SPP) XE(4)=SNP*DELT WRITE(6.*) XE(4) WRITE(6.»)  XMIN2.XE(2).INDEX2  WRITE(6.»)  XMIN1.XE(1).INDEX 1  1 CONTINUE STOP END T h i s s u b r o u t i n e e x t r a p o l a t e s the s t a t e matrix SUBROUTINE REAL  XUPDAT(XU.XE.FS) XU(9),XE(9),FS(9.9)  190  X(I).  Listing 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 47 1 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 49 1 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510  KALMAN4 at C  10:07:47 on DEC 13. 1987 for CCid=>SEIS  INTEGER  OR  OR • 9 C C  CALL GMATV(FS.XE.XU.OR.OR.OR) RETURN ENO  C C C This subroutine propagates and updates the Kalman gain matrix K. C which is a f u n c t i o n of P and S-wave a r r i v a l s . C SUBROUTINE GAIN(HV.PU.RW.KS) • C REAL RW,PU(9,9),KS(9),HV(9).TMP1(9) REAL D0T.TMP2(9) INTEGER OR C OR = 9 C CALL GVMAT(HV.PU.TMP1.OR.OR.OR) C DOT » GVV(TMP1.HV.OR) • RW C W R I T E ( 6 . » ) DOT C CALL GMATV(PU.HV.TMP2.0R,0R.0R) C DO 1 1=1.9 KS(I)»TMP2(I)/DOT 1 CONTINUE  c  c c c c c  c c c c c c  RETURN END T h i s subroutine gives us s t a t e estimate updates. SUBROUTINE NWSTAT( XE,KS,XU.MV,KT,Z) REAL REAL INTEGER  DELTZ.Z(200O).XU(9).XE(9) KS(9).MV KT  DELTZ(KT)«Z(KT)-H(1)»X1(KT)-H(2)•XE(2)-H(3)»X3(KT)- H(4)»XE(4) 1-H(5)*X5(KT)-H(6)*X6(KT)-H(7)»X7(KT)-H(8)«X8(KT)-H(9)»X9(KT) DELTZ = Z(KT)-MV(XU(7)+XU(9) ) WRITE(6.»)  KS(1).KS(2).KS(3).KS(4).KT  XE( 1 )=XU(1)+KS(1)*DELTZ XE(2)=XU(2)+KS(2)»DELTZ XE(3)=XU(3)+KS(3)»0ELTZ XE(4)=XU(4)+KS(4)'0£LTZ X£(5)=XU(5)+KS(5)«DELTZ XE(6)=XU(6)+KS(6)«DELTZ XE(7)=XU(7)+KS(7)'DELTZ  191  I i  L i s t i n g of KALMAN4 at 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 54 1 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568  10:07:47 on DEC 13. 1987 f o r  CCld=SEIS  XE(8)«XU(8)+KS(8)»DELTZ XE(9)-XU(9)+KS(9)»DELTZ C RETURN END C C C C C C C  ^ S i n c e the t r a n s i t i o n m a t r i x , F, Is n o n l i n e a r w i t h r e s p e c t to the s t a t e parameters X2 and X4, i t Is necessary to l i n e a r i z e It. T h i s Is accomplished In the f o l l o w i n g s u b r o u t i n e , by a p p l y i n g L e i b n i t z ' s r u l e and the extended Kalman f i l t e r equations. SUBROUTINE EXTEND(FL,CI,C2.C3,C4,XU.KK.N1,N2.WS.WP.OELT,N3.N4.FT 1)  C INTEGER REAL REAL  MS.SL.SH,N2,N4.FT MM,NN.C1.C2.C3.C4.DELT,N1,N3,KK XU(9).WP(200).WS(200).FL(9.9)  C C C  MM=XU(2)+N1 WRITE(6.*)KK.XU(2).MM IF (KK.LT.XU(2).OR.KK.GT.MM) THEN  C C C  FL(8.2)=0.0 FL(9.2)-0.0 ELSE SH-N2+1  C C C C C C C C  FL(9,2)=C1*(C3»WP(N2)-C4»WP(SH))*XU(1) FL(8.2)=C1»C2*WP(N2)»XU(1) W R I T E ( 6 , « ) F L ( 9 , 2 ) . F L ( 8 . 2 ) . C 1 . C 2 , C 3 . W P ( N 2 ) . X U ( 1) END IF NN=XU(4)+N3 IF  (KK.LT.XU(4).0R.KK.GT.NN)  THEN  FL(8.4)=0.0 FL(9.4)=0.0 C ELSE C SL=N4+1 C FL(8.4)=C1»C2»WS(N4)*XU(3) FL(9.4)»C1»(C3'WS(N4)-C4»WS(SL))»XU(3) C WRITE(6.«) 104 END IF  FL(8.4 ) .FL(9.4)  C RETURN  192  Listing 569 570 57 1 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 61 1 612 613 614 615 616 617 6 18 619 620 62 1 622 623 624 625 626  r KALMAN4 at C C C  C C C  10:07:47 on DEC 13. 1987 for CCId'SEIS  END This subroutine e x t r a p o l a t e s  the e r r o r  covariance matrix  SUBROUTINE EXTRAP(PU.PS.FL) REAL PS(9.9).PU(9.9).FL(9.9).FLT(9.9) REAL TMPI(9.9).TMP2(9.9) INTEGER OR OR = 9 CALL GTRAN(FL.FLT,OR,OR.OR.OR) CALL GMULT(FL.PS.TMP1.OR.OR.OR.OR.OR.OR)  C CALL GMULT(TMP 1 . FLT . PU, OR . OR . OR, OR .OR. OR) C C C C  C C C  C C C C C  C  RETURN END This subroutine updates the e r r o r c o v i r l a n c e matrix P. SUBROUTINE PUPDAT(HV.KS.PU.PS) REAL KS(9),PU(9,9).HV(9),ID(9.9).TMP1(9.9) REAL TMP2(9.9),PS(9.9) INTEGER OR OR = 9 CALL UNIT(ID.OR.OR) DO 1 1=1.9 DO 2 d=1.9 TMP1(1.J)=KS(I)*HV(J) 2 CONTINUE 1 CONTINUE CALL GSUBUD.TMP1.TMP2,OR.OR.OR.OR.OR) CALL GMULT(TMP2,PU,PS.OR,OR,OR.OR,OR,OR) RETURN END SUBROUTINE  MlNSN(LX.X.XMIN.INDEX)  REAL X(2) XMIN=X(1) DO 10 I = 1.1.* 10 XMIN=AMIN1(XMIN.X(I)) DO 20 J=1.LX 1NDEX=J IF(X(d)-XMIN) 20.30.20 20 CONTINUE 30 RETURN END  193  L i s t i n g of KALMAN4 at 10:07:47 on DEC 13. 198*7 f o r C C i d - S E I S 627 628 629 630, 631 \ 632 633 634 635 636 637 638  SUBROUTINE MAXSN(LX.X.XMAX.INDEX) C REAL X(2) XMAX-X(l) DO 10 1-1.LX 10 XMAX-AMAXKXMAX.X(I)) DO 20 J - 1 . L X 1NDEX-J IF(X{JJ-XMAX) 2 0 . 3 0 . 2 0 • 20 CONTINUE 30 RETURN END  194  

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}]}"
                            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:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0062674/manifest

Comment

Related Items