- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Application of digital filtering techniques for reducing...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Application of digital filtering techniques for reducing and analyzing in-situ seismic time series 1988
pdf
Page Metadata
Item Metadata
Title | Application of digital filtering techniques for reducing and analyzing in-situ seismic time series |
Creator |
Baziw, Erick John |
Publisher | University of British Columbia |
Date Created | 2010-09-09T21:12:41Z |
Date Issued | 2010-09-09T21:12:41Z |
Date | 1988 |
Description | The introduction of digital filtering is a new and exciting approach in analyzing in-situ seismic data. Digital filters are also in the same spirit as the electric cone which replaced the mechanical cone in CPT* testing. That is, it is desirable to automate CPT testing in order to make it less operator dependent and increase the reliability and accuracy. In CPT seismic cone testing seismic waves are generated at the surface and recorded downhole with velocity or acceleration transducers. The seismic receivers record the different seismic wavelets (e.g., SV-waves, P-waves) allowing one to determine shear and compression wave velocities. In order to distinguish the different seismic events, an instrument with fast response time is desired (i.e., high natural frequency and low damping). This type of instrument is characteristic of an accelerometer. The fast response time (small time constant) of an accelerometer results in a very sensitive instrument with corresponding noisy time domain characteristics. One way to separate events is to characterize the signal frequencies and remove unwanted frequencies. Digital filtering is ideal for this application. The techniques of digital filtering introduced in this research are based on frequency domain filtering, where Fast Fourier, Butterworth Filter, and crosscorrelation algorithms are implemented. One based on time domain techniques, where a Kalman Filter is designed to model'the instrument and the physical environment. The crosscorrelation method allows one to focus on a specific wavelet and use all the information of the wavelets present averaging out any noises or irregularities and relying upon dominant responses. The Kalman Filter was applied in a manner in which it modelled the sensors used and the physical environment of the body waves and noise generation. The KF was investigated for its possible application to obtaining accurate estimates on the P-wave and S-wave amplitudes and arrival times. The KF is a very flexible tool which allows one to model the problem considered accurately. In addition, the KF works in the time domain which removes many of the limitations of the frequency domain techniques. The crosscorrelation filter concepts are applied by a program referred to as CROSSCOR. CROSSCOR is a graphics interactive program which displays the frequency spectrums, unfiltered and filtered time series and crosscorrelations on a mainframe graphics terminal which has been adapted to run on the IBM P.C. CROSSCOR was tested for performance by analyzing synthetic and real data. The results from the analysis on both synthetic and real data indicate that CROSSCOR is an accurate and user friendly tool which greatly assists one in obtaining seismic velocities. The performance of the Kalman Filter was analyzed by generating a source wavelet and passing it through the second order instrumentation. The second order response is then fed into the KF with the arrival time and maximum amplitude being determined. The filter was found to perform well and it has much promise in respect that if it is finely turned, it would be possible to obtain arrival times and amplitudes on line resulting in velocities and damping characteristics, respectively. * Cone Penetration Test |
Subject |
Electric Filters, Digital Kalman Filtering |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | Eng |
Collection |
Retrospective Theses and Dissertations, 1919-2007 |
Series | UBC Retrospective Theses Digitization Project [http://www.library.ubc.ca/archives/retro_theses/] |
Date Available | 2010-09-09T21:12:41Z |
DOI | 10.14288/1.0062674 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
URI | http://hdl.handle.net/2429/28364 |
Aggregated Source Repository | DSpace |
Digital Resource Original Record | https://open.library.ubc.ca/collections/831/items/1.0062674/source |
Download
- Media
- UBC_1988_A7 B39.pdf [ 8.94MB ]
- Metadata
- JSON: 1.0062674.json
- JSON-LD: 1.0062674+ld.json
- RDF/XML (Pretty): 1.0062674.xml
- RDF/JSON: 1.0062674+rdf.json
- Turtle: 1.0062674+rdf-turtle.txt
- N-Triples: 1.0062674+rdf-ntriples.txt
- Citation
- 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 Universi ty 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 th is thesis 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 this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. The University of British Columbia 1956 Main Mall Vancouver, Canada Department V6T 1Y3 nF.firt/ft-n 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 exc i t ing approach in analyzing i n - s i t u seismic data. D ig i ta 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 test ing 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 test ing seismic waves are generated at the surface and recorded downhole with ve loc i ty or accelerat ion transducers. The seismic receivers record the d i f fe rent 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 is t inguish the d i f fe rent seismic events, an instrument with fast response time is desired ( i . e . , high natural frequency and low damping). This type of instrument is charac te r i s t i c of an accelerometer. The fast response time (small time constant) of an accelerometer resul ts in a very sens i t i ve 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 ig i ta l f i l t e r i n g is ideal for th i s app l i ca t ion . The techniques of d i g i t a l f i l t e r i n g introduced in th is research are based on frequency domain f i l t e r i n g , where Fast Four ier , Butterworth F i l t e r , and crosscorre lat ion 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 crosscorre lat ion 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 re ly ing 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 appl icat ion to obtaining accurate estimates on the P-wave and S-wave amplitudes and a r r i v a l times. The KF is a very f l e x i b l e tool which allows one to model the problem considered accurately . In add i t ion , the KF works in the time domain which removes many of the l imi ta t ions of the frequency domain techniques. * Cone Penetration Test i i The crosscorre lat ion f i l t e r concepts are applied by a program referred to as CROSSCOR. CROSSCOR i s a graphics in teract i ve program which displays the frequency spectrums, unf i l te red and f i l t e r e d time ser ies and crosscorre lat ions on a mainframe graphics terminal which has been adapted to run on the IBM P.C. CROSSCOR was tested for performance by analyzing synthet ic and real data. The resul ts from the analysis on both synthetic and real data indicate that CROSSCOR is an accurate and user f r iend ly tool which greatly ass is ts 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 into the KF with the a r r i v a l time and maximum ampli - tude 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 ine resu l t ing in v e l o c i t i e s and damping charac- t e r i s t i c s , respect ive ly . i i i TABLE OF CONTENTS Page ACKNOWLEDGEMENTS xv ABSTRACT 1i LIST OF TABLES vi LIST OF FIGURES v i i i I. INTRODUCTION 1 I I . DESCRIPTION OF PHYSICAL PROBLEM 8 I I I . APPLICATION OF FREQUENCY DOMAIN TECHNIQUES 22 A. Discussion of Crosscorrelation F i l t e r Formulation 22 B. F ie ld Program and Discussion of Results 48 B . l Research Sites 51 B.2 Cross-Over Method Vs. CROSSCOR with Analogue F i l t e r Applied and Different Seismic Probes Used 61 B.3 CROSSCOR with F i l te red (Analogue) and Unf i l tered ( i . e . , Removing P-Waves and Noise with CROSSCOR) Sei smi c Traces . 72 B.4 CROSSCOR with Nonpolarized Sources 80 B.5 CROSSCOR in Obtaining Compression Wave Ve loc i t ies 89 B.6 Bandwith Considerations when Applying CROSSCOR... 98 B.7 Conclusions 107 IV. KALMAN FILTER FORMULATION FOR ESTIMATING AMPLITUDES AND ARRIVAL TIMES 109 A. Kalman F i l t e r Results I l l B. Description of Continuous Version of the Kalman F i l t e r 123 C. Description of Discrete Version of the Kalman F i l t e r 127 D. Derivation of Transfer Function for Geophones 131 iv TABLE OF CONTENTS (Cont'd) Page E. KF Equations for Estimating Amplitudes and Ar r i va l Times 136 F. Comments on Computing Estimation Error Covariance Matr ix , P 147 G. 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 ; 165 B. Discrete Model of Second Order System 178 C. Program L i s t i n g of the KF Fomulation 181 v LIST OF TABLES No. T i t l e Page B.2.1 Calculated Veloci ty 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 B.2.2 Calculated Veloci ty 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 iezoe lec t r i c Bender Accelerometers 65 B.2.3 Calculated Veloci ty 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 Appl ied. In th i s Set of Data i t was not Possible to Obtain Accurate Cross- Over Picks 68 B.2.4 Calculated Veloci ty 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 B.3.1 Calculated Veloci ty 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 Unf i l tered Seismic Data. The Data was Acquired with an Accelerometer 75 B.3.2 Calculated Veloci ty P r o f i l e s from Ti lbury S i te Comparing Shear Wave Ve loc i t ies Obtained with CROSSCOR from Seismic Data Acquired by Two Difference Accelerometers. In th i s Set of Data i t was not Possible to Obtain Accurate Cross-Over Picks 78 B.4.1 Calculated Veloci ty P r o f i l e s of Lower Langley (232) Comparing Shear Wave Ve loc i t ies Obtained from the Hammer Shear Source and the Buffalo Gun Source. The Data was Acquired with an Accelerometer (Unf i l tered) 83 B.4.2 Calculated Veloci ty P r o f i l e s for Annacis P i l e S i t e , where Ve loc i t ies 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 (Unf i l tered) : 85 vi LIST OF TABLES (Cont'd) No. T i t l e Page B.4.3 Calculated Veloci ty P r o f i l e s for Annacis P i l e S i t e , where Ve loc i t ies 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 (Unf i l tered) 87 B.5.1 Congressional and Shear Ve loc i t ies in Sediments and Rocks 9 2 B.5.2 Calculated Veloc i ty P r o f i l e from Lower Langley, where an Accelerometer was Used for Data Acqu is i t i on . The Compression Wave Ve loc i t i es were Obtained from the Buffalo Gun Source and by Applying a Bandpass of 200 to 1000 Hz 9 3 B.5.3 Calculated Veloci ty P r o f i l e from Lower Langley where an Accelerometer was Used for Data Acqu is i t i on . The Compression Wave Ve loc i t i es were Obtained from the Buffalo Gun Source and by Applying a Bandpass of 750 to 1500 H z . . . . 96 vi i LIST OF FIGURES No. T i t l e Page 1.1 Schematic of an Acceleration Measuring Device 3 1.2 Unit Step Response Curve Showing Transient Response Spec i f icat ions 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 Veloci ty Calculat ion 7 2.1 Schematic of Seismic Transmission P r o f i l e 8 2.2 CPT Seismic Equipment Layout for 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 in Wave Type O c c u r r i n g . . . . 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 F i l t e r s 16 2.7 Ef fect of Di f ferent F i l t e r Cut -off Frequencies on Accelerometer Responses 19 2.8 Seismic Data Acquired from Annacis Island Vibro- Compaction S i te on May 18, 1988. Data Recorded with Accelerometer having a 300 Hz Low-Pass F i l t e r Appl ied. Diagram I l l u s t r a t e s D i f f i c u l t y in Obtain- ing Cross-Overs Due to Signal being Masked by Many Low Frequencies 20 2.9 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 Appl ied. The Diagram I l l u s t r a t e s Ideal Signals for Determining Cross-Overs 21 3.1 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 vi i i LIST OF FIGURES (Cont'd) No. T i t l e Page 3.2 Sampling and A la i s ing Dif ferent 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 ie lds 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 inc ( .5at ) 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 I l l u s t r a t i o n of the A la i s ing Problem when F i l t e r i n g 31 3.10 Relationship Between the Deformed Angular Frequencies wd, for 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 32 3.11 An I n f i n i t e l y Long Sinusoidal and i t s ' . Fourier Transform 34 3.12 A Truncated 60 Cycle Signal and i t s ' Fourier Transform. The Signal was Switched on for #/2 Seconds 35 3.13 A Data Window Used to Taper a Discrete Function 36 3.14 Autocorrelat ion of a Process 37 3.15 Calculat ing the Crosscorrelat ion of Two Functi ons 37 3.17 Padding Signals with Zeros at the End of the i r Time Series 38 3.18 Autocorrelat ion of Box Car Function 39 i x LIST OF FIGURES (Cont'd) No. T i t l e Page 3.19 Signal with D.C. Offset of -0 .20 40 3.20 Signal with D.C. Offset of -0 .06 40 3.21 Crosscorrelat ion 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 Crosscorrelat ion Value for D.C. S h i f t Corrected Signals 42 3.25 Superimposed Signals to give Visual Aid for Estimating Time Lag 43 3.26 Processing Synthetic Data with CROSSCOR. (a) Input Trace (b) Frequency Spectrums (c) Bandpassed Traces (40 to 60 Hz) (d) Crosscorrelat ion Function of the Fi l tered Traces 47 B . l General S i te Location Map : 49 B . l . l Cone Penetration P r o f i l e for McDonald Farm S i te 52 B.1.2 Cone Penetration P r o f i l e for Grant McConachie Way S i te 54 B . l . 3 Cone Penetration P r o f i l e for T i lbury Natural Gas Plant S i te 56 B.1.4 Cone Penetration P r o f i l e for Annacis P i l e Research S i te 58 B . l . 5 Cone Penetration P r o f i l e for Lower Langley 232 St . S i te 60 B.2.1 Seismic Data Acquired from Annacis Island Vibro- Compaction S i te on May 18, 1988. Data Recorded with Accelerometer Having a 300 Hz Low-Pass F i l t e r Appl ied. Diagram I l l u s t r a t e s D i f f i c u l t y in Obtain- ing Cross-Overs due to Signal being Masked by Many Dominant Low Frequencies. 52 x LIST OF FIGURES (Cont'd) No. T i t l e Page B.2.2 (a) Veloc i ty P r o f i l e from McDonald Farm Comparing the Cross-Over Method and CROSSCOR. The Data was Obtained with a Se 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 B.2.3 (a) Veloc i ty 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 Appl ied, (b) Cone Bearing P r o f i l e from McDonald Farm 67 B.2.4 (a) Veloci ty 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 Appl ied. In th is Set of Data i t was not Possible to Obtain Accurate Cross- Over P icks , (b) Cone Bearing P r o f i l e from the Laing Bridge S i t e . . . 69 B.2.5 (a) Veloci ty 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 Appl ied, (b) Cone Bearing P r o f i l e from Lower Langley S i te 71 B.3.1 (a) Seismic Trace Recorded at Grant McConachie. Way with the P-Plate Source (Unf i l te red) . (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 B.3.2 (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) Correspond- ing Crosscorre lat i on Function 74 B.3.3 (a) Veloc i ty P r o f i l e from Lower Langley Comparing F i l t e r e d (300 Hz Low-Pass) and Unf i l te red Seismic Data. The Data was Acquired with an Accelerometer. (b) Cone Bearing P r o f i l e from Lower Langley 76 XT LIST OF FIGURES (Cont'd) No. T i t l e Page B.3.4 (a) Calculated Veloc i ty P r o f i l e s from Ti lbury S i te Comparing Shear Wave Ve loc i t i es Obtained with CROSSCOR from Seismic Data Acquired by Two Dif ferent Accelerometers. The Veloc i ty 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 Pro- f i l e . In th is Seismic Data i t was not Possible to Obtain Accurate Cross-Over P icks , (b) Cone Bearing P r o f i l e from Ti lbury Island S i te 79 B.4.1 Seismic Section from Lower Langley where a Buffalo Gun was Used as the Source and an Unf i l tered Accelerometer Recorded the Data. From the Figure one can Clearly Distinguished the Compression and Shear Waves 81 B.4.2 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 B.4.3 Calculated Veloc i ty P r o f i l e s from the Annacis P i l e Research S i t e , where Ve loc i t ies Obtained from Both Sides of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer (Unf i l te red) . (b) Cone Bearing P r o f i l e of the Annacis Bridge P i l e Research S i t e . . . 86 B.5.1 (a) Seismic Trace from Grant McConachie Way where a P-Plate Source was Used and an Accelerometer (Unf i l tered) 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 B.5.2 Veloc i ty 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 B.5.3 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) in Order to Extract the Compression Wave Responses. The Resulting Compression Wave Responses are Highly Correlated 95 xi i LIST OF FIGURES (Cont'd) No. T i t l e Page 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 B.6.1 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 Unf i l tered Accelerometer (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response at 50 Hz 99 B.6.2 Seismic Trace Recorded at Lower Langley where a Hammer Shear Source was Used, (a) Seismic Trace at 3.7 Metres Acquired with an Unf i l tered Accelerometer (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response at 50 Hz 100 B.6.3 (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 Veloc i ty of 103 m/sec. (c) Traces in (a) F i l t e r e d Between 30 to 70 Hz Result ing in a Veloci ty of 103 m/sec. 101 B.6.4 Seismic Trace Recorded at Lower Langley where a Buffalo Gun Source was Used, (a) Seismic Trace at 3.7 Metres Acquired with an Unf i l tered 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 B.6.5 Seismic Trace Recorded at Lower Langley where a Buffalo Gun Source was Used, (a) Seismic Trace at 4.7 Metres Acquired with an Unf i l tered Accelero- meter (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 B.6.6 (a) Seismic Traces Recorded at Depths of 3.7 and 4.7 Metres, (b) Trace in (a) F i l te red between 20 to 200 Hz Resulting in a Veloc i ty of 94 m/sec. (c) Traces in (a) F i l te red between 60 to 180 Hz Resulting in a Veloci ty of 94 m/sec 105 B.6.7 Trace in Figure B.6.6 (a) F i l t e r e d between 100 to 200 Hz Resulting in a Veloci ty of 94 m/sec 106 xi i i LIST OF FIGURES (Cont'd) No. Title Page 4.1 Block Diagram of System, Measurements, and Kalman Filter 110 4.2 Schematic Illustrating the Interaction of Input Noise, Body Wave and Recording Instrument n ? With Resulting Output 4.3 The Effect An Increasing Time Constant Has On White Noise 113 4.4 Low-Pass Noise With Corresponding Geophone Response as the Time Constant Increases 114 4.5 Geophone Response Due to P and S-Wave Impulses and Noise With Increasing Variance 117 4.6 Illustrating the Synthetic (a) Input and (b) Geophone Response 119 4.7 Updating the States (a) x(2) (Arrival Time) and (b) x(l) (Amplitude) 121 4.8 Mean Square Error for (a) State x(2) and (b) State x(l) 122 4.9.a Block Diagram of State & Measurement Equations (4.1) & (4.2) 126 4.9.b 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 132 4.12 Block Diagram of Second Order Geophone Transfer Function 136 xiv ACKNOWLEDGEMENTS The work outlined i n t h i s thesis would not have been possible without the funding and guidance of Dr. Carapanella. Dr. Ken W h i t t a l l from the Geophysics Department was fundamental i n helping me with the geophysical aspects of th i s t h e s i s . In addition, Drs. Robertson, Dunbar and Slauson provided considerable advice i n preparing t h i s t h e s i s . The rest of the i n - s i t u group's help was invaluable. Don G i l l e s p i e and John Sully were always available for assistance i n formulating my ideas. I would also l i k e to thank Jim Greig and Michael Ehling for t h e i r help with the in t e r a c t i v e graphics. I would l i k e to express my thanks to my buddies at the Savage Mansion, Damon, Steve and Bear for keeping my s p i r i t s up at a l l times. Las t l y , my parents were invaluable i n helping me prepare my research and typing i t up (Thanks, Mom!). xv I. INTRODUCTION The objective of th i s research is 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 is considerable interest in these s o i l charac te r i s t i cs because they provide ins ight into the response of s o i l to imposed loads such as bu i ld ings , heavy equipment, and dynamic loads from earthquakes and explosions. Richart et a l . (1970), Mooeny (1974) and Borm (1977) have wri t ten extensive l i t e r a t u r e on the re la t ion 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 test was f i r s t investigated at UBC by Rice (1984) where he found the appl icat ion of the technique could provide a rapid and accurate method for carrying out a downhole shear wave ve loc i ty survey. In his work, Rice determined that a uniaxial geophone was s u f f i c i e n t for ident i f y ing shear waves, i f the instrument was orientated properly. Rice also investigated several d i f fe rent 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 ampl i f ied at the p ick -up , Rice stated "Below 30 metres the strength of indiv idual energy impulses was attenuated such, that c lear 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 invest igate the appl icat ion of the seismic cone in i n - s i t u test ing was Laing (1985). Laing focused her attent ion on sources and receivers with the seismic cone for both on-shore and of f -shore app l icat ions . Laing came up with the same conclusion as Rice with respect to the use of the geophone as a seismic rece iver , she s ta tes , "accelerometers tend to show a response which more c lose ly represents the response of the s o i l . Thus, damping character i s t i cs of s o i l should theore t i ca l l y be attainable from accelerometer responses." Ideal ly one would l i k e an instrument to t e l l the whole seismic story . Figure 1.1 i l l u s t r a t e s a typ ica l accelerat ion and/or ve loc i ty measuring device. The dynamics behind the Veloci ty and accelerat ion measurements can be expressed by the fol lowing two second order equations: 1 acce lerat ion : X 0 ( S ) / s 2 X i = X 0 (s )/X i = k/((s2/a)^)+2^s/u)n+l) (1.1) v e l o c i t y : x Q ( s )/sx i = x Q (s )/x . = s/(s 2+2^ ns+a) 2) (1.2) x 0 r e l a t i v e displacement of mass M x. = absolute displacement of mass M u n =.natural frequency X = damping k = proport ional i ty constant s = Laplace transform var iable The response of the instrument is described by w n and c In general , the higher the value of x, and smaller value of w n the greater the reduction in the s e n s i t i v i t y of the instrument. This corresponds to an instrument with a large time constant resu l t ing in slower transient response s p e c i f i c a t i o n s , th is is i l l u s t r a t e d in Figure 1.2. In th is 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 typ ica l frequency response curves for geophones and accelerometers. As is 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 for lower frequency inputs , and becomes unstable for 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 bet- ween 50 to 300 Hz, i t is desirable to have an instrument which has a f l a t res - ponse curve below 2000 Hz and is 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 accelero- meter with natural frequency at approximately 4000 Hz. The desire of using sens i t i ve accelerometers in i n - s i t u test ing c a l l s for 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 po lar i t y of shear body waves and the appl icat ion of analogue f i l - t e rs . In Section I I , th is technique is out l ined 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 B K D N d l R a c(t) - mass of suspended magnet - flux density between magnets - stiffness of spring which supports the magnet - viscous damping between magnet and the case - number of turns in coil - diameter of coil - length of co i l , ir-d-N - 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 3 TRANSIENT RESPONSE SPECIFICATIONS 1. Delay time td 2. Rise time tr 3. Peak time tp 4. Maximum overshoot Mp 5. Settling time t, 4 0) c o o a> u o •a: d •rt C e Q 8 10 15 2? 25-3040 5060 80 100 150200 300 403fj 600 1000 2000 3000 4000 025 6 8 100 200 300 400 500 1000 Frequency, Hz 2000 3000 4000 Figure 1.3. Accelerometer and Geophone Frequency Response Curves. The Geophone I l lus t ra ted has a Natural Frequency of 28 Hz, and the Accelerometer has a Natural Frequency.of 500 Hz. The dynamic behind source wavelet generation and f i n a l ve loc i ty deter- mination i s as fo l lows . A source, such as a hammer blow or seismic cap, generates potent ia ls at the surface which correspond to SV*, SH**, and P*** waves. These i n i t i a l potent ials can be considered to be Dirac del ta functions ( i . e . , impulses) having a broad frequency spectrum. As these potent ia ls propa- gate through the earth as seismic waves, the i r frequency spectra are bandlimited ( i . e . , higher frequencies are attenuated). This i s due to the earth act ing as a low-pass analogue f i l t e r . At a cer ta in depth, the generated s ignals 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 superposit ion of many s igna ls . 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 is app l ied , the time of fset between two s i m i l a r waves ( e . g . , SH) recorded at successive depths is determined and interval v e l o c i t i e s are ca lcu lated . The fol lowing diagram (Figure 1.4) i l l u s - t rates the essent ial r e l a t i o n between source, ear th , recorder, f i l t e r , and interval ve loc i ty c a l c u l a t i o n . Referring to Figure 1.4, the major goal was to develop d i g i t a l f i l t e r s . These f i l t e r s should give c lear s ignals without d i s to r t ing the o r ig ina 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 th is thesis are the one based on the frequency domain f i l t e r i n g (Section I I I ) , where Fast Four ier , Butterworth, and crosscorre lat ion 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 hor i zonta l l y t rave l ing 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 hor izonta l ly t ravel ing shear wave so polar ized that the p a r t i c l e motion is a l l hor i zonta l . *** Compression wave consist ing of p a r t i c l e motion with a l ternat ing condensations and rarefact ions during which adjacent pa r t i c les of the s o l i d are c loser together and farther apart during successive half cyc les . The motion of the par t i c les are always in the d i rec t ion of wave propagation. 6 Depth 1 I Noise P-Wave S-Wave I 1 Source P o t e n t i a l s Earth Model Recording Equipment Band-Pass F i l t e r Depth 2 SV SH 5^1 I Noise P-Wave S-Wave I < 5 > V AV = 6d Superimposed Signal V e l o c i t y C a l c u l a t i o n s Figure 1.4. Block Diagram of S o u r j ^ . J i a r t h , Recorder, F i l t e r and Veloci ty Ca lcu lat ion . 7 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 ve loc i t ies of a strat igraphic sect ion. The shear wave ve loc i t ies are d i rect ly related to the maximum shear modulus (Rice 1984) of the so i l by the following mathematical expression. 2 max i 's G_... = f • V where ^ - i s the mass d e n s i t y o f the s o i l V - is the shear wave ve loc i ty s J G - i s the maximum shear modulus a t s t r a i n s l e s s t h a n max 10 percent. The signif icance of the above strain level is that at strains of 10" per- cent and less , the shear modulus is maximum and l i t t l e effected by s t ra in amplitude. The maximum shear parameter is an important so i l parameter; shear modulus reduction curves can be generated from i t . The shear modulus reduction curves are very important in dynamic so i l analysis,because they relate the change in the shear modulus with induced s t ra ins . In seismic cone test ing a response signal is recorded at one metre interval depths. The wavelet is 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 . recording depth to < recording sledge hammer bl ow cone inserted here d' - d" - ! j < * i v s i * T l " 1 i 1 , d 2 Vs2 s l 1 ! ' d ' d " • — i d 3 _ i ... »S3 T —rr - 1 " * T 2 - arrive to f i r s t seismic recei ver time for shear to arr ive 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 Nicolet 4094 d i g i t a l osc i l loscope with a CRT screen and a floppy disk storage c a p a b i l i t y . A complete descr ipt ion of th i s osc i l loscope is given by Rice (1984). The seismic cone is pushed into the s o i l using the U.B.C. In - s i tu Test- ing vehic le . 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 cy l inders with which the cone is pushed into the s o i l . To provide the reaction force needed to push the cone, the truck is raised on two hydraulic pads as shown in Figure 2 .2 . The hydraulic pads also act as SH- Wave sources. SH-Wave sources are desired because i t i s easy to reverse the po la r i t y of the SH-Wave. This energy reversal resul ts in two oppositely polar ized wavelets which are recorded on the Nicolet osc i l loscope as is 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 ca lcu la t ion of the r e l a t i v e travel time dif ference of shear waves between successive depths which allows the ca lcu la t ion of shear wave v e l o c i t i e s . The above technique is thoroughly discussed by Rice (1984) and Laing (1985), and w i l l be referred to as the reverse po la r i t y method and/or crossover method. The crossover times discussed above are related to the ve loc i ty p r o f i l e by the fol lowing equations: d II - d ' V corr corr (2.2a) s3 where d corr ii (2.2b) d corr (2.2c) and d , d " , T, , and T 9 are defined in F igu re .2 .1 . 9 Figure 2.2 CPT Seismic Equipment Layout for Mechanical Source, (after Rice, 1984) 10 one-way travel time (msec) f i r s t a r r i v a l s 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 ra t io of the interval shear v e l o c i t i e s i s of the order of uni ty . These conditions w i l l general ly 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 for 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 is not separated much from the rece iver , 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 neg l ig ib le 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 re f ract ion (Figure 2.4b). The same can be said for the P-Wave. An incident P-Wave subjected to these boundary conditions w i l l transmit most of i t s energy, at an in te r face , 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 ra ight 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 ta tes , "For most p rac t i ca l purposes a s t ra ight l i n e t ravel path assumption would be quite s a t i s f a c t o r y . " As was previously discussed, the interval t imes, T̂ and J^, are determined by taking the di f ference of the crossover times for the successive wavelets. When only noise free shear waves e x i s t , the above procedure is trouble f ree . In many cases, one is not dealing with an ideal wavelet, but has a seismogram with corrupting P-Waves and noise. These par t i cu la r s ignals make i t d i f f i c u l t to pick a r r i v a l and crossover times. The method which is current ly used to obtain c lear 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 - to 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 lect ron ic 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, un f i l te red signals were generally preferred and were used in the i r study. A numerical example of the possible phase s h i f t or corresponding time delay, for an analogue f i l t e r , fo l lows. 12 a7 ' E i n a l 2 3 4 Pi 1.036 1.07 1.0'J 1.10 1.00 1.103 1.07 •7T <*l °2 0.24 0.24 1.09 0.25 0.26 1.14 1.30 0.28 0.16 1.286 1.2tl6 0.25 0.25 0«tai l i doubtful 20° 40° 60° Angle of Incidence 80' p - density of s o i l a - shear wave veloci ty B - compression wave veloc i ty a - attenuation • {*>) S K . ijji • i .1 II ' ]i /I' s v c t / <H • ^ - s \ </ Ii J i L ' — i — Figure 2.4. Ratio of Reflected or Transmitted Energy to Incident Energy, "V E o u t / E - j n > M i t h N o Change in Wave Type Occurring (Ewing, Jardetz , and Press, 1957) Figure 2. 5. Ratio of Reflected or Transmitted Energy to Incident E n e r g y , V E Q U t / E i n , When Waves Change Type (Ewihg, Jardetz, and Press, 1957) The accelerometers used in th is invest igat ion were bimorph transducers, which were 0.5 inches square in dimension and made of G1195 piezoceramic bender mater ia 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 th is accelerometer can be obtained in technical papers from Piezo E l e c t r i c Product's, Inc . , New Jersey. Figure 2.6 i l l u s t r a t e s the phase response curve for a typ ica l f i r s t - o r d e r low-pass analogue f i l t e r s i m i l a r to f i l t e r s within the Nicolet recorder. Equation (2.3a) i l l u s t r a t e s the interact ion of the input wavelet, f i l t e r , and output wavelet. |AU)|e input 1 |Y(«o)|e ' f i l t e r |A(M)||Y( u)|e output j(*1+<j>2) (2.3a) define A(UJ) = |A(w)|e 1 Y(co) = |Y(o.)|e J4>< G(u>) = |A(to) | | Y(ui) | e Input F i l t e r Output (2.3b) (2.3c) (2.3d) angular frequency If i t i s assumed that the f i l t e r has unity ga in , then Equation (2.3d) reduces t o , j(+1+<|)2) G(u>) = |A(u)|e (2.3e) The Fourier transform, for an aperiodic wavelet, of Equation (2.3'e) i s CO i r Jwt g(t) = ~- J G(o))e dto .... (2.3f) In terms of the frequency,•f, we have 15 ,.25 .3 FREQUENCY .6 .7 . 8 . 91 .0 2.0 3.0 4.0 -90u (a) Low-Pass Response K - Time constant S - Laplace transform var iable 0° FREQUENCY 25 .3 .4 .5 .6 .7 . 8 . 9 1.0 3 -45°- < X a. -90°- • out KS • —— - in -> • * 2.0 3.0 4.0 (b) High-Pass Response Figure 2 .6. Phase Response of f i r s t Order Low-Pass and High-Pass F i l t e r s (Lancaster, 1982) 16 r j2irft g(t) = / G(f)e df (2.3g) Subst i tut ing Equation (2.3e) into Equation (2 .3g) , gives / |A(f)|e j(27Tft+Y(f)) g(t) df (2.4a) -co where Y(f ) <l'1(f)+'t'2(f) For an actual waveform, g(t) i s real and hence Equation (2.4a) reduces to If in Equation (2.3a) i t i s assumed that <j>̂=0 ( i . e . , there is no phase lag in the input s i g n a l ) , then We next consider varying the frequency of the input wavelets and determine actual time s h i f t s by re fer r ing to Figure 2 .6 , which is normalized to a cutoff 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 with in these ranges. If our input wavelet has a frequency of 50 hz, then from Figure 2 . 6 , we have <t>?=-25°. From Equation (2 .4c ) , we also have At maximum amplitude, the cosine term in Equation (2.4b) should equal 1 .0 ; (2.4b) - co y ( f ) = <l>2 (2.4c) Y ( f ) = -25° = -0.436 rad 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° = - .960 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 ference 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 for d i f fe rent applied c u t - o f f frequencies. An addit ional consideration is that the reverse po la 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 inves t iga t ion . But i t has been found that in numerous cases th is 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 f i l t e r was appl ied) . Ideal wavelets for determining crossovers are i l l u s t r a t e d in Figure 2 .9 . 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 crosscorre lat ion 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 ign i f i cance on dominant responses), th is method should provide better and more r e l i a b l e resul ts than the reverse po la r i t y method. Detai ls of th is 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 shear wave a r r i v a l s " •I C u t - o f f . frequency = 20 Hz (dela y due Co phase s h i f t ) j 1 1 ^ % I .IWVIAAA. 600 Hz 500 H: 1000 Hz 20 40 GO 30 Time (msec) 100 120 140 igure 2.7. E f f e c t of D i f f e r e n t F i l t e r Cut-off Frequencies on Accelerometer Responses (Laing, 1985) 19 I 1 1 1 1 1 1 I L_ 30.0 60.0 90.0 Time (msec) Figure 2.8. 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. Dia- gram 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 . 44.61 Oppositely Polarized SH-Waves at 2.1m 48.89 53.77 Oppositely Polarized SH-Waves at 3.1m Oppositely Polarized SH-Waves at 4.1m 20 Figure 2.9. 40 60 100 110 120 140 80 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 lus t rates Ideal Signals for Determining Crossovers I I I . APPLICATION OF FREQUENCY DOMAIN FILTERING A. Discussion of the Crosscorrelat ion 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 lags , and time s h i f t s , one should have a fundamental idea of what is 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 is corrupted with many d i f fe rent components. Most f i l t e r s permit the select ion of the upper and the lower l i m i t s of the passband. This concept is i l l u s t r a t e d in the fol lowing diagram: A(w) Y D X(u)) OJ +-> «=c f f f f T o T l r 2 T 3 Passband f o f l Frequency Figure 3.1 . 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 usual ly possible to select the sharpness of the cutoff (the rate 9 at which the ga in , |Y(")| , decreases as i t leaves the passband). Typical f i l t e r response curves are speci f ied by the 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 is c r i t i c a l l y dependent on the sampling rate of data a c q u i s i t i o n . The sampling theorem defines the minimum sampling interva l one can apply data acqu is i t ion without los ing a wavelet's ident i t y ( a l i a s i n g ) . In general , no information is l os 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 is i l l u s t r a t e d in Figure 3 . 2 , where we have a 4 ms sampling interval for 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 faster than the s i g n a l . 22 (c) Figure 3.2. Sampling and A la is ing Dif ferent 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. Geldart , 1983) Half the sampling frequency is 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 indist inguishable from the lower frequency fN-Af . In summary, one can state that wavelet g(t) can be recovered exact ly from the sampled values provided g(t) does not contain frequencies higher than the Nyquist frequency. When taking the Fourier transform of g ( t ) , G ( f ) , i t w i l l only contain information about g(t) for 23 frequencies within the interval - fN to fN. The Fourier ser ies w i l l repeat G(f) in each interval 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. F i l t e r s pass the desired information from a source wavelet. Usually f i l t e r s attempt to remove spec i f ied frequency components of a sampled wavelet, which sometimes resul ts in the f i l t e r a l t e r - ing the amplitude and/or phase spectra of s ignals which pass through i t . This a l te r ing resul ts because a f i l t e r i s defined by i t s own amplitude and phase spectra just 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 interact ion of the wavelet and f i l t e r . This output i s then defined by i t s own amplitude and phase spectra. This process is 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(w)|e Input Y(ai)=|YU) |e j ^ 2 - X(cu) '• Output where X(o>) |AH||Y(<o)|e' Output Phase Spectra i= l ,2 A( w ) 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 Frequency f i l t e r s are general ly c l a s s i f i e d as low-pass, high-pass band-pass, respect ive ly , as they discr iminate against frequencies above or below a certa in 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 is one which has a gain of unity (|g("0| = 1) f ° r frequencies within the passband and zero for those within the stopbands. Moreover, i t i s advantageous that the f i l t e r produce zero phase s h i f t , so that the phase spectra, <f>2(f), i s zero for a l l frequencies. " Ideal" f i l t e r s of these types are the fol lowing Low-pass G L ( W ) = +1, |u)|<|a>0| 0, | u |<| U o | High-pass G H U ) = 0, |co|<|a> | + 1, |u|>|a>0|- Band-pass G g U ) = +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 w0. o l 2 The low-pass f i l t e r i s obtained by applying a box car function with a l i n passband oi . The i t s associated transform. certa box car function is i l l u s t r a t e d in Figure 3 .4 , with _-«/:» \ sine { -8»/o - 6 tia 2* la O Hn/a I where s inc( .5at ) = 2sir i ( ;5at) at Figure 3.4. Example of a Box Car Function, G(ai), and Its Transform, Sine ( .5at) (R. E. Sher i f f & L. P. Geldart , 1983) 25 The low-pass f i l t e r i s defined by the fol lowing expression: G L(u) = box2w0(w) (u0/Tr)sinc(ai0t) = g L ( t ) (3.1) where «- -»• denotes the Fourier transform. For d i g i t a l funct ions, provided |u o | < = ir/At, = Angular Nyquist frequency L §\ = (1/ir) £ wo s inc(nw 0At) (3.2) n=-L where L -*• « It can mathematically be proven that the expression of the high-pass f i l t e r is essent ia l l y the same as that for the low-pass f i l t e r , where g H ( t ) = - g L ( t ) , pro- vided both f i l t e r s have the same cutoff frequency <d . In deriv ing the passband f i l t e r one just manipulates the high and low-pass f i l t e r s to interact where the passband l i e s within the desired l i m i t s . The resul t of using a f i n i t e series for gj; is to introduce r i p p l e s , both within and without the passband, the ef fect i s espec ia l l y noticeable near the cutoff frequency. This i s referred to as the Gibbs' Phenomenon. Gibbs' Phenomenon describes the d is tor t ions which occur at d iscont inu i t ies within f i l t e r e d s igna ls . This phenomenon can be conceptualized in the time domain. A signal truncated due to a short recording t ime, i s analogous to mul t ip ly ing the signal with a box car function which results in d is tor t ions or " r ing ing" in the frequency domain. As the s i g n a l ' s recording time increases " r ing ing" becomes less prevalent because more information of the signal i s transformed and the f i n i t e series expressed in Equation (3.2) i s carr ied out fur ther . These concepts are i l l u s t r a t e d in Figure 3.5(b) where a 100 hz cosine signal i s being multipled with a box car funct ion. Figure 3.5(a) i l l u s t r a t e s the "r inging" ef fect in 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 wave- forms. In these f igures a rectangular wave i s constructed by adding a suc- cessive number of harmonics to the fundamental frequency. As is 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 funct ion. 26 (a) 100 ZOO 300 FREQUENCY, Hertz 400 500 (b) 1.0 0.5 to g 0.0 -0.5 1.0 /SO 100 , 150 200 TIME (N = 200,DELT=.00l) NOTE:FREQUENCY = 100 HZ. Figure 3 .5 . 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 0. W. Stimson, 1983). 28 Having discussed the problems associated with designing f i l t e r s , one should now explore ex is t ing f i l t e r s which lessen or remove the sever i ty of undesired f i l t e r i n g e f fec ts . A f i l t e r which has many advantageous q u a l i t i e s i s the Butterworth f i l t e r . Butterworth F i l t e r 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,/ U o) 2 N } (3.3) where u> i s the "cutoff" frequency and N determines the sharpness of the cutof f , o p Curves of |G(w)| for various values of N are shown in Figure 3 .7 . 0 1 2 3 W.'oJ0 normalized frequency Figure 3 .7. Amplitude Response of a Butterworth F i l t e r (R. E. Sher i f f & L. P Geldart , 1983) 29 The advantages associated with the Butterworth f i l t e r are as fo l lows : o Their t ransfer 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 is down by 3 dB (factor of 1/2) at the cutoff frequency. (The cutoff frequency determines the half-power point of the f i l t e r ) . The value of N spec i f ies the rate of attenuation where a larger value of N gives a greater rate of attenuation and " r ing ing" e f fect ( i . e . , Gibbs' Phenomenon). An example of a ca lcu la t ion to determine the value of N necessary i s given below. If the required attenuat ion, measured in units of power, at 2 w o is 48 db down from attenuation at oj o=0, what value of N. is necessary for a Butterworth f i l t e r ? 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 Q (1/(1+W 2 N ) ) W = 0 /(1/(1+W 2 N ) N = 2 = 48 so that (1/(1+2 2 N)) = 1 0 " 4 - 8 Approximately then (1/2 2 N ) = l O " 4 - 8 and N = 4.8/2 log 2 = 4.8/2(.301) = 8 Therefore N=8 is 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 resu l t when the 30 t ransfer function of Equation 3.3 i s appl ied. 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 is noted from Figure 3.8 that the Fourier Transform of th is function is per iodic with period 2TT/W^. A Figure 3.8. Fourier Transform of a Signal with Nyquist to If we t ry to f i l t e r frequencies higher than the Nyquist ( i . e . , apply a cut- off frequency higher than the Nyquist) , higher frequencies w i l l fo ld 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 . CD -a 3 43 1.0 - - Butterworth Window u i -Cuto f f Frequency to (0 N C to Angular Frequency Figure 3.9. I l l u s t r a t i o n of the A l ias ing Problem When F i l t e r i n g ~ ~ 31 In Figure 3 . 9 , frequencies within w* are incorporated into the s i g n a l . In order to avoid th is problem, a b i l i n e a r Z transform i s appl ied . The b i l i n e a r Z transform converts a continuous t ransfer function into one which can be used on sampled data. Furthermore i t el iminates a l i a s i n g errors by applying the fol lowing frequency transformation: w d = | r tan ^ -Tr/At < (o < Ti/At (3.4) where bid = transformed or deformed angular frequency used to ca lcu late the b i l i n e a r Z transform to = angular frequency in the o r ig ina l t ransfer function Equation (3.3) This frequency transformation, wd, causes warping of the input frequencies when they approach the Nyquist frequency. This warping is i l l u s t r a t e d in Figure 3.10. 1 1 tut • 1 / 1 1 i / 1 / ) 1 1 / / , 1 1 1 * I 1 -ws 1 | 1 ' / I 1 1 ' / 1 r / i / i / 1 1 1 deformed frequencies actual frequencies a f te r appl icat ion of b i l i n e a r Z transform Nyquist frequency jd - d i sc re t i zed frequency Figure 3.10. Relationship Between the Deformed Angular Frequencies md, for the B i l i n e a r Z Transform and the Actual Frequencies, m, D ig i ta l F i l t e r (E. R. Kanasevich, 1981). of the 32 The Butterworth f i l t e r , when constructed in the time domain, i s a phys i - c a l l y r e a l i z a b l e , recursive f i l t e r . Phys ica l l y rea l i zab le (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. The Butterworth f i l t e r having the t ransfer function given by Equation 3.3 has the fol lowing time-domain re lat ionsh ip N N Y = y " c X + B Y (3.5) n iTo n n " m t l m n " m where the index N determines the order of the Butterworth f i l t e r . The Butterworth t ransfer 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) Gl ( f ) Butterworth f i l t e r s . For instance, an eighth order Butterworth f i l t e r would be defined by the fol lowing 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 error r e l a t i v e to a brute force ca lcu la t ion of G(f) for each order (Reference 14). Thomson and Chow (1980) discuss how the Z.-domain t rans - fer 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 for the f i r s t order (N=l) and second order (N=2) Butterworth f i l t e r s are, respect ive ly . Y n = d o X n + d l X n - l + e l Y n - l < 3' 7) Y n = c o X n + c l X n - l + c l X n - 2 + b l Y n - l + b 2 Y n - 2 where X n = input Y n = f i l t e r e d output 33 The Butterworth f i l t e r has zero phase lag because a waveform i s front fed and then back fed into the f i l t e r . C lea 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 exact ly . Another consideration to take into account when processing the data, i s the . f i n i t e length of the time ser ies . This f i n i t e length can be viewed as a t runcat ion, . where we mult ip ly an i n f i n i t e length ser ies by a box car funct ion . This m u l t i p l i c a t i o n resul ts in Gibbs1 Phenomenon in the frequency domain, as was previously discussed. These concepts are again i l l u s t r a t e d in Figures 3.11 and 3.12. 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 th is signal gives ideal Dirac deltas at -60 and +60 hz in the frequency domain. Figure 3 .11. An I n f i n i t e l y Long Sinusoidal and I t s ' Fourier Transform 34 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 -Jr- time, t 60 frequency, f Figure 3.12. 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 "r inging" e f f e c t , the data is 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 is usual ly 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. w. (1 + cos 1 (1 + cos ir(n+L) Tr(n-L) )/2 )/2 (L+M)<n<-L L i n l L L < n < L+M (3.9.) where M determines the period of the cosine be l l ( i . e . , the r o l l - o f f of the window), and 2L is 80% of the length of the time ser ies . Hamming (1977) recommends that M be about 10% of the ex ist ing data with 80% in the f l a t part of the window. The data can then be padded with zeros. In th is way, there are no d iscont inu i t ies 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 i I AX j a k A i A" -(L+M) -L L discrete time L+M Figure 3.13. A Data Window Used to Taper a Discrete Function Crosscorrelat ion The f i r s t step in determining the interval t imes, once seismic wavelets are adequately f i l t e r e d , i s to apply c rosscor re lat ions . The crosscorre lat ion of the successive wavelets i s defined as where x i s the time s h i f t between the two time s e r i e s . The crosscorre lat ion provides information between the two wavelets. In th is case we are looking for the a r r i v a l time differences between the dominant responses of the recorded wavelets. With the above considerat ion, one can postulate that when the value of $ X y( . T ) in Equation 3.10 is maximum, we have time s h i f t which is representative of the time interva 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 di f ference in the a r r i v a l times of the dominant responses of the wavelets. Therefore t h e i r c rosscorre lat ion i s essent ia l l y an autocorrelat ion 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 in Figure 3.14. (3.10) 36 Process Autocorrelation Figure 3.14. Autocorrelation of a Process An important aspect to consider when applying crosscorrelat ions to seismic wavelets i s D.C. s h i f t s . D.C. s h i f t s occur when the recorded signal i s not centred at zero mean. Equation 3.10 defined the crosscorre lat ion of two data sets X t and where t is the displacement of Y t r e l a t i v e to X^. This equation i s general ly applied to the signals 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 cor re lat ing the two funct ions, X t ={l , -1, -1/2} and Y t = { 1 , 1 / 2 1 / 2 } , shown in Figure 3.15. X! Z3 Q. E -a 4-> CL E (a) Sh i f t of -2 +1/2 ' 0 - I t QJ 4-> -fa- Time +1 + 1/2 QJ -a Time -3 •1/2 ^ E <f>xy(-2)=0+0+l/2 = 1/2 (b) Sh i f t of - 1 i [•1/2 0 Time • - 1 * +1 " +1/2 0 Time -1/2 (c) Sh i f t of 0 QJ 3 CL E > +1 k +1/2 0 V Time - - 1 QJ T3 3 +-> CL E i >+l i L +1/2 0 1 'Time , +1/2 P x y ( - l ) = 0 - l + l / 4 =-3/4 4, (0) = l - l / 2 - l / 4 = +1/4 Figure 3.15. Calculat ing the Crosscorrelat ion of Two Functions 37 At the end of the signals the functions are padded with zeros in order that they w i l l provide no contr ibut ion to cor re lat ion values. Let us now consider the two s i g n a l s , shown in Figure 3.16, with nega- t i ve D.C. o f fset d, and d 0 , respect ive ly . CD -o Q. E Time - A A M A - (a) -a 3 CL E Time I — W W 1 (b) T Figure 3.16. Signals with D.C. Sh i f ts d^ and d,. Before taking the crosscorre lat ion of these s i g n a l s , they are f i r s t padded with zeros, which is i l l u s t r a t e d in Figure 3^17. 4 CD -a 3 CL Padding (a) ime CD -a E Padding / Time (b) Figure 3.17. Padding Signals with Zeros at the End of Their Time Series 38 When taking the crosscorre lat ion of the funct ions, shown in Figure 3 .17, we e s s e n t i a l l y have the autocorrelat ion of a box car funct ion , 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 for the signals to give representative c rosscor re la t - ion values, they should f i r s t be adjusted so that they are centred at zero o f f s e t . Another way of looking at the D.C. s h i f t problem is to consider Equa- t ion .3.10. An of fset can be represented by s h i f t i n g one signal re la t i ve to the other, that i s X^=Xk+C where C is a constant. In th i s case Equation 3.10 becomes VT> = £Xk Yk+x ( 3 - l l a ) = £ ( x k + c ) Y k + x ( 3- l l b) V T ) s V T ) + ?CYi<^ ( 3 - l l c ) Equation (3.11c) c l e a r l y shows that the D.C. s h i f t would resu l t in a mis- representative crosscorre lat ion value. For a pract ica l example, consider the two wavelets shown in Figures 3.19 and 3.20 with D.C. s h i f t s of -0 .06 and - 0 . 2 , respect ive ly . 39 -I I I I I 0 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 crosscorre lat ion of these two s i g n a l s , we obtain the t r iangle function in Figure 3.21 (with pos i t ive time s h i f t shown). 40 60 80 100 POSTIVF: TIME SHIFT UNITS... Figure 3.21. Crosscorrelat ion of Offset Signals If we correct for D.C. o f f s e t s , the signals take the form shown in Figure 3.22 and 3.23, and the i r crosscorre lat ion function i s shown in Figure 3.24. 60 80 TIME UNIT Figure 3.22. Signal 1 Corrected for 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 for 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 for D.C. Shi f t Corrected Signal? 42 The maximum of th is crosscorre lat ion 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 ve loc i ty p r o f i l e can thus be summarized as fo l lows: (1). Record data (2) Remove D.C. s h i f t (3) Taper time-domain data with cosine be 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 f i l t e r (7) Take inverse FFT (8) Take crosscorre lat ion of f i l t e r e d time-domain s ignals (9) Determine t , time of fset (10) 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 in teract i ve program which displays the frequency spectra , un- f i l t e r e d and f i l t e r e d t races , and crosscorre lat ion plots on a mainframe graphics 43 terminal and has been implemented on the IBM - P.C. The fol lowing subroutines 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 Addit ional subroutines implemented in CROSSCOR are: (10) SUBROUTINE NORME NORME normalizes a data ser ies by i t s RMS energy, in order to bring the contr ibut ion of each trace equal in magnitude when c rosscor re la t ing . (11) SUBROUTINE P0W.2 P0W2 determine M in N=2**M for a s p e c i f i c trace length ( i . e . , L) and then pads the rest of the t race , L to N, with zeros. The value of N is 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 SUBROUTINE FXFFT This subroutine saves computing time,.because i t allows for only one c a l l to FASTF. In FASTF one trace is input as a real component and second trace is input as an imaginary components. Reference 5 further expands on th is (Chapter 3, Section 3.6) 45 CROSSCOR was tested for performance by superposit ioning several s inusoidal waves and f i l t e r i n g out the desired frequencies. In Figure 3.26 there i s an ident ica l traces containing fequencies of 30, 50, and 100 Hertz. The s ignals are o f fset from each other by 20 time un i t s . The signals are then bandpassed between 40 and 80 Hertz resu l t ing in only the 50 Hertz signal remaining as i l l u s t r a t e d in Figure 3.26(c) . The crosscorre lat ion function is i l l u s t r a t e d in Figure 3.26(d) where negative and pos i t i ve time s h i f t s are i l l u s t r a t e d . In determining the a r r i va l time differences between the successive wavelets, i t is only necessary to obtain the maximum pos i t i ve cor re la t ion value at a pos i t i ve time s h i f t . This is due to the fact that the crosscorre lat ion function is com- paring a r r i v a l time differences of the dominant responses ( i . e . , maximum pos i t i ve cor re la t ion value) and the wavelet from the deeper layer w i l l always have a longer a r r i v a l time (pos i t ive time s h i f t ) . The crosscorre lat ion 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 un i t s . 46 Time Units (b) g H •J T O.Ql 1 v 100 200 300 400 500 Frequency, Hz Time Units 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 (cl Bandpassed I races (4U to 80 Hz) (d) Crosscorrelation hunction of the F i l te red Traces. 47 B. F ie ld Program and Discussion of Results In th is section CROSSCOR is 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 ive research s i tes are given, and the ve loc i ty prof i les from d i f fe rent data and methods of data reduction are presented and compared. The research s i tes are i l l u s t r a t e d in Figure B . l and are as fo l lows : 1) McDonald Farm Site 2) Grant McConachie Way Si te at Laing Bridge 3) T i lbury Island Natural Gas Plant S i te 4) Annacis Bridge P i le Research S i te 5) Lower Langley 232nd St . Site on Trans Canada Highway The McDonald Farm Site data was acquired on June 25 and July 2, 1987. In th is invest igat ion f i ve sets of data were compiled and two di f ferent sensors were used. The f i r s t sensor was the Hogentogler Super Cone with a Geospace* ve loc i ty transducer and the second sensor was the UBC Seismic Cone Pressure- meter, which had two (P iezoelectr ic 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 polar i ty method and CROSSCOR. In th is research s i t e , CROSSCOR was found to be very helpful in respect to reducing noisy accelerometer data where picking cross - overs was very d i f f i c u l t or impossible. The Grant McConachie Way Si te data was acquired on November 17, 1987, where the UBC #7 Seismic Cone was used with a P iezoelect r ic bender element accelero- meter sensor i n s t a l l e d . Two di f ferent types of sources were used in th is invest igat ion , 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 accelero - meter observed. The f i r s t response was due to the shear wave (S-wave) occurring at 100 Hz, second response was a compression wave, (P-wave), at 900 Hz , the th i rd response resulted from the resonanting accelerometer at 3000 Hz. The Ti lbury 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 th i s research s i t e i t was not possible to use the reverse po lar i ty method due to the very noisy responses. However, CROSSCOR proved to be extremely helpful in reduc- ing and analyzing the time ser ies . * See Laing (1985). 48 U.S.A. Figure B. l . General Site Location Map 49 The Annacis P i l e Research Si te data was acquired on January 8 , 1988 and February 10, 1988. On January 8, 1988 the UBC Seismic Cone Pressuremeter was used for data a c q u i s i t i o n , and on February 10, 1988 the UBC #8 Seismic Cone was appl ied. In these investigations shear ve loc i t i es were obtained using shear sources. Since both sides of the shear source were used, i t was poss i - ble to obtain two ve loc i ty estimates for one depth increment using CROSSCOR as opposed to the reverse po lar i ty and visual pick method which allows for only one ve loc i ty estimate. The Lower Langley 232nd St. data was acquired on November 19, 1987, where the UBC #7 Seismic Cone was used for data acqu is i t ion . In this invest igat ion f i l t e r e d (analogue) and unf i l tered reverse polarized hammer shear source data was obtained. In add i t ion , 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 fol low in Section B i s as fo l lows: Section B . l : In th i s section the research s i tes investigated are presented along with typ ica l subsurface stratigraphy and cone penetration p ro f i l es which give some indicat ion 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 igni f icance of p ro f i l es presented, one should refer to the numerous UBC i n - s i t u Group reports. Section B.2: In th is section CROSSCOR i s compared to the established reverse po lar i t y method in obtaining ve loc i ty estimates with low^-pass ana- logue 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 th i s invest igat ion i s to prove that CROSSCOR works as well as the reverse po la r i t y method when reducing clean s igna ls , and CROSSCOR works better than the reverse po lar i t y method with signals containing many dominant low frequencies Section B.3: In th is section CROSSCOR is investigated for i t s performance in removing P-waves and noise from unf i l te red accelerometer data. The ve loc i t i es 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 th is section CROSSCOR i s tested for i t s performance in obtaining ve loc i ty 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 fe rent types of sources and giving more ve loc i ty estimates for a downhole seismic p r o f i l e as opposed to the reverse po lar i ty method. Section B.5: In th is section CROSSCOR is tested for i t s a b i l i t y to extract compression wave events from seismic traces and give corresponding com- pression wave v e l o c i t i e s . Section B.6: In th is section specifying the Bandwidths ( i . e . , cutoff f r e - quencies) when applying CROSSCOR is 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 si l t The groundwater table is generally found to l ie lm below ground surface. Additional information can be obtained from Campanella et a l , 1981. * Cone Penetration Test 51 UBC I M SIT U T • E: s ~ r I MC3 S l t a Locat ion! McOONALD FARM CPT Data i 88/09/25 Pago Noi 1 / 1 On S i t e Loci 577-86-2 Cons Unodi HOC SUPER STD PP C o B « a n t « i PIEZO CONE SNO FRICTION RATIO Rf (X) 0 5 SLEEVE FRICTION (ba-) 0 t (A L ai •P Q) E 0_ LU a CONE BEARING Ot (bar) PORE PRESSURE U Ob of »at« r ) 0 100 DIFFERENTIAL P.P. RATIO AU/Ot - .2 0 .8 0 10 20 30 «d soft or rparnc ay loose-d coarse f ine with hd t soft no consoli clayey nse and mal ly lated ; i l t Depth Incramont • . 025 m Max U*pth i 24. OS • Figure B . l . l . Cone Penetration P r o f i l e for McDonald Farm S i te Laing Bridge S i te The Laing Bridge S i te is 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 predict ing settlements. A typical CPT p r o f i l e from the s i te is i l l u s t r a t e d in Figure B.1.2 and has the corresponding strat igraphy: 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 s i l t The groundwater table is generally found to l i e between lm and 1.5m below ground surface, with f luc tuat ion due to t i d a l inf luence. 53 U1 u. 0) E a UBC I M s i n - U T E S " r i NG S i to L o c a t i o n . L a i n g Br- idga S. CPT Data i 10/22/B7 09i 30 Pago Noi 1 / 1 On S l t a L o c i McConach la 0/P Cons Usadi UBC*7 S t d / B F S PP Commantei C P T U - 6 DIFFERENTIAL P.P. RATIO AU/Qt - . 2 0 .8 PORE PRESSURE U (a. of aalar) 0 175 4f> CONE BEARING Qt (bar) SLEEVE FRICTION (bar) 0 ' 1.5 40 FRICTION RATIO Rf (X) 0 5 INTERPRETED PROFILE 40 20 sandy, s i l t y c l a l medium dense to very dens sands normally consolida clayey s i ed t D e p t h Incrament • . 0 2 5 m Max Depth i 54. 775 m Figure B.1.2. Cone Penetration P r o f i l e for Laing Bridge Site Ti lbury Island Natural Gas Plant S i te The Ti lbury Island Natural Gas Plant is located on Ti lbury Is land, 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 typical of a r i ve r delta complex. John Howie (1987-88) conducted extensive work with the UBC Seismic Cone Pressuremeter at th is s i t e . A typical CPT p r o f i l e from the s i t e is i l l u s t r a t e d in Figure B . l . 3 and has the fol lowing strat igraphy: 0-2m F i l l 2-25m Sand 25-45m Interbedded s a n d - s i l t layers The groundwater table is generally found to l i e lm below ground surface, with f luc tuat ion due to t i d a l inf luence. 55 UBC I M S I T U TES B T I M G S l t o L o c a t i o n ! TILBURY IS -LNG CPT D a t a i 31/03/B7 Pago Noi I / 1 On S l t o L o c i TILBURY if 1 Cono Usodi UBC SCP1 Commontoi PPi BT&BS I/) L 01 01 E 0_ Ld a PORE PRESSURE SLEEVE FRICTION U ( » . of wotsr) (bar) 0 50 0 1.5 .0 CONE BEARING Qc (bar) FRICTION RATIO Rf U> 200 0 5 n DIFFERENTIAL P.P. RATIO AU/Oc - . 2 0 .6 15 30 45 INTERPRETED PROFILE interbedded sand -s i l 1ayers D e p t h Inc romant i . 0 2 5 m Max Dopth i 39. 9 m Figure B .1 .3 . Cone Penetration P r o f i l e for Ti lbury Island Si te Annacis P i l e Research S i te The Annacis P i l e Research S i te is located on the extreme east end of Lulu Island which is 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 Min i s - t ry of Transportation and Highways (1984). The s u r f i c i a l geology of the Lulu Island region is s i m i l a r to that of Sea Island and Ti lbury Island ( i . e . , de l ta i c d i s t r ibu tary channel f i l l and marine sediments). A typ ica l CPT p r o f i l e from the s i t e is i l l u s t r a t e d in Figure B . l . 4 and has the corresponding strat igraphy: 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 th in sand layers 57 CO 19 L a a E a. LU a UBC I M SITU TESTING CPT Data i 850813 MO AS DV Poga Noi 1 / 2 Con. U » a d . U0C0 STO TIP C o « « « n t . i NEAfl CASING S l t a Lacat lom ANNA PLT On 3 l t « Loci CPT PR1 PORE PRESSURf SLEEVE FRICTION U U. of «<rtjr) Hw) 0 100 .0 CONE BEARING Ot (bar) 10- 20 2.3 0 0 FRICTION RATIO DIFFERENTIAL P.P. Rf CD RATIO lU/Ot 200 0 3 . _0 . ' INTERPRETED PROFILE . • 10" - 20- 30 SAND f i l l soft oraanic s i l t y CLAY med. dense SAND minor s i l t y SAHp l e n s e s spe o v e r Dopth I n c r « » « n t i .025 n Max Dap th i 35.075 - Figure B.1 .4. Cone Penetration P r o f i l e for Annacis Bridge P i l e Research Si te Lower Langley 232nd St . The Lower Langley 232nd St . S i te i s located adjacent to the northern ex i t road on the westbound lane of the main Trans Canada Highway out of Vancouver which follows the Fraser Va l ley . The e x i t road is the of f ramp to 232nd Street . This s i t e has been investigated for geotechnical and geological parameters by John Sul ly (1987) and James Greig (1985). The s u r f i c i a l geology of th is s i t e is typ ica l of the Fraser Lowland area. A CPT p r o f i l e from the Lower Langley 232nd St . S i te is i l l u s t r a t e d in Figure B . l . 5 and consist predominantly of normally consolidated sens i t i ve s i l t y clay with occasional sand lense. The near surface material is overconsolidated due to dess icat ion . 59 UBC IN SITU TEST IMG S i t e L o c a t i o n . R G C - E 8 - A B On S i t e L o c i L0WER232 LANCLEY CPT D a t a i 1 1 - 1 9 - B 7 1 6 i J 5 Cong Usedi HOC SUPER STD u Page Not 1 / I C o m m e n t « i C 7 7 - 8 7 1 3 5HHF1LT FRICTION RATIO Rf (Jt) 0 5 SLEEVE FRICflON (bar) 0 2 CONE BEARING Oc (bar) PORE PRESSURE DIFFERENTIAL P.P. U (a. of later) RATIO lU/Qc INIERPREtEQ PROFILE Depth Inc rement i . 0 2 5 m * Overconsolidated Max O e p t h • 29 . 655 m f i l l o .c . c lay N.C. c lay N.C. clay with f ine interbeddn N.C. c lay Interbedqed clay and snad sand **Normany consolidated Figure B . l . 5 . Cone Penetration P r o f i l e for Lower Langley 232nd St . S i te B.2 Cross-Over Method vs. CROSSCOR with Analogue F i l t e r e d Applied and Dif ferent Seismic Probes Used In th is section the well establ ished reverse po la r i t y or Cross-Over method is 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 iezoe lec 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 I I , the Cross-Over method demands that in picking a Cross-Over, a consistent method must be used throughout the seismic p r o f i l e . 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 is futher i l l u s t r a t e d in Figure B.2.1 with data recorded at Annacis Is land. In the discussion to fo l l ow , 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 ident ica l v e l o c i t i e s to that of CROSSCOR with clean s ignals ( 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 is from McDonald Farm, where the Hogentogler Super Cone ve loc i ty transducer was used for data acqu is i t ion 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 Cross-Over method and CROSSCOR. As is shown, the v e l o c i t i e s are very comparable due to the clean signals the slow responding ve loc i ty transducer gives. One must keep in mind that since CROSSCOR uses a l l the information contained in the seismic t races , 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 th is set of data is 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 in ve loc i ty with increase in bearing). The second set of data presented is also from the McDonald Farm S i t e , where the UBC Seismic Cone Pressuremeter (SCPM) was used for data acqu is i t ion on June 25, 1987. In th is invest igat ion the signals recorded had several dominant f r e - quencies present even with the 300 Hz low-pass f i l t e r appl ied . This character - 61 Time(msec) Figure B . 2 . 1 . Seismic Data Acquired from Annacis Island Vibro-Compaction S i te on May 18, 1988. Data Recorded with Accelerometer having a 300 Hz Low- Pass F i l t e r Appl ied. 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. Average Depth (m) 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 63 Velocity ( m/s ) 100 200 • t i I i — i — i — I — h Cone Bearing Qt (bar) (b) 0 100 200 «~~~~ CROSSOVER CROSSCOR CL a • i i i — i — i — l Figure B . 2 . 2 . (a) Veloci ty 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 Sel f F i l t e r i n g (Natural Frequency = 28 Hz) Geo- phone (b) Cone Bearing P r o f i l e from McDonald Farm. 64 i s t i c of these seismic traces is typical of the accelerometer transducer due \ to i t s high s e n s i t i v i t y . The ve loc i t i es determined from the Cross-Over method and CROSSCOR are given in Table B .2 .2 . The SCPM has accelerometers in both the cone portion and the pressuremeter portion of the instrument. The cone accelero- meter (SCPM CN) gives noticeable discrepancies between the two methods of ve lo - c i t y ca lcu la t ion . But when these ve loc i t ies are compared to the SCPM pressure- meter accelerometer (SCPM PM), the ve loc i t ies obtain from CROSSCOR are more s i m i l a r . The ve loc i t i es obtained from th is research s i t e are plotted in Figure B.2.3(a) . 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 is shown, the v e l o c i t i e s in Figure B.2.3(a) vary with the cone bear- ing p r o f i l e . Laing Bridge S i te The th i rd set of data presented is from the Laing Bridge S i t e , where the data was acquired on November 17, 1987 with the UBC #7 Seismic Cone. For th i s data i t was also not possible to accurately determine consistent Cross-Overs, due to the sensit ive nature of the accelerometer (as discussed in Section I ) . The shear wave ve loc i t i es 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 ve loc i ty determination (note depths 6 . 3 , 10.3, 11.3 and 13.3m), but when compared to the ve loc i t i es obtained from the P-Plate source at 10.3 and 11.3m, CROSSCOR gives the most repeatable values. The ve loc i t ies obtained from th is s i te are plotted in Figure B.2.4 along with the cone bearing p r o f i l e . The las t set of data presented in th is section was obtained from the Lower Lang ley 232nd St . on November 19, 1987, where the UBC #8 Seismic Cone was used for data acqu is i t ion . In th is invest igat ion very clean signals were recorded with pre- dominantly one S-wave frequency of 80 Hz present. The calculated ve loc i t ies from the Cross-Over method and CROSSCOR are given in Table B .2 .4 . As is i l l u s - t rated, the ve loc i t i es obtained are in very close agreement. The ve loc i t ies are also plotted in Figure B.2.5 along with the cone bearing p r o f i l e . 65 Table B.2.2 Calculated Velocity P ro f i l es 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 iezoelect r ic Bender Accelerometers. Cross-Over CROSSCOR CROSSCOR June 25, 1987 June 25, 1987 June 25, 1987 Average Depth (m) Depth (m) Interval SCPM - CN Vs (m/sec) SCPM - CN Vs (m/sec) SCMP - PM Vs (m/sec) 3.58 3.08-4.08 150 164 4.58 4.08-5.08 168 158 5.45 4.95-5.95 154 6.58 6.08-7.08 95 185 7.45 6.95-7.95 214 7.58 7.08-8.08 153 188 8.45 7.95-8.95 189 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 66 Velocity ( m/s ) (t>) Cone Bearing Qt (bar) 0 100 200 0 100 200 Figure B . 2 . 3 . (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 Appl ied, (b) Cone Bearing P r o f i l e from McDonald Farm. 67 Table B.2.3 Calculated Velocity P ro f i l es from Laing Bridge, 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 Appl ied. In t h i s Set of Data i t was not Possible to Obtain Accurate Cross-Over P icks . Seismic Cone Seismic Cone Seismic Cone Shear Wave Source Shear Wave Source P-Plate Source Average Cross-Over CROSSCOR CROSSCOR Depth (m) Vs (m/sec) Vs (m/sec) 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 11.3 210 150 12.3 145 144 13.3 120 197 14.3 193 183 15.3 197 224 16.3 246 189 152 68 Figure B .2 .4 . Veloci ty 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 Appl ied. In th is Set of Data i t was not Possible to Obtain Accurate Cross-Over P icks , (b) Cone Bearing "Prof i le from the Laing Bridge"sTte. 69 Table B.2.4 Calculated Veloci ty P ro f i les 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 Appl ied. Cross-Over CROSSCOR Hammer Source Hammer Source Average Depth (m) F i l t e r On F i l t e r On Depth (m) Interval Vs (m/sec) Vs (m/sec) 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 1 1 1 I 1—1 30 Figure B.2.5. Veloci ty 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 Appl ied. (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 te red (Analogue) and Unf i l tered ( i . e . , Removing P-Waves and Noise With CROSSCOR) Seismic Traces In th is section CROSSCOR is applied to unf i l te red 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 f i r s t response is the shear wave located at 100 Hz, the second response is the com- pression wave at 800 Hz and the th i rd corresponds to the accelerometer resonating at 3000 Hz. From th is t race , i t was possible to i so late the shear wave by apply- ing a bandpass of 80 to 120 Hz, as is i l l u s t r a t e d in Figure B . 3 . 2 . The above extract ion of the shear wave was possible with most of the unf i l te red seismic traces recorded and the se lect ion of the bandwidth was not a c r i t i c a l c r i t e r i o n . 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 with in the shear wave. In most of the data reduced, determination of the ve loc i ty was rout ine , because one only had to specify the same bandwidth throughout the p r o f i l e . Lower Langley 232nd St . The f i r s t set of data is from the Lower Langley S i t e . In th i s invest igat ion the shear wave was located at 85 Hz, and the applied bandwidth was 70 to 100 Hz. The e f fect of time s h i f t s resu l t ing from the 300 Hz low-pass analogue f i l t e r was determined by c rosscorre lat ing the f i l t e r e d and unf i l te red t races. From th is procedure, i t was found that time s h i f t s were n e g l i g i b l e . The calculated ve lo - 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) i s a plot of the v e l o c i t i e s calculated from CROSSCOR and the reverse po la r i t y method, and the cone bearing p r o f i l e is given in Figure B.3.3(b) . The Figure i l l u s t r a t e s that the unf i l te red and f i l t e r e d resul ts are nearly ident ica l and vary along with the cone bearing p r o f i l e . Thus, in th is set of data CROSSCOR separated out the shear waves from the P-waves and noise, 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 lbury Island Natural Gas Plant S i te The second set of data is from the T i lbury Natural Gas Plant S i t e , where i t was not possible to pick Cross-Overs. Therefore, CROSSCOR was found to be a invaluable tool in obtaining the shear v e l o c i t i e s . In th is case the shear wave 72 1.0 3.8 Metre Seismic Trace (3_ - 1 . 0 0.0 TIME (msec) (a) Seismic Trace Recorded at Laing Bridge with the P-Plate Source (Unf i l te red) . 1.0 S (100 Hz) M j = 3 i—• P (800 Hz) !3 , | •CC \L J 3000 FREQUENCY HERTZ 5000 (b) 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 6.8 metre •1.0 3.8 metre TIME (msec) j TIME (msec) F I L T E R E D (80-120 Hz) 6.8 metre FILTERED (80-120 Hz) ' ' 1.01 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 te red with a Bandpass of i 80 to 120 Hz Inorder to Extract the Shear Wave. 0.0 CROSSCORRELATION 200 TIME SHIFT (msec) (b) Corresponding Crosscorrelat ion 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 te red (300 Hz Low-Pass) and Unf i l tered Seismic Data. The Data was Acquired with an Accelerometer. Cross-Over Hammer Source Average Depth (m) F i l t e r On Depth (m) Interval Vs (m/sec) 2.2 1.7-2.7 82 3.2 2 .7 -3 .7 92 4.2 3 .7 -4 .7 103 5.2 4 .7 -5 .7 104 6.2 5 .7 -6 .7 96 7.2 6 .7 -7 .7 98 8.2 7 .7 -8 .7 108 9.2 8 .7 -9 .7 104 10.2 9.7-10.7 109 11.2 10.7-11.7 117 12.2 11.7-12.7 119 13.2 12.7-13.7 118 14.2 13.7-14.7 128 15.2 14.7-15.7 128 16.2 15.7-16.7 137 17.2 16.7-17.7 127 18.2 17.7-18.7 134 19.2 18.7-19.7 149 20.2 19.7-20.7 140 21.2 20.7-21.7 139 22.2 21.7-22.7 155 23.2 22.7-23.7 162 24.2 23.7-24.7 177 25.2 24.7-25.7 174 26.2 25.7-26.7 184 27.2 26.7-27.7 213 28.2 27.7-28.7 189 * => not possible to pick Cross-Overs 75 Cross-Over CROSSCOR CROSSCOR Hammer Source Hammer Source Hammer Source F i l t e r Off F i l t e r On F i l t e r Off Vs (m/sec) Vs (m/sec) Vs (m/sec) * 83 * 101 101 * 104 109 * 105 99 * 96 96 * 99 100 * 110 110 * 104 111 * 103 104 * 113 120 * 123 120 * 123 120 * 124 127 * 130 133 * 134 138 * 130 124 * 131 141 * 138 134 * 150 160 * 135 131 * 160 160 * 151 155 * 172 172 * 184 185 * 171 168 * 192 199 * 191 198 because of varying s ignals . Figure B . 3 . 3 . (a) Velocity P r o f i l e from Lower Langley Comparing F i l te red . (300 Hz Low-Pass) and Unfi l tered 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 resul t ing ve loc i t ies 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 ve loc i t i es vary with the cone bearing p r o f i l e except for the ve loc i t i es at 3.66 and 6.64 metres. This maybe due to poor responses being recorded at these depths. 77 Table B.3.2 Calculated Veloci ty P r o f i l e s from Ti lbury S i te Comparing Shear Wave Ve loc i t ies Obtained with CROSSCOR from Seismic Data Acquired by Two Different Accelero- meters. In th is Set of Data i t was not Possible to Obtain Accurate Cross- Over P icks . SCPM-CN SCPM-PM Seismic Cone Pressuremeter Average Accelerometer Accelerometer Depth (m) Vs (m/sec) 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 Velocity ( m/s ) 0 80 160 i i i I i i i I (a) Cone Bearing Qc ( bar ) ***** SCPM CN ooooo SCPM PM o o * J L CL Figure B .3 .4 . Calculated Velocity P ro f i l es from Ti lbury Site Comparing Shear Wave Veloc i t ies Obtained with CROSSCOR from Seismic Data Acquired by Two Different Accelerometers. The Veloci ty P ro f i l es are also I l lus t ra ted with the Cone Bearing P r o f i l e . In th is Seismic Data i t was not Possible to Obtain Accurate Cross-Over P icks , (b) Cone Bearing P r o f i l e from the Ti lbury Island S i t e . 79 B.4 CROSSCOR with Nonpolarized Sources In this section the appl icat ion 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 . , r ight and l e f t ) . The Buffalo Gun source is thoroughly discussed in Reference 2 by Laing (1985) with regard to design and use. The Buffalo Gun is desirable because i t can be used both on and o f f shore to obtain both shear and compression wave ve loc i t ies w i th - out the problems associated with a seismic cap or explosive source. CROSSCOR is applied to data obtained from s t r i k i n g each Hammer Shear Source in one d i rec - t i o n , therefore, i t allows one to obtain two veloci ty 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 ight s ide) . This is opposed to the reverse po lar i t y method which allows for only one ve loc i ty 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 is i l l u s t r a t e d in Figure B .4 .1 . From th is f igure , one can c lear l y d ist inguish the compression and shear waves. The shear waves.were obtained by bandpassing the data from 90 to 150 Hz. The resul t ing shear wave ve loc i t i es 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 . Some of these values are comparable, while others deviate considerably. This deviation was found to resu l t from tr igger ing problems with the Buffalo Gun. The Buffalo Gun source is tr iggered with a geophone which may resul 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 resul t ing 2.71 msec s h i f t was recorded. For a Vs of 100 m/sec the interva l time is 10 msec over a 1 metre distance and 2.71 msec error represents +27% error on measured Vs. If the t r igger error were a constant 2.72 msec then the potential maximum error increases as Vs increases. The measured t r igger 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 igger ing 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 imi la r to the shear p late ) . 80 0 20 40 60 80 100 120 140 TIME {m sec) Figure B . 4 . 1 . Seismic Section from Lower Langley Where a Buffalo Gun was Used as the Source and An Unf i l tered Accelerometer Recorded the Data. 81 00 cd C C/3 i - l O > X) V •1-1 1-1 0. low-pass analogue f i l t e r a pplied 0 20 40 60' 80 TIME (msec). 100 120' 140 Figure B .4 .2 . Two Buffalo Gun Seismic Traces from Lower Langley I l l u s t r a t i n g t r igger ing Prop ems with the Buffalo Gun Source. The Data was Acoui reT^Tr i r t r , 1 Accelerometer at a Depth of 13 Metres. The Result ing Time W f f e r i H F e was Table B.4.1 Calculated Veloci ty P ro f i les of Lower Langley (232) Comparing Shear Wave Veloc i t ies Obtained from the Hammer Shear Source and the Buffalo Gun Source. The Data was Acquired with an Accelerometer (Unf i l te red) . CROSSCOR CROSSCOR Hammer Source Buffalo Gun Depth (m) F i l t e r Off F i l t e r Off % Difference Interval Vs (m/sec) Vs (m/sec) (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 igger error creates two incorrect successive measurements of Vs. 83 Annacis Bridge P i l e Research S i te The second set of data is from the Annacis P i l e Research s i t e , where velo - 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 ve loc i t i es 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 ve loc i t i es are expected for th is very soft organic s i l t and tend to l inear increase with depth as the i n - s i t u stresses increase ( i . e . , essent ia l l y constant cone bearing p r o f i l e which has been corrected for overburden stresses) . The higher ve loc i ty values at 12.0 metres are ind icat ive of the variable ( in depth) sand unit at approximately 15.0 metres. The ve loc i t i es 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 p r o f i l e . A special note must be made with respect to the cone bearing p ro f i l es shown in Figures B.4.3 and B .4 .4 . These two prof i les are the same and are from a d i f ferent invest igat ion ( i . e . , d i f ferent CPT locations) than the ve loc i t ies shown in Figures B.4.3 and B.4 .4 . Thus, the pro f i les 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 fe rent CPT locat ions ) . The ve loc i t ies shown for the Annacis P i l e Re- search S i te are very comparable in each set of data and give greater confidence in the resu l t s . 84 Table B.4.2 Calculated Veloci ty P r o f i l e s for Annacis P i l e S i t e , where Ve loc i t ies Obtained from St r ik ing each Side of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer (Unf i l te red) . January 8 , 1988 January 8 , 1988 Depth (m) Pressuremeter Right Side Vs (m/sec) Pressuremeter Left Side Vs (m/sec) % Difference (Reference-Riqht Side) 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 Velocity ( m/s ) Cone Bearing Qt (bar) 100 (b) 0 * * * * * CROSSCOR S C P V - P U Rtciht S M « acsaB CROSSCOR S C P M - P M un SJ6« Figure B .4 .3 . (a) Calculated Veloci ty Prof i les from the Annacis P i l e Research S i t e , where Ve loc i t ies Obtained from both Sides of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer (Unf i l te red) . (b) Cone Bearing P r o f i l e of the Annacis Bridge P i l e Research S i te 86 Table B.4.3 Calculated Veloci ty P ro f i l es for Annacis P i l e S i t e , where Ve loc i t ies Obtained from St r ik ing each Side of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer (Unf i l tered) . CROSSCOR CROSSCOR January 8 , 1988 January 8 , 1988 Cone #8 Cone #8 Right Scale Left Scale % Difference Depth (m) Vs (m/sec) Vs (m/sec) (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 Prof i les from the Annacis P i l e Research S i t e , where Ve loc i t ies Obtained from both Sides of the Hydraulic Pads of the Hammer Shear Source are Compared. The Data was Acquired with an Accelerometer (Unf i l tered) . (b) Cone Bearing P r o f i l e of the Annacis Bridge P i l e Research S i te 88 B.5 CROSSCOR in Obtaining Compression Wave Ve loc i t ies In th is section the appl icat ion of CROSSCOR to obtaining compression wave v e l o c i t i e s is addressed. As was discussed previously , the seismic traces obtained from unf i l te red accelerometer data contain information of both shear and compression body waves. The extent of th is information large ly depends on sources used. 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 so la te the indiv idual body wave responses, i t is necessary to obtain the frequency spec- trum 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-Plate source at Laing Bridge S i te at a depth of 6.8 metres. From the frequency spectrum obtained, one can ident i f y three dominant responses. The f i r s t response is the shear wave occurring at 150 Hz, the second response is the compression wave at 900 Hz, and the th i rd 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 largely 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 is i l l u s t r a t e d in th 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 ve loc i ty must be at least greater than that for water or about 1500 m/sec. Table B.5.1 shows that typ ica 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 St . 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 th is data (as was shown in Section B.4) was that t r igger ing was inconsistent . The compression wave v e l o c i t i e s obtained from the Buffalo Gun data are shown in Table B.5.2 where a bandpass of 200 to 1000 Hz was appl ied . The nonsensical values obtained are most l i k e l y a resul t from t r igger ing problems. 89 INPUT HAUELET o.o TIME (msec) (a) Seismic Trace from Laing Bridge where a P-Plate Source was Used and an Accelerometer (Unf i l tered) Recorded the Data at 6.8 Metres. 1.0 FREQUENCY SPECTRUM OF TRACE •-3 3000 FREQUENCY HERTZ 5000 (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 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 . Veloci ty Versus Density for 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 Ve loc i t ies in Sediments and Rocks (Milton B. Dobrin, 1976) Compressional velocity Shear velocity Material and Source m/s ft/s 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 Sandstone Sandstone conglomerate, Australia 1400-1300 2400 4620-14,200 7920 Limestone: Soft Solenhofcn, Bavaria Argillaceous, Tex. Rundle, AJberta 1700-1200 5970 6030 6060 5610-13,900 19,700 19,900 20,000 2880 3030 9500 10,000 Anhydrite, U.S. Midcontinent, Gulf Coast 4100 13,530 Clay Saturated Loose sand Saturated 1100-2500 1800 3630-8250 5940 500 1,650 SOURCE: Sydney V. Clark, Jr. (ed.), "Handbook of Physical Constants," rev. ed., Geol. Soc. Am. Mem. 97, 1966. 92 Table B.5.2 Calculated Veloci ty P r o f i l e from Lower Langley, where an Accelerometer was Used for Data Acqu is i t i on . The Compression Wave Veloc i t ies were Obtained from the Buffalo Gun Source and by Applying a Bandpass of 200 to 1000 Hz. Depth (m) CROSSCOR Buffalo Gun F i l t e r Off Vp (m/sec) 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 ve loc i t ies due to poor and/or uncorrelated responses. 93 I l lus t ra ted in Figure B .5 .3 . are traces obtained from the Buffalo Gun data along with the isolated compression wave responses. This f igure i l l u s t r a t e s that the compression wave responses are c lear l y i d e n t i f i a b l e and correlated ( i . e . , s i m i - l a r ) . Thus, one can conclude that i f t r igger ing was repeatable, accurate com- pression wave ve loc i t ies could be obtained. The next set of data was from Lower Langley where the Hammer Shear was used. The compression waves ve loc i t i es obtained from this set of data (unf i l tered accelerometer) are shown in Table B.5.3 and were obtained with a bandpass of 750 to 1000 Hz appl ied. The nonsensical ve loc i t ies are most l i k e l y due to the fact that c lear compression wave responses were not obtained. Figure B.5.4 i l l u s - trates that the iso lated compression wave responses are very noisy and poorly corre lated. 94 4.7 metres 5.7 metres 0.0 TIME (msec) TIME (msec) 1.0 4.7 metres (200-10Q0 Hz) CD 3 CL E - 1 . 0 5.7 metres (200-1000 Hz) TIME (msec) TIME (msec) 1 0 0 11.7 metres TIME ( msec) 12.7 metres 0,0 TIME (msec) 1 0 0 1.0 11.7 metres (200-1000 Hz) CD CL E <: •1.0 j P-Wave 12.7 metres (200-1000 Hz) 0.0 IME(msec) 1 0 0 T T UV 1 TIME (msec) 1 0 0 Figure B . 5 . 3 . Seismic Traces Obtained at Lower Langley where a Buffalo Gun Source was Used. The Traces were F i l te red (200 to 1000 Hz) in order to Extract the Compression Wave Responses The resu l t ing Compression Wave Responses are Highly Correlated. g 5 Table B.5.3 Calculated Velocity P ro f i l es from Lower Langley, where an Accelerometer was Used for Data Acqu is i t ion . The Compression Wave Veloc i t ies were Obtained by Applying a Bandpass of 750 to 1500 Hz. CROSSCOR Hammer Source F i l t e r Off Depth (m) Vp (m/sec) 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 2448 10.7-11.7 * 11.7-12.7 * 12.7-13.7 4937 13.7-14.7 * 14.7-15.7 1651 15.7-16.7 * 16.7-17.7 * 17.7-18.7 * 18.7-19.7 * 19.7-20.7 * 20.7-21.7 * 21.7-22.7 1244 22.7-23.7 * 23.7-24.7 * 24.7-25.7 * 25.7-26.7 1246 26.7-27.7 * 27.7-28.7 * * => nonsensical compression wave ve loc i t ies are due to poor and/or uncorrelated responses Q ( . 1.0 cu -o 3 +-> - 1 . 0 9.7 metres T i n t (msec) 0.0 TIKE (msec) 1.0 8.7 metres (750-1000 Hz) g 7 m e t r e s ( 7 5 0 _ i 0 0 0 Hz) cu -o 3 -t-> CL E - 1 . 0 °-° TIKE (msec) 200 •o 3 -t-> Q- E - 1 . 0 °-° TIKE (msec) 200 14.7 metres 15.7 metres Figure 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 Result - ing Compression Waves are Poorly Correlated. 97 B.6 Bandwidth Considerations when Applying CROSSCOR In Section IIIA, i t was stated that the crosscorrelat ion function provides information about the s i m i l a r i t y and time difference between two successive wavelets. In th is case we are looking for the a r r i va l time difference between the dominant responses of the recorded wavelets. The main idea behind CROSSCOR is to i so late spec i f i c body waves in the frequency domain and then crosscorre- late the f i l t e r e d wavelets in the time domain. The shear and compression body waves are ident i f i ed by the 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 - ter ing the most important considerations are to : (1) Include dominant frequency character iz ing wavelet, because the other contributing frequencies do not effect the dominant responses s i g n i f i - cant ly . (2) The bandwidths spec i f ied must be consistent throughout the seismic p r o f i l e in order that the crosscorrelated dominant responses are s imi la r in transient response speci f icat ions ( 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 this frequency spectrum the dominant shear wave response is i den t i f i ed 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 th is seismic wavelet, the dominant shear wave res- ponse is also ident i f i ed 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 the 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 band- passed 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- cor re lated , the two d i f ferent sets of f i l t e r e d data gave the same time s h i f t and corresponding ve loc i ty 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 is 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. From Figure B.6.4 the dominant shear wave responses are from 50 to 200 Hz. Figure B.6.5 i l l u s t r a t e s 98 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 ( msec) (a) Seismic Trace at 3.7 Metres Acquired with An Unf i l tered Accelerometer. 50 500 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 . 2 . Seismic Trace Recorded at Lower Langley Where a Hammer Shear Source was Used. 100 2.7 metres 3.7 metres •0 TTUP i \ 200 0.0 T T i i t v \ 200 TIME (msec) TIME(msec) 2.7 Metres (20-200 Hz) 3.7 Metres (20-200 Hz) - 1 . 0 °-° Tiyp / \ 200 0.0 TTUP/' ^ 200 l int (msec) TIME (msec) Figure B . 6 . 3 . (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 Veloci ty 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 ( m s e c ) (a) Seismic Trace at 3.7 Metres Acquired with an Unf i l tered Accelerometer. S-Wave 50 200 300 FREQUENCE HERTZ 400 500 (b) 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. Figure B .6 .4 . Seismic Trace Recorded at Lower Langley Where a Buffalo Gun Source was Used 102 4.7 Metres 0.0 TIME (msec) (a) Seismic Trace at 4.7 Metres Acquired with an Unf i l tered Accelerometer. S-Wave 200 FREQUENCY HERTZ 500 (b) Frequency Spectrum I l l u s t r a t i n g Dominant Shear Wave Response Situated Between 50 to 200 Hz. Figure B .6 .5 . 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 the 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- cor re lat ing these two d i f fe rent sets of f i l t e r e d data resulted in the same time s h i f t and veloc i ty of 94 m/sec. Figure B.6.7 shows the same Buffalo Gun Traces bandpassed between 100 to 200 Hz. The resul t ing veloc i ty (93 m/sec) obtained from th is set of data was negl ig ib ly d i f ferent from those previously obtained. 104 3.7 Metres 1.0 4.7 Metres m •a 4-> CL E - 1 . 0 1.0 cu -a 13 •4-> 1 CL E - 1 . 0 0.0 TIME (msec) 160 0.0 TIME (msec) 160 1.0 3.7 Metres (60-180 Hz) 4.7 Metres (60-180 Hz) cu -a CL E •1.0 °-° TIME (msec) 160 0.0 TIME (msec) Figure B .6 .6 . (a) Seismic Traces Recorded at Depths of 3.7 and 4.7 Metres, (b) Trace in (a) F i l te red Between 20 to 200 Hz Resulting in a Veloci ty of 94 m/sec. (c) Traces in (a) F i l te red Between 60 to 180 Hz Resulting in a Veloci ty of 94 m/sec. 105 3.7 metres (100-200 Hz) 4.7 metres (100-200 Hz) °'° TIME (msec) 1 6 0 0-0' T T M r , ^ 160 l int (msec) Figure B .6 .7 . Trace in Figure B.6.6(a) F i l te red Between 100 to 200 Hz Resulting in a Veloci ty of 93 m/sec. 106 B.7 CONCLUSIONS The resu l ts from the previous discussions suggest that CROSSCOR i s a very benef ic ia 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 - gat ions, 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 fe rent 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 reduct ions, i t was found that the d i f fe rent 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 is 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 t imes, 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 fe rent 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 com- pression 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 igger ing was not consistent thorough out a seismic p r o f i l e resu l t ing in erroneous ve loc i t y c a l c u l a t i o n s . It was suggested that the t r igger ing problem could be resolved by using an accelero- meter (having a fast response time) instead of the geophone as well as gound- ing the Buffalo Gun f i r i n g pin to create a contact closure ( s imi la r to the shear p la te ) . The performance of CROSSCOR is summarized in the fo l lowing l i s t of advantages/disadvantages with comparisons made to the reverse po la r i t y method in the l i s t of advantages and comparisons made to the fol lowing 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 extract ing shear and compression waves. (2) The crosscorre lat ion function in CROSSCOR makes use of a l l the information in the s ignals (averaging out i r r e g u l a r i t i e s and putting s ign i f i cance or dominant responses) as opposed to the reverse po la r i t y method which r e l i e s only on one cross-over point to determine interva l t imes. 107 (3) CROSSCOR works as well as the reverse po la r i t y method when reducing clean s i g n a l s , and CROSSCOR works better than the reverse po la r i t y method with s ignals containing many dominant low frequencies ( e . g . , Figures 2.8 and B .2 .1) . (4) CROSSCOR can obtain ve loc i ty 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 po la r i t y method in working with d i f fe rent types of sources and giv ing more ve loc i ty estimates for a downhole seismic p r o f i l e as opposed to the reverse po la 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 end ly program which has f l e x i b i l i t y in speci fy ing the cutoff frequencies i f there is a repeatable t r igger of fas t response and s u f f i c i e n t accuracy and reso lu t ion . Pi sadvantages (1) CROSSCOR cannot give s o i l damping propert ies . (2) CROSSCOR cannot give any ind icat ion 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 br ief descript ion 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 ar r i va l times. More deta i ls of the KF equations, with examples and appl icat ions to geophysical problems, can be found in Reference 24. - The KF should be investigated for i t s possible appl icat ion to obtaining accurate estimates of the P-wave and S-wave amplitudes and a r r i va l times, for the par t i cu lar 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 imi tat ions of the frequency domain approaches that were discussed in Section I I I . The Kalman F i l t e r is 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 th is 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 condit ions, 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 er rors , and s t a t i s t i c a l information about the i n i t i a l condi - t ions. Figure 4.1 i l l u s t r a t e s the basic re lat ion between the system, the measurements, and the Kalman F i l t e r , for appl icat ion 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 er rors , and the apr ior 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 inear time-varying systems, and with non-stationary system and measurement s t a t i s t i c s . Problems with non l inear i - t ies can sometimes be handled by l inear i z ing the system and measurement equa- t ions . Furthermore, the Kalman F i l t e r is readi ly applied to est imation, smoothing, and pred ic t ion . One of the most important aspects that is real ized when applying the KF to solving a spec i f ic problem, is that a considerable e f fo 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 or difference equations. In addi t ion , the s t a t i s t i c a l properties of the system 109 and measurements need to be modelled by first and second order statistics (expected values and covariances). It should also be noted that satisfactory results are not always obtained by a specific formulation of a problem. It may be necessary to try different approaches before satisfactory results are obtained. For example, different smoothing techniques may need to be considered to improve the estimates. System Measurement Apriori Error Sources Error Sources Statistics System State ' x(t) , Observations 1 Estimate of System State System MaaCM romant z(t) Kalman x(t) Filter P(t) x(t) i 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 equations for a typical acceleration measuring device. 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 charac te r i s t i c natural frequency and damping. The modelling of a seismic receiver resu l ts in 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 ca l ib ra t ion when applying his Acoustic Sub-Seabed inter rogator , he s t a t e s , "The spectral response of the transmitter and receiver is also being pursued to ref ine the i r overal l 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 for rapid measurements of attentuation in marine sediments." The seismic receiver acts as a transform function in respect to the input of noise and body waves propagating through the earth. This in teract ion is i l l u s t r a t e d in Figure 4.2 when we have no ise , state x (5) , and signal (amplitude ^x( l ) and a r r i v a l time ^x(2)) in teract ing with the instrument which is repre- sented by the second order system F(s) with natural frequency w , and damping ? . The output i s a superposition of the no ise , x (5 ) , passing through the instrument and giv ing 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) resu l t in the second order responses x(8) and x(9) . Therefore, knowing the output, transform function ( instrument), i t i s desired to determine the characte r i s t i cs of the input signal ( 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 in Appendix C, The noise, 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 ef fect 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 is 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 ef fect an increasing time constant has on a white noise input. As shown, an increase in Tc tends to smooth the signal ind icat ing 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.. I l l Acceleration Input N o i s e % P - w a v e * \ ^ x^=arrival time X2=arrival time Geophone G(s) = ~2 5 + 2 * V + u o Figure. 4.2. Schematic I l lustration of the States, x̂ for the KF Formulation 112 SYNTHETIC DATA:INPUT NOISE I I I 1 i _ 0.0 0.5 1.0 1.5 2.0 Time (sec) Figure 4 .3 . The Effect An Increasing Time Constant Has On White Noise (Note: T = .003, 005, . 0 1 , .04, .07 , .09 , . 1 , .2) c 113 Time(.sec) f*° Time (sec) Input Noise T = .01 Geophone Response Time (sec) Time (sec) Figure 4 .4 . 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 is i l l u s t r a t e d , the impulses resu l t in an i n i t i a l displacement than decay which is a function of the damping and natural frequency of the instrument. Also shown is 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 is i l l u s t r a t e d in Figure 4 .6 . The impulse was generated by the fol lowing 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 Input Noise Geophone Response Figure 4 .5 . Geophone Response Due to P and S-Wave Impulses and Noise with Increasing Variance 117 Input Noise Geophone Response Time (sec) Time (sec) Figure 4 . 5 . 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 SO 100 150 200 TIME (N-200.DELT-.01,TC=. 10.VARlANCE-.09.R0-.ia) SYNTHETIC DATA: Z 50 100 ISO TIME (N=200.DELT=.0i.TC=. 10,WO=8) 200 Figure 4 .6 . I l l u s t r a t i n g the Synthetic (a) Input and (b) Geophone Response 119 The signal i s defined by a 50 Hertz cosine wavelet exponentially decay- ing at .01 and having a maximum amplitude of 2500. The signal is i n i t i a t e d at 66 time units and ends at 76 time units with the to ta l trace last ing for 200 time uni ts . The second order response is 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 amplitude, and x(2) , the signal a r r i v a l time. In the KF i t i s necessary to give apr ior i estimates of the above two s ta tes , these i n i t i a l estimates were x(l)=2500 and x(2)=75 time un i ts . 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 un i t s ) . In Figure 4.7 we have eight i te rat ions or passes made through the f i l t e r , with each subsequent i te ra t ion having i t s i n i t i a l state value given by the last i te ra t ion where the mean square e r r o r , Figure 4 . 8 , was minimum. The resul ts show that the a r r i v a l t ime, state x (2 ) , approached the actual value (66 time units) and reached i t af ter 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 resul ts are s i g n i f i c a n t , because if. th is f i l t e r can be f ine l y tuned , i t would be possible to obtain a r r i va l times and amplitudes on l i ne resul t ing 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 , respect ively . In add i t ion , i t would also be possible to model and ca l ib rate 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 fo r t needs to be made to ver i f y the numerical results and models. In a d d i t i o n , i t w i l l be necessary to invest igate smoothing techniques to obtain reasonable estimates of the P-wave and S-wave parametres x ( l ) to x(4). Since the ef fects of the P-wave and S-wave wavelets pass through measure- ment interval i s a r e l a t i v e l y short t ime, recognition of th i s fact to derive a suitable smoothing version of the KF needs to be made. In addi t ion , further invest igat ion needs to be made to optimize the method of working mult iple passes of the data through the KF. 120 UPDATED STATES XZ (X0(2)-75 . Iterations - 3) 400 900 1200 DISCRETE UNIT TIME (sampling interval".01) 1600 UPDATED STATES XI (XO(I)-ZSOO. Iterations - B) 400 BOO 1200 DISCRETE UNIT TIME (sampling interval-.01) 1600 Figure 4 .7 . Updating the States (a) x(2) (Arr iva l Time) and (b) x ( l ) (Amplitude) 121 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 1200 DISCRETE UNIT TIUE (sampling xntmal-01) ISOO Figure 4 . 8 . Mean Square Error for (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, is 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_. The KF f i l t e r i s appl icable for systems that can be described by a f i r s t order d i f f e r e n t i a l equation in x. and a l inear (matrix) equation in z. The state and measurement equations assumed to be in the fol lowing 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-vector , w i s a p-vector , z and v_ are m-vectors, and F, G and H are known matr ices. 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) T) = Q( t )6( t - r ) (4.3b) E(v ( t )v (x ) T ) = R ( t ) « ( t - r ) (4.3c) It is 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 funct ion . Note: The KF can be modified to treat the case when equations (4.3c) and (4.3d) are not t rue ; however, the notation is s impl i f ied considerably when these assumptions are made. 123 It i s further assumed that the i n i t i a l s ta te , x(t )HX^, i s a zero mean, Gaussian random process and that the covariance of i s speci f ied as a pos i t ive semi -def in i te nxn matrix Po. Thus we have E(x(t Q))=0 andE (x ( t o ) x ( t Q ) T ) = P q (4.4) 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 interval t <x<t. It i s convenient — O — Q— — to make the fol lowing d e f i n i t i o n s : D e f i n i t i o n : 1. If t '> t , x_(t'|t) is a predicted estimate. 2. If t '=t , xjt'|t) is 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 is defined as x ( t '| t ) = x ( t ' ) - x ( t ' | t ) . (4.5a) and the performance index of the estimation is given by I (x ( t '| t ) ) = E(e(x(t '|t ) ) ) (4.5b) where E ( . ) i s an error function of the estimation error . Although more general error functions can be assumed, we shal l assume the e is the square of 7 ( t ' | t ) ; i . e . , I (x ( t '|t ) ) = E ( x U ' | t ) x ( t ' | t ) T ) (4.5c) The estimation problem can now be stated as fo l lows: 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 interval t ^ r < t , determine an estimate of £ ( t ' | t ) such that I (7( t '| t ) ) is minimized. The solut ion to th is problem was given by R. E. Kalman (Reference 11) in 1960. This solut ion i s speci f ied by the fol lowing 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 ) T R(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) Note: The equations for the predicted and smoothed estimates are not given at th is t ime; however, i t appears that smoothed estimates should be used to obtain sat i s fac to ry resul ts the our par t i cu la r problem. A block diagram of the system, measurement, and f i l t e r e d estimates is shown in Figure 4.9.. with P ( t Q ) = P Q , given 125 state noi se input matrix W(t) O S ( t ) x(t) x(t) state matrix F(t) measurement matrix > H(t) state vector x(t) i n i t i a l state measurement noi se v(t) measurement z(t) Figure 4 .9a. Block Diagram of State and Measurement Equations (4.1 ) and (4.2 ). The Sol id Lines Indicate Vector Quantities measurements •2(t) Kalman Gain x(t) state estimate x(t) F(t) P(t) covariance calcu lat ions a p n o r i estimate x —o measurement estimates z(t) a p n o r i covariance P Figure 4.9b. Block Diagram of Continuous Kalman F i l t e r , Equations (4.6 ) and (4.7 ) 126 C. Description of Discrete Version of the Kalman F i l t e r In most cases i t is more pract ica l to formulate the KF in discrete form; th is is espec ia l l y true for computational purposes. Sometimes the system and measurement equations are natura l l y in discrete 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 iscrete state and measurement equa- t ions, can be writ ten as x k + 1 = •(k+l,k)x k + r (k+l ,k)w k (4.11) z k + 1 - H ( k + l ) x k + 1 + v k + 1 " (4.12) where x^ + 1 i s an n-vector, wk is a p-vector, and z ^ + 1 and v_ +̂̂ are m-vectors. The nxn matrix $(k+l,k) i s the state t rans i t i on matr ix ; th is matrix can be obtained from the continuous case (Equation 4.1) as the solut ion to the homogeneous (zero input) solut ion to x = F x. The nxp matrix r (k+l ,k) i s the input t rans i t i on matr ix . The subscript k refers to the discrete time t = t k , k=0, 1 The random (vector) processes w^ and are assumed to be zero mean, white noise sequences, so that E ^ ) = 0 and E ( v k + 1 ) = 0 (4.13a) E ^ w / ) = Qk * j k (4.13b) E ( v j + 1 v j + 1 ) - R k + 1 * j k (4.13c) It also assumed that wk and v_k are s t a t i s t i c a l l y independent of each other, and that the i n i t i a l state x.(t )=x^ is a zero mean, Gaussian random process and that the covariance of x i s spec i f ied as a pos i t ive semidefinite nxn matr ix . Thus we have E(ljwJ) = 0 (4.13d) E(x.) = 0 and E ^ ) = P Q (4.14) 127 As in the continuous case, we use the notation x.(k|j) to denote an estimate of the state at time t^ , t ^ t , based on measurements 2̂ . on the in te rva l t ^ t ^ t j . As in the continuous case, i t i s convenient to make the fo l lowing d e f i n i t i o n s : f i l t e r e d est imate. smoothed est imate. The est imation er ror i s defined as x(k|j) * X|< - x(k|j) (4.15a) and the performance index of the est imation i s given by I (x(k|j) = E(e(x(k|j))) (4.15b) where E ( - ) i s an er ror funct ion of the est imation e r r o r s . We sha l l choose e to be the square of 7(k | j ) , i . e . , I ( £ ( k | j ) ) = E (x (k| j )x (k| j ) T ) (4.15c) The d i sc re te est imation problem can be stated as f o l l o w s : Estimation Problem: Given the d isc re te 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 in te rva l t 0 . f t^_ f t j , determine an estimate x_(k|j) such that I(x_(k|j)) i s minimized. Est imation Equations for the F i l t e r e d Est imate, k=j State Estimation equation £ k + 1 = •(k+l,k)xk + K k + i d k + r H k + l * ( k + 1 > k ) i k ) with (4.15) 128 Kalman Gain Matrix Kk+1 = p(k+l|k)H(k+l)T(H(k+l)P(k+l|k)H(k+l)T + R(k+l))_ 1 (4.17) Estimation Error Covariance (Matrix Ricatti) equation P(k+l|k) = *(k+l,k)P(k|k)*(k+l,k)T + r(k+l,k)Q r(k+l,k) T (4.18) P(k+l|k+l) = (l-K k + 1H 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 ^+1=_x^+i~^<+i• 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 Q s p e c i f y ^ , P Q and Q Q, and compute *(1,0), H^, and R .̂ • At t j , compute the projected estimate of the covariance matrix P(l|0)=*(l,0)P o*(l,0) T+Q o . 0 Compute the gain matrix K1=P(l| 0 ) H^(H 1 P ( 1 , 0 ) H|+R 1 ) " 1 • Using the measurement at t=t^, the best estimate of the state at t 1 is given by lx = •d.0)x ( J + K ^ - H ^ l . O ) ^ ) t The estimation covariance matrix at t^ is given by 129 P ( l | l ) = P(1|0) - K 1H 1P(1|0) • At t = t 2 , a new measurement i s obtained, and the computational cyc le i s repeated. 130 D. Derivat ion of Transfer Function for Geophone An electromagnatic device whose output i s proportional to the input acce lerat ion of i t s case i s shown schematical ly in Figure 4.10. When the case of the device moves, in the d i rec t ion of i t s sens i t i ve 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 vo l tage, e Q , as a funct ion of t ime. Figure 4 .10. Schematic of an Accelerat ion Measuring Device (Reference 12 ) In Figure 4 . 1 , the fo l lowing parametres are def ined: M - mass of suspended magnet B - f lux densi ty between magnets 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 131 d I R a c(t) diameter of coil length of co i l , ir-d-N 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 The dynamic equations for the device can be derived by summing the voltages induced in the electrical measuring circuit . These two situations are shown schematically in Figure 4.11. - mi I U/(</c<t) + x'l r.T i Kx Eft Jilx Ac 4- (c) forces on M (b) measuring circuit 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 - Mac ^ 2 1 a ^ and, summing the voltages in the circuit gives L dT + R c 1 + e o = "B** eo = V (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+B2£2/R)/M + K x/M = - a c . Taking the Laplace Transform of these equations gives (Ms2+Ds+K)X(s) = -Bt I(s)- MAc(s) (4.22a) (Ls+Rc+Rr)I(s) = -Bi sX(s) (4.22b) E0(s) = Rr I(s) (4.22c) Since our goal is to obtain the transfer function of the output, E Q(s), over the input, A (s), that is T(s) = E0(s)/Ac(s) it is necessary to eliminate X(s) and I(s) in Equations (4.22). This is accomplished by substituting X(s) from (4.22b) into (4.22a). After some algebraic manipulations we obtain the desired result 133 T(s) = E0(s)/Ac(s) = Rr BAMs/P(s) (4.24a) where P(s) = MLs3 + (DL+MR)s2 + (DR+KL+B2*2)s + KR (4.24b) R = R + R (4.24c) c r 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 Fi l ter . It is noticed that the current in the coil circuit (or voltage across the load resister eQ=i*Rr) is used as the output variable. 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 q0(t) = J e 0(t ' )dt' (4.25a) Then eQ(t) = qQ(t) and EQ(s) - sQQ(s) (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 sQo(s) K rs Qn(s) K; KJJT = rur> o r orr = ^T ( 4 - 2 6 ) 134 where s R f B£M 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) = MRs2 + (DR+B2*2)s + KR (4.27a) so that the transfer functions can be written as where Q (s) K r • °, . = _ E (4.27b) c^ s ' s +2C<" s+« 0 0 Kc = BARr/R w 2 = K/M (4.27c) o c = (D+B2^2/R)/2MwQ In Section II, Table 2.1, the geophone specifications are given as u0~2?*28.0 and 5=0.18. 135 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 s , and arrival times, Tp and T s , at the geophone location are 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 it is seen that they are in the same form, with the following relationship or K c a c ( t ) = w o w { t ) w(t) = K c a c ( t )/^ (4.28) 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>Q and ; are constants scaled acceleration w(t) input s 2 + 2;u» s + a)2 o o w(t) = Kc a c(t)/" magnet position y(t) 9- output Figure 4.12. Block Diagram of Second Order Geophone Transfer Function 136 Remember now that PQ(t) is an output of the instrument that is related to the integral of the voltage across a load resistor R r in the coil c i r - cuit. This could be the voltage across some sort of a capacitance circuit it depends on how a particular instrument is designed. One could assume that the measurements (for the KF formulation) are cfU k)+v(t k), where v(t k) is white noise which could be due to computer roundoff, quantization, circuit noise, etc. Actually it might be better to assume that q(t k) 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) = qQ(t) (4.29a) then we can write y(t) + 2 C « 0 y ( t ) y(t) = K. a c(t) (4.29b) Next define j y : (t ) = y(t), y 2(t) = y(t), i ={yl,y2) > and (4.29c) u(t) = a c(t) 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 (first order correlation will be assumed) Xg = geophone output due to background noise (x,-) x7 = *6 Xg = geophone output due to P&S waves x 9 = x 8 NOTE: x, & x7 and x„ & xQ 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. Actually, the P and S waveforms can be more general than "impulses". 138 "x6(k+D" a l l a12 "x6 (k)" V = x5(k) _x 7(k+l)_ _a21 a 22. _x? (k)_ - b 21. (4.30a) where ^ t ^ j " ^ a n c ^ ^ s t h e 9 e o P n o n e input noise. The constants a ĵ b i j a n later. b.- are defined by Equation (4.6f) of Reference 4, and will be discussed The variables x̂ and Xg are governed by the same discrete equation, except that it is clearer to write the outputs as a function of tp and ts< Thus (assuming tp<t$). "Xg (k+1)" "o" _Xg (k+l)_ It _0_ (4.30b) "x8 (k+1)' " b n " _xg (k+l)_ - b 12- x1 (k+l)/A (4.30c) NOTE: X£(k+1) = tp and x^(k+l) = ap. We divide by A to approximate a 5-function for the P-wave. Generalizations are discussed later. "Xg (k+1)" a l l a12 ' x g ( k ) ' _Xg (k+l)_ _a21 a22- (4.30d) I f V l s t s 139 "x8 (k+1)" a l l a12 "x8(k)" " b u " — + x,(k)'/4 _xg (k+l)_ - a21 a 22. Xg(k)_ . b 2 1 . 0 NOTE: x 4 (k+l ) = t s and x3 (k+1) = a s Thus the "nonlinearity" occurs in the times t and t $ ( i . e . , x2(k+l) and x4(k+l)). 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 rst order correlated, i . e . , i t is assumed that x 5 + Y X 5 = Y W a (4.31a) where Y=l/T , and T c is the time constant of the input noise to the instru- ment, and w is a white noise process with 3 Exp(wa) = 0 and Cov(wa) = q a a a (4.31b) This noise model is discussed in more detail in a separate section. If w3 is assumed constant over a sampling interval, A , then a discrete model a for the input noise can be written as x5(k+l) = awx5(k) + bw wa(k) (4,31c) where a w = e " Y A and b„ = ( l -e" T U ) - Y A , (4.31d) 140 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 1 ,x 2 , . . . , x 9 ) T (4.32a) One can write the discrate state recurrence equation as - I (JVV + (4.32b) where \ + l = (0,0,0,0,b w ,0,0,0)w a(k+l) (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) ( i . e . , t and t $) to see how it 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)my fav is a scalar) (4.33b) cov (v(k+l))= rra and exp (v(k+l)) = 0 (4.33c) i . e . , we have a scalar measurement mv(x^+Xg) with uncorrelated measurement noise having variance rm- 141 We wr i te 1(*K>V = F * k (4.34a) where where F = I = c B s o A ; 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 (4.34b) an i d e n t i t y matrix to represent the 4 constant state var iables b l l a l l a 12 L b 2 1 a 2 1 q 2 2 . the a . , and b , . spec i f ied above 11 '21 0 b i i o" a l l a 1 2 ~ and A' = 0 b 2 1 0 s _ a 21 a 2 2 _ It 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 2 & x, at t=0, and c a l l them x 2 (0) & x 4 ( 0 ) . The KF w i l l not change these estimates u n t i l t k + 1 i x 2 ( 0 ) , and x . ( k + l ) = x . (0) (i=2,4) for t k + 1 < x 2 ( 0 ) . Note: One can see the d i s t i n c t i o n between generating the geophone data with the models, and estimating the parameters with the models. In one case, x ^ O ) , i=l to 4 are known, and in the other case, the estimates x\(0) 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 Fin terms of matrices I . A , B , and c n s A because when you propagate ahead in the KF you need to make s T operations likeFPF , and, by partitioning F, one can make simpli- fications over the intervals ( O ^ ) , ^ J X ^ ) , a n d (X^,T), where T is the end time of the geophone trace. 143 The Kalman Filter Equations for Velocity Measurements We have the state equation, at time t k + 1 =(k+l) , • KvV + \ n ( 4 - 3 5 a ) where w | ( + 1 = (0,0,0,0,b w ,0,0,0,0) Tw a(k+l) cov(wa(k+l)) = q f l and the measurement equation z k + l " H i i k + 1 + v m^ k + 1 ) c o v ( v m ( k + 1 ) = rm ( 4 ' 3 5 b ) where H = (0 0 0 0 0 0 my 0 m ) 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 init ial estimate of the state and estimation matrix, x̂ and P , are available. Step 1 Compute $>=3f/3_x i f x 2(k+l) <_ t k + 1 <_ x 4 (k+l) + N pA, or if $4(k+l) < t L + 1 <_x4(k+l) + N sA. Otherwise, *=F. Step 2 Propagate the estimate of the state ahead one step, i . e . , » V i • I<VV < 4 - 3 6 a ) Step 3 Propagate Covariance matrix ahead P'k+1 = *Pk * T + Q (4.36b) 2 where all the elements of Q are zeros except Q(5,5) = b q W Q 144 Step 4 Compute the updated measurement noise covariance - » » ( P 7 7 + P 7 9 t P 9 7 t P 9 9 ) * rm a s c a 1 a r The p ĵ refer to elements of matrix P'̂ +j Step 5 Compute the Kalman Gain matrix K, T P ' . •. H T (4.36c) k+l • P i ? + P i 9 " »26 + P 2 9 L P 9 7 + P 9 9 j Step 6 Compute measurement residuals A z k + l = z k + l " H i ' k + l = z k + r m v ( x ' 7 ( k + l ) + x g ' ( k + l ) ) Step 7 Compute new State Estimate *k+l = x ' k + l + K k + l A z k + l Step 8 Update State Estimation Error Covariance matrix P k + l = P ' k + l " K k + l H P ' k + l (4.36d) (4.36e) (4.36f) (4.36g) 145 N o t e : S i n c e H P ' K + 1 = j n v ( p i 7 + P I 9 ' - - - ' P Q 7 + P Q Q ) ' a n d 97 H 9 9 ' r K l l ( P l 7 + p l 9 } K n ( p Q 7 + PQQ^ K k + 1 H P k + 1 1 1 ^ 9 7 H 9 9 J - K 9 1 ( p 1 7 + P 1 9 } K 9 1 ( P 9 7 + P 9 9 ) A c o n s i d e r a b l e n u m b e r o f t h e c o m p u t a t i o n s a r e s a v e d by u s i n g t h e f a c t t h a t H P j ^ + 1 i s a 1 x 9 v e c t o r a n d K ^ + 1 i s a 9 x 1 v e c t o r . 146 F, Comments on Computing the Estimation Error Covariance Mat r i x , P Since a considerable port ion of the c a l c u l a t i o n s required to obtain the KF estimates are re lated to updating the Estimation Error Covariance Mat 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 read i l y implemented. The essent ia l computational e f f i c i e n c i e s are r e - lated to the structure of the state t r a n s i t i o n matr i x , <t>, and to the rather simple structure of the measurement mat r i x , H. That i s , the upper l e f t 4x4 port ion of the * matrix i s an i d e n t i t y mat r i x , the upper r ight 4x5 matrix i s a nul l matr ix , 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 fo l lowing equation: P' Kk+1 * P k *' + Q (4.37a) where, as prev iously noted, a l l elements of the state noise covariance matrix Q are zero , except for the Q(5,5) element. The state t r a n s i t i o n matrix « i s of the fo l lowing form: * = *44 °45 B 54 A 5 5 (4.37b) where denotes a 4x4 i d e n t i t y matrix O^g denotes a 4x5 nul l matrix 55 a w 0 0 0 0 b l l a l l a 12 0 0 b 2 1 a 21 a 22 0 0 0 0 0 a l l a 12 0 0 0 a 2 1 a 22 and 147 '54 0 0 0 0 0 0 0 0 0 0 0 0 11 c l l b h c h 21 c 2 1 b 2 1 C 2 1 The elements of the and B^4 matrix have been prev iously def ined. It i s noted that i f the time t k + 1 f a l l s into any of the fo l lowing in te rva l s t k + 1 < x 2 ( k + l ) or or x 2 (k+l) + N p .A < t k + 1 <x 4 (k+l) (4.38) t k + 1 > x 4 ( k + l ) + N s -A then Bg 4 = 0 ^ 4 , i . e . , a 5x4 nul l matr ix , The updated P matrix for time t k + ^ i s given by the fo l lowing equation: P k + 1 •• P M " K M H P M ^ 3 9 > where K k + ^ i s the Kalman Gain matrix given by Kk+1 = Pk+1 " W + l (4.40) where r k + l = H P k + 1 + rm (4.41) 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 = mv(0 0 0 0 0 0 1 0 1), for the velocity measurement case, it is evident that P'H T = m. Pi7 + Pi7 L.P97 + p99 (4.42a) and r' = ^ ( P ) 7 + P ) 9 + P 97 +P99) + rm (4.42b) where the p'. • are elements of P'. The gain matrix has the fair ly simple form 1 T m i r P ' H 1 =A r r Pi7 + Pig' P97 + P99 (4.43) It is next necessary to partition the P and P' matrices into the same structure as the partitioning of the * matrix. Thus we define P ^ a t L b c -1 (4.44a) where 149 P , denotes a 4x4 symmetric matrix a P. denotes a 4x5 matrix b P c denotes a 5x5 symmetric matrix With these definitions i t is seen that the following matrix product is obtained * P * = P A P B T + P . A T a a b BP + A P J BP B T +BP . A T + A PIB T + A P A T _ a b c b b c - i (4.44b) For convenience of notation, the following matrix products are defined: X = B P K A ' b Y = P B T + P. A T a b (4.44c) and Z = B P ^ 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 fal ls into one of the intervals defined above, then the computations are reduced considerably since X = 0 Y = P. A b (4.44d) 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 it is seen that •YT Z+X+XT+APCAT- + Q (4.45) In the final step of updating the P matrix, it is convenient make a different partitioning of P'. This partitioning is based on the structure of the H matrix. In particular, it is readily seen that P = P' P" (4.46a) where p« = i , . P' (HTH)P' (4.46b) Since T 2 H ' H = mc v ro 1 o 2 L ° 2 '3 151 where denotes a 6x6 nu l l matrix 0 2 denotes a 6x3 nu 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' in to the same s t r u c t u r e , i . e . , P' P i P 2 P l P 3 (4.46c) where P^, P 2 and P^ are dimensioned the same as 0 ^ , 0 2 , and Î , r e s p e c t i v e l y . With t h i s p a r t i t i o n i n g , i t i s seen that H THP' = °1 °2 L P i T P 3 J (4.46d) and P ' ( H ' H ) P 1 P 2 P 2 T P 2 P 3 L P ^ P ^ (4.46e) where, i t i s noted that Pĵ and P^ are symmetric matr ices. 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 Equations" Section VI. 152 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 . , f ( t . . ) d t - y | i ( t . . ) d t + f U l , . , _ 1 . f ( + 2 > c , ) _ i (4.47) In the present case, ^ ( a ^ x ^ n A - A and <j>2(a)=x2+nA, n=l to N . Thus iand_ j_ i = i (4.48a) dx 2 dx 2 Also f ( t,a) i s not a funct ion of a ( i . e . , x 2 ) thus 3 f ( t > o t ) = Q (4.48b) 3a Thus 4 - / fir, J x2+nA -r- , f ( t ) d t = f ( x „ + n A - A ) - f ( x 0 + n A ) (4.48c) d a x 2+riA-A 2 1 Since the present problem is time invar iant ( i . e . , the so lut ion for the second order system depends only on the d i f ference x 2 + A - x 2 = A , we need only consider the l i m i t s 0 and A . NOTE: One needs to be careful on what the values of a (t) (and a ( t ) ) are over t h i s i n t e r v a l . That i s , the l i m i t s 0 and A d o not change, but the values of ap(t) do. Thus the l i m i t s of in tegrat ion are <j>2(a) = x 2 + n A (4.48d) ^ ( a ) = x 2 + ( n - l ) A n = l , 2 , . . . , N p . 153 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 arguments are also true for the S-wave). ' P * W = W ' w h e r e t h e w p a r e predefined scalars x2+NpA T—l »—f-ff^ r- U time, t. Ap-wp(k) According to Equation (4.6b) of Reference 4, we have A F - (A - t ) e 01 2 ~C A-A _<-' t r - s in p ' (A - t ) s inp' (A- t )+cosp*(A-t ) dt -a p ( t k ) Define f(t) = e r I - s inp ' (A-t ) 01 - c s in p'(A-t) + pcosp' (A-t ) ap(V (4.49a) then f ( 0 ) - r - sin p 0) .-c sin p"+pcosp"_ y v (4.49b) and f(A) = e 0 P a P < W (4.49c) 154 Thus f ( 0 ) - f ( A ) = - sinp" _COS(p"+a)-pe" A -(w (t. )-w (t. A .)) p v pv k' pVk+1'' (4.50a) where p" =p'A = p-cu-A C" =5 A = ?-oj-A (4.50b) k = x2/A to x2/A + Np Note: x2 and x4 should be selected/estimated so that they are multiples of the sampling time A. Therefore _d_ dx. f f(t)dt aie sinp" • ap(t") cos(P"+a) • ap(t')-paj 2e" 2 ?"ap(t +) a) e - sin P " 01 ap(t ) -r" + cos(p"+a) • ap(t ) - p e ^ ap(t ) (4.51a) where t = x2 + n A - A t = x2 + n A (4.51b) n = 1 to Np Alternatively, if we use Equation (D.6.e) of Appendix D we have 155 / .2 r " 4 cosp" I sinp" sin(p"+a) COS(p"+a) i icosp't ii - c ' t I sinp't dt-w. . I 4 . I K where a = sin c Applying similar notation, i .e . , letting f ' ( t ) = e sinp't COSp't (4.52a) we obtain f'(0) = e " 0 ( J ) -p(O (4.52b) f ( A ) = e" c ,. /sinp"\ \cosp" / p (4.52c) Thus 0 aj2e i sinpM .a p(t')-e' c"(-cosp"sinp"+sinp"cosp")a p(t +) cos(p"+a).a n(t")-e" ;"(sinP+a)sinp"+cos(p"+a)cosp") a. OJ e i s i n p » . a p ( t - ) cos(p"+a)-a p (t ' ) - P e' ? "a p (t + ) (4.53a) which is the same result as above. Since this result holds for both the P-wave and the S-wave, we can write, for i=2 and 4, 156 A. J L f J 2 c f(t)dt = S_i x^+nA-A sin P".w. , ( t ) OJO 1-1 cos(P"+a)-w._1(t")-pe"?"w i_1(t+)_ (4.53b) where w^(t") = x^(t")-Wp(t") t" = x 2+n-A-A, n=l to Np w (t~) = x-,(t~)-w (t~) t" = x.+nA-A, n=l to N (4.53b) Note: When the ŵ 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 2 ( t " -A ) , x 3 ( t _ - A ) , and x 4 ( t ' -A ) (4.53c) are used. If the state vector is defined by 3F , . then F is used for updating x^, and $ = is used for updating the estima- tion error covariance matrix P .̂ It is noted that $ = F for t^<x2 and t | c > x ^ + N s - A and *=F for x2<_t^x4+NA except for the following elements: 3F 8 2 /3x 2 , 3F g 2 /3x 2 , 3F 8 4 /3x 4 , and 3F g 4 /3x 2 In summary, these elements can be written for i=2 and 4, as 8i L/9iJ = c. c 2 • w . ^ t ) c 3 • w-.^t") - c 4 • W ^ t ) 157 (4.54a) c 0 = — sin p" p" = pcuA (4.5b) C, = COS(p"+a) a - s i n ' 1 5 X i - 1 ^ 0 = best estimates of at t", i =2 and 4 w x(t _) = wp(x2+nA-A) n = 1 to Np w^t*) = w^t'+A) (4.54c) w3(t") = wp(x4+nA-A) n = 1 to N$ w3(t+) = w3(t_+A) 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 th is 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 ve loc i t ies is 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 po lar i t y of shear waves to determine a r r i va l times of shear waves at s p e c i f i c depth interva ls to calculate shear wave ve loc i t ies over each depth i n t e r v a l . The reverse po la 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 va 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 is 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 va l time from the oppositely polar ized waves. The determination of the cross-over can be d i f f i - cu 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 ass 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 to r t the signal and erroneously a l t e r a r r i v a l times. An addit ional consideration i s that the reverse po lar 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 invest igat ion . But i t has been found that a s ingle c r i t e r i o n may not be appropriate throughout the seismic p ro 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 layer - ing. The two f i l t e r i n g methods considered in th is thesis are the one based on frequency domain f i l t e r i n g , where Fast Four ier , Butterworth, and crosscorrelat ion 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 crosscorre lat ion 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 crosscorrelat ion funct ion. The crosscorre lat ion function is applied to two successive waves (corresponding to the successive depths) where the maximum value of the cross- cor re lat ion function i s d i r e c t l y related to the difference in travel times 159 between the waves. With th is a r r i v a l time d i f ference, a ve loc i ty can be calculated. The method i s very advantageous because i t allows one to focus on a spec 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 rely ing upon dominant responses as opposed to the reverse po lar i t y method where there is dependence on one point ( i . e . , the cross-over) to define ar r i va l time dif ferences. 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 cross- corre lat ion 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 . s h i f t s . The-crosscorrelation f i l t e r was coded in Fortran language, and the program is referred to as CROSSCOR. CROSSCOR is a graphics interact ive program which displays the frequency spectra , unf i l te red and f i l t e r e d times series and cross- correlat ions 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. CROSSCOR was 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 tes are given, and the ve loc i ty p ro f i l es from di f ferent methods of data acquis i t ion 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 Site in Richmond, B.C. In th is case study f i ve sets of data were compiled and two d i f fe rent 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 Pressure- meter, which had two accelerometers i n s t a l l e d . The seismic data was reduced by the reverse po lar i t y method and CROSSCOR, and compared with cone bearing p r o f i l e s . The second s i te 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 ferent types of sources were compared in th is inves t igat ion , i . e , the shear and compression waves sources. The th i rd research s i t e analyzes data acquired at the T i lbury National Gas Plant S i te on Ti lbury Island where the UBC seismic cone pressuremeter was used for data acqu is i t ion . The fourth s i te was the Annacis Bridge P i l e Research S i te where the UBC Seismic Cone Pressuremeter and UBC #8 Seismic Cone were used for data acqu is i t ion . In th i s invest igat ion shear ve loc i t i es were obtained using shear sources. Since both sides of the source were used, i t was possible to obtain two veloci ty estimates for one depth increment using CROSSCOR as opposed to the reverse po lar i t y method which allows for only one ve loc i ty estimate. The las t research s i t e was the Lower 160 was the Lower Langley 232nd St . S i t e . In th is invest igat ion f i l t e r e d (analogue) and unf i l te red reverse polarized Hammer Shear Source data was obtained. In add i t ion , 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 tes indicate that CROSSCOR is a very benef ic ia 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 ferent 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 - t i e s . The procedure in extract ing the desired instrument responses proved to be routine and user f r iend ly due to the d i f fe rent responses being consistent throughout a cone p r o f i l e . In add i t ion , since CROSSCOR i s a d i g i t a l f i l t e r i t removes most human bias ( e . g . , determining cross-overs in the reverse po la 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 es 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 ferent sources in obtaining seismic wave v e l o c i t i e s . The problem with tr igger ing when applying the Buffalo Gun Source may be resolved by using an accelerometer (having a fast response time) as the t r igger as well as grounding the Buffalo Gun f i r i n g pin to create a contact closure (s imi lar 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 appl icat ion to obtaining accurate estimates on the P-wave and S-wave amplitudes and a r r i v a l times. The KF is a very f l e x i - ble tool which allows one to model the problem considered accurately. In addi - t i o n , the KF works in the time domain which removes many of the l imi tat ions of the frequency domain techniques. 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. Appl icat ion of th is 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 di f ference equations which, with i n i t i a l condit ions, 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 character is t ic 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 is than fed into the Kalman F i l t e r with the a r r i va 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 ne l y tuned, i t would be possible to obtain a r r i v a l times and amplitudes on l i ne resul t ing in ve loc i t i es and damping c h a r a c t e r i s t i c s , respect ively . 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 is a very important parameter in earthquake engineering). The damping obtained from the KF would be a de f in i te advantage over conventional techniques which re ly on power spectrum measurements to determine the damping character is t ics 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 races) , there i s greater chance of error than the Kalman F i l t e r formulation. 162 VI. SUGGESTIONS FOR FUTURE RESEARCH The research in th is thesis indicates that d i g i t a l f i l t e r s are extremely benef ic ia l in i n - s i t u tes t ing . As was out l ined , there are many d i f ferent types of f i l t e r s one can implement; therefore, i t i s important to become f luent in a l l the d i f ferent f i l t e r i n g techniques and apply the ones which best su 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 analys is . Future research at UBC should concentrate having the i n - s i t u group becoming more f luent 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 tes t ing . One of these parametres is damping of the s o i l . So i l damping and the shape of the seismic wavelets are very important parametres in earthquake engineering. Another area of interest 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 . , velo- c i t y and damping) that the waves have t rave l led through. Therefore, i f one establishes accurate methods for determining P and S-wave ve loc i t i es 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 into 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 for seismic data reduction. 163 LIST OF REFERENCES 1. Richart , F. E . , J r . , H a l l , J . R., J r . , and Woods, R. D., (1970), "Vibrations of So i l s and Foundations," P r e n t i c e - H a l l , Inc . , Englewood C l i f f s , New Jersey, 414 pp. 2. Mooney, H. M. (1974), "Seismic Shear Waves in Engineering," Journal of the Geotechnical Engineering D iv i s ion , ASCE, Vol . 100, No. GT8, Aug. Proc, Paper 10745, pp. 905-923. 3. Borm, G. W., (1977), "Methods from Exploration Seismology: Ref lec t ion , Refract ion, and Borehole Prospecting/' Proceedings of DMSR 77, Karlsruhe, 5-16 Sept. 1977, Vo l . 3 , pp. 87-114. 4. Rice, 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 Test ," 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. 6. Close & Fredr ick, Modelling and Analysis Dynamic Systems, Houghton Mif1 in Co. , 1978. 7. Ogata, J . , Discrete Process Control , Cambridge University Press, 1987. 8. Ewing, Jardetz , 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 Affect ing In -Si tu Seismic Measure- ments," Proceeding of Conference on Earthquake Engineering and Soi l Dynamics, ASCE Goetechnical Engineering Div is ion 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. Sher i f f and L. P. Geldart , Exploration Seismology, Volume 2, Data Processing and Interpretat ion, Cambridge University Press, 1983. 12. George W. St inson, Introduction to Airborne Radar, Hughes A i r c r a f t Co. , 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 D ig i ta l F i l t e r s : with Appl icat ion to Data Processing on the UNIVAC 1106 Computer," P a c i f i c Marine Science Report 80 -9 , Inst i tute of Ocean Sciences, Sidney, B.C. 15. E. A. Robinson, "Co l lect ion of Fortran II Programs for F i l t e r i n g and Spectral Analysis of Signal Channel Time Ser ies , " Chapter 22 of textbook by Unknown. 16. U.B.C. P lo t t ing Programs, Unpublished Notes. 17. Campanella, R. G . , Robertson, P. K., and G i l l e s p i e s , D., 1982, "Cone Pene- t ra t ion Testing in Delta ic S o i l s , " Canadian Geotechnical Journal , Vo l . 20, February. 18. Howie, J . A . , 1989, Ph.D. Thesis ( in preparat ion) , University of B r i t i s h Columbia, Vancouver, B.C. 19. Davies, M. , 1987, "Predict ing A x i a l l y and Latera l ly and P i l e Behaviour Using In -S i tu Testing Methods," M.A.Sc. Thesis, Department of C i v i l Engineer- i n g , University of B r i t i s h Columbia, Vancouver, B.C. 20. Greig , 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," Ph.D. , University of Bath-England, 1986. 23. Asten, M. W., "Theory and Pract ice of Geophone Cal ibrat ion Ih -S i tu Using Modified Step Method," IEEE Trans. Geosc., Vo 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 Appl ication to Seismic Time Series Ana lys is , " 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 Predict ion Problems," J . Basic Energy, Series D, Vo l . 45, pp. 25-36, March 1960. 165 APPENDIX A Program l i s t i n g of CROSSCOR 166 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27, 1988 f o r C C i d = S E I S on G 1 REAL CI 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 ) 2 REAL X X ( 2 0 4 8 ) , Y Y ( 2 0 4 8 ) , F C 1 , F C 2 , F C 3 , F C 4 . D 1 ( 8 ) , G 1 3 REAL A V R G X , A V R G Y , X 2 ( 2 0 4 8 ) , Y Z ( 2 0 4 8 ) , D 2 ( 8 ) , G 2 , D E L T 4 REAL D L T , G ( 5 0 0 0 ) , G A ( 1 0 0 0 0 ) , C M A X , E X ( 1 0 0 0 0 ) 5 INTEGER LG,LA,LMAX,HMAX 6 REAL X 1 ( 2 0 4 8 ) ,Y 1 ( 2 0 4 8 ) ,FN 7 c 8 REA D ( 3 , * ) N1,DELT 9 DLT=1000.0*DELT 10 FN=0.5/DELT 1 1 DO 100 1=1,N1 12 R E A D ) 1 , * ) X ( I ) 13 R E A D ( 2 , * ) Y ( I ) 14 100 CONTINUE 15 c 16 c 17 CALL NORME(N1,X) 18 CALL NORME(N1,Y) 19 CALL REMAV(N1,X,AVRGX) 20 CALL REMAV(N1,Y,AVRGY) 21 CALL P0W2IN1.X.Y.N) 22 c 23 DO 102 1=1.N1 24 X X ( I ) = X ( I ) 25 Y Y ( I ) = Y ( I ) 26 102 CONTINUE 27 c 28 CALL F A S T F ( X , Y , N ) 29 CALL F X F F T ( X 1 , Y 1 , X , Y , N ) 30 CALL SMOOTH(XT.N) 31 CALL SMOOTH(Y1,N) 32 CALL RAD(W,N,DELT) 33 CALL NORME(N,X1) 34 CALL NORME(N,Y1 ) 35 NN=N/2 36 DO 101 1=1,NN 37 WRITE! 7 , * ) W ( I ) , X 1 ( I ) 38 WRITE18,*) W ( I ) .Y 1 ( I ) 39 101 CONTINUE 40 c 4 1 CALL PL(X1,W,NN) 42 . c 43 W R I T E ( 6 , 4 ) 44 4 F O R M A K ' W h a t f r e q u e n c i e s do y o u want f i l t e r e d ? ' ) 45 R E A D ( 5 , * ) F C 1 , F C 2 46 c 47 c 48 CALL PL(Y1,W,NN) 49 c 50 W R I T E ( 6 , 4 ) 51 R E A D ( 5 , " ) F C 3 , F C 4 52 c 53 c 54 CALL BNDP'S(FC1,FC2,DLT,D1 ,G1 ) 55 CALL F I L !. ,(XX,N1,D1 ,G1 , IG1 ) 56 CALL B N D P A S ( F C 3 , F C 4 , D L T , D 2 , G 2 ) 57 CALL F I L T E R ( Y Y , N 1 , D 2 , G 2 , IG2 ) . 58 CALL NORME(N1,XX) 167 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27. 1988 f o r C C i d = S E I S on G 59 CALL N0RME(N1,YY) 60 c 61 DO 201 I=1,N1 62 W R I T E O , * ) X X ( I ) 63 W R I T E M C , ' ) Y Y ( I ) 64 201 CONTINUE 65 C 66 LG=2*N1-1 67 c 68 CALL CROSS(N1,XX,N1,YY,LG,G) 69 C 70 LA = 2 * LG 71 DO 30 J=1.LG 72 I J = J - 1 73 G A ( J + L G ) = G ( J ) 74 30 CONTINUE 75 C 76 CALL C R0SS(N1,YY,N1.XX.LG.G) 77 LH=LG+1 78 DO 40 11=1,LG 79 G A ( I I ) = G ( L H - I I ) 80 40 CONTINUE 81 C 82 c 83 DO 50 1=1.LG 84 E X ( I ) = I - L G 85 50 CONTINUE 86 c 87 LGG=LG+ 1 88 c 89 DO 53 I=LGG,LA 90 E X ( I ) = I - L G 91 53 CONTINUE 92 c 93 DO 109 1=1,LA 94 WRITE( 1 1 ,*) E X ( I ) , G A ( I ) 95 109 CONTINUE 96 c 97 CALL MAXSN(LA,GA,CMAX,LMAX) 98 c 99 HMAX=LMAX-LG- 1 100 c 101 CALL P L A ( G A , E X , L A ) 102 c 103 W R I T E ( 6 , 3 0 0 ) HMAX 104 300 FORMAT ('The Maximum Time S h i f t i s : ' ,17) 105 c 106 W R I T E ( 6 , 3 0 3 ) FN 107 303 FORMAK'CORRECT I F NOT F I L T E R E D PAST NYQUIST Fn = 108 c 109 STOP 1 10 END 111 c 112 SUBROUTINE F A S T F ( F R . F I , N ) 1 1 3 c 114 c N i s t h e number o f d a t a p o i n t s = 2 1 15 c FR i s t h e r e a l d a t a s e t 1 16 c F I i s t h e i m a g i n a r y p a r t o f d a t a s e t ( = 0 . 0 i f o n l y 168 L i s t i n g 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 117 C F i r s t c o m p u t e M 118 C 119 REAL F R ( N ) , F I ( N ) , G R , G I , E R , E I . E U . E Z 120 M=0 121 KD=N 122 1 KD=KD/2 123 M=M+1 124 I F ( K D . G E . 2 ) GO TO 1 125 ND2 = N/2 126 NM1 = N-1 127 L=1 128 C 129 C S h u f f l e i n p u t d a t a i n b i n a r y d i g i t r e v e r s e o r d e r 130 C 131 DO 4 K=1,NM1 132 I F ( K . G E . L ) GO TO 2 133 GR = F R ( L ) 134 GI = F I ( L ) 135 F R ( L ) = F R ( K ) 136 F K L ) = F K K ) 137 F R ( K ) = GR 138 F K K ) = GI 139 2 NND2 = ND2 140 3 I F (NND2.GE.L) GO TO 4 141 L = L-NND2 142 NND2 = NND2/2 143 GO TO 3 144 4 L = L + NND2 145 P I = 3 . 14 159265 146 C 147 C F i r s t a r r a n g e a c c o u n t i n g o f M s t a g e 148 c 149 DO 6 J= 1 , M 150 N J = 2"* J 151 NJD2 = N J / 2 152 EU = 1 .0 153 EZ = 0.0 154 ER = C O S ( - P I / N J D 2 ) 155 E I = SINK - P I / N J D 2 ) 156 C 157 C Compute F o u r i e r t r a n s f o r m i n e a c h M s t a g e 158 C 159 DO 6 IT=1,NJD2 160 DO 5 IW=IT.N.NJ 161 IWJ = IW + NJD2 162 GR = F R ( I W J ) * E U - F I ( I W J ) " E Z 163 GI = F I ( I W J ) * E U + F R ( I W J ) * E Z 164 F R ( I W J ) = F R ( I W ) - GR 165 F I ( I W J ) = F I ( I W ) - GI 166 F R ( I W ) = F R ( I W ) + GR 167 5 F I ( I W ) = F K I W ) + GI 168 SEU = EU 169 EU = S E U ' E R - E Z ' E I 170 6 EZ = EZ * ER + S E U ' E I 17 1 RETURN 172 END 173 C 174 C F X F F T c o r r e c t s t h e a m p l i t u d e s o f F A S T F , a s o u t l i n e d by K a n a s e w i c h 169 L i s t i n g 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 175 C 176 SUBROUTINE F X F F T ( X 1 , Y 1 , X , Y , N ) 177 REAL X ( N ) , Y ( N ) , X 1 ( N ) , Y 1 ( N ) , A ( 3 0 0 0 ) , B ( 3 0 0 0 ) , C ( 3 0 0 0 ) 178 REAL D ( 3 0 0 0 ) , E ( 3 0 0 0 ) , F ( 3 0 0 0 ) , G ( 3 0 0 0 ) , H ( 3 0 0 0 ) 179 DO 10 1 = 1 ,N 180 A ( I ) = . 5 * ( X ( I ) + X ( N + 2 - I ) ) 181 B ( I ) = . 5 * ( Y ( I ) - Y ( N + 2 - I ) ) 182 C ( I ) = A ( I ) ' A ( I ) 183 D ( I ) = B ( I ) *B( I ) 184 X 1 ( I ) = S Q R T ( C ( I ) + D ( I ) ) 185 E ( I ) = . 5 * ( X ( I ) - X ( N + 2 - I ) ) 186 F ( I ) = . 5 * ( Y ( I ) + Y ( N + 2 - I ) ) 187 G ( I ) = E ( I ) * E ( I ) 188 H ( I ) = F ( I ) " F ( I ) 189 Y 1 ( I ) = S Q R T ( G ( I ) + H ( I ) ) 190 10 CONTINUE 191 RETURN 192 END 193 C 194 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 ) , 195 C 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 . 196 C 197 SUBROUTINE POW2(L,X,Y,N) 198 DIMENSION X ( L ) , Y ( L ) 199 REAL TT,ST,SS,S 200 S S = 0 . 3 0 1 0 3 201 T T = F L 0 A T ( L ) 202 S T = L 0 G 1 0 ( T T ) 203 S=ST/SS+1.0 204 Q = I F I X ( S ) 205 N=1 206 DO 20 1 = 1 ,Q 207 N=2'N 208 20 CONTINUE 209 DO 10 I=L,N 210 X ( I ) = 0 . 0 211 Y ( I ) = 0 . 0 212 10 CONTINUE 213 RETURN 214 END 215 C 216 C NORME i s n o r m i 1 i z a t i o n o f an a r r a y by i t ' s RMS e n e r g y . 217 C 218 SUBROUTINE NORME(LL,F) 219 REAL F ( 4 0 0 0 ) , F S , 2 220 INTEGER FF 221 . C I F ( L L ) 3 0 , 3 0 , 1 0 222 CALL M A X S N t L L , F , F S , F F ) 223 W R I T E ( 6 , * ) FS 224 Z=FS 225 DO 20 1=1,LL 226 20 F ( I ) = F ( I ) / Z 227 RETURN 228 END 229 C 230 C DOT I S DOT PRODUCT 231 SUBROUTINE DOT ( L , X , Y , A N S ) 232 I M P L I C I T REAL'8 (A-H.O-Z) 170 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27, 1988 f o r C C i d = S E I S on G 233 REAL * 8 X ( 2 ) , Y ( 2 ) ,ANS 234 ANS=0.0 235 I F ( L ) 3 0 , 3 0 , 1 0 236 10 DO 20 1 = 1 ,L 237 20 ANS = A N S i-X(I)'Y(I) 238 30 RETURN 239 END 240 C 241 C REMAV r e m o v e s t h e a r i t h m e t i c a v e r a g e i n o r d e r DC s h i f t s a r e r e m o v e d 242 C 243 SUBROUTINE REMAV(LY,Y,AVERG) 244 DIMENSION Y ( 2 ) 245 S=0. 246 DO 10 1=1,LY 247 10 S=S+Y(I) 248 AVERG= S/FLOAT(LY) 249 DO 20 1 = 1 ,LY 250 20 Y ( I ) = Y ( I ) - A V E R G 251 RETURN 252 END 253 C 254 SUBROUTINE RAD(W,N,DELT) 255 DIMENSION W(N) 256 REAL SS.ST 257 c 258 DO 10 1=1,N 259 S S = F L 0 A T ( I ) 260 ST=FLOAT(N) 26 1 W (I ) = S S/(ST'DELT) 262 10 CONTINUE 263 RETURN 264 END 265 c 266 c S m o o t h c o m p u t e s s m o o t h e d s p e c t r u m f r o m t h e c o s i n e t r a n s f o r m 267 c u s i n g T u r k e y + H a m m i n g f o r m u l a . The s u b r o u t i n e i n p u t s a r e 268 c S P E C T = u n s m o o t h e d s p e c t r u m 269 c L S = i t s l e n g h t 270 c The s u b r o u t i n e o u t p u t i s 27 1 c S P E C T = s m o o t h e d s p e c t r u m 272 c 273 SUBROUTINE SMOOTH(SPECT,LS) 274 c 275 DIMENSION S P E C T ( L S ) 276 MM=LS- 1 277 A = . 5 4 * S P E C T ( 1 ) + . 4 6 * S P E C T ( 2 ) 278 B = . 5 4 * S P E C T ( L S ) + . 4 6 * S P E C T ( M M ) 279 S J = S P E C T ( 1 ) 280 S K=SPECT(2) 281 DO 10 J=2,MM 282 S I = S J 283 SJ=SK 284 S K = S P E C T ( J + 1 ) 285 S P E C T ( J ) = . 5 4 * S J + . 2 3 ' ( S I + S K ) 286 10 CONTINUE 287 SPECT(1 ) = A 288 S P E C T ( L S ) = B 289 RETURN 290 END 171 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27, 1988 f o r C C i d = S E I S on G 29 1 C 292 SUBROUTINE P L ( A M P , Z , L ) 293 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 ) 294 INTEGER L.FSMPL 295 1 N=C 296 DO 100 I-1 . L 297 s s m = z m 298 T T ( I ) = A M P ( I ) 299 I F ( I . E Q . L ) GO TO 200 300 100 N=I 301 C 302 C 303 C I F EQUAL TO ZERO TERMINATE PROGRAM 304 C 305 200 I F ( N . E O . O ) GO TO 300 306 C CALL SUBROUTINE TO DO THE PLOTTING 307 CALL SUMFUN(SS,TT,L) 308 C AND RETURN FOR MORE DATA 309 C PROGRAM COMES HERE WHEN F I N I S H E D . 310 300 CONTINUE 31 1 C TERMINATING PLOTTING THEN STOP 312 CALL PLOTND 313 RETURN 314 END 315 C 316 SUBROUTINE SUMFUN(X,Y,N) 317 C SUBROUTINE DRAWS A L I N E THROUGH N POINTS 318 REAL X ( N ) , Y ( N ) 3 19 C REAL A M P S { 3 0 0 0 ) , Z S ( 3 0 0 0 ) 320 321 C 322 C AMPMIN=AMP(1) 323 324 C DO 10 1=1,N 325 C A M P S ( I ) = S N G L ( A M P ( I ) ) 326 C Z S ( I ) = S N G L ( 2 ( 1 ) ) 327 C AMPMIN=AMIN1(AMPMIN,AMPS(I)) 328 C 10 CONTINUE 329 c F I R S T SCALE THE POINTS 330 CALL S C A L E t X , N , 1 0 . 0 , X M I N . D X , 1 ) 331 CALL SCALE!Y,N,10.O.YMIN.DY,1) 332 CALL A X I S t O . , 0 . , ' F R E Q U E N C Y ' , - 9 , 1 0 . 0,0. .XMIN.DX) 333 CALL A X I S I O . ,0. . 'AMPLITUDE' ,9, 10.0 ,90. ,YMIN,DY) 334 C F I N A L L C PLOT THE L I N E S 335 CALL L I N E U ,Y ,N,1) 336 c MOVE THE O R I G I N AND RETURN 337 RETURN 338 END 339 c 340 SUBROUTINE P L A ( A M P , Z , L ) 34 1 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 ) 342 INTEGER L,FSMPL 343 1 N=0 344 DO 100 1=1,L 345 S S ( I ) = Z ( I ) 346 T T ( I ) = A M P ( I ) 347 I F ( I .EQ.L) GO TO 200 348 100 N=I 172 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27, 1988 f o r C C i d = S E I S on G 349 C 350 C 351 C I F EQUAL TO ZERO TERMINATE PROGRAM 352 C 353 200 I F ( N . E Q . O ) GO TO 300 354 C CALL SUBROUTINE TO DO THE PLOTTING 355 CALL S U M F N 1 ( S S , T T , L ) 356 C AND RETURN FOR MORE DATA 357 C PROGRAM COMES HERE WHEN F I N I S H E D . 358 300 CONTINUE 359 C TERMINATING PLOTTING THEN STOP 360 CALL PLOTND 361 RETURN 362 END 363 C 364 SUBROUTINE SUMFN1(X.Y.N) 365 C SUBROUTINE DRAWS A L I N E THROUGH N POINTS 366 REAL X ( N ) , Y ( N ) 367 C REAL A M P S ( 3 0 0 0 ) , Z S ( 3 0 0 0 ) 368 369 C 370 C AMPMIN=AMP(1) 371 372 c DO 10 1=1,N 373 C A M P S ( I ) = S N G L ( A M P ( I ) ) 374 C Z S ( I ) = S N G L ( Z ( D ) 375 C AMPMIN=AMIN1(AMPMIN,AMPS(I)) 376 C 10 CONTINUE 377 C F I R S T SCALE THE POINTS 3/8 CALL S C A L E ( X , N , 1 0 . 0 , X M I N , D X , 1 ) 379 CALL S C A L E ( Y , N , 1 0 . 0 , Y M I N , D Y , 1 ) 380 CALL AX I S ( 0 . , 0 . , 'TIME S H I F T ' , - 9, 1 0 . 0 , 0 . ,XMIN,DX) 381 ' CALL A X I S 1 0 . ,0. , 'CROSS CORRELATION VALUE' , 9, 10 . 0 , 9 0 . ,YMIN,DY) 382 C F I N A L L C PLOT THE L I N E S 383 CALL L I N E ( X , Y,N,1) 384 C MOVE THE O R I G I N AND RETURN 385 RETURN 386 END 387 c 388 SUBROUTINE BNDPAS(F1,F2,DELT,D,G) 389 c 390 c S u b r o u t i n e by Dave G a n l e y on M a r c h 5, 1977. 391 c 392 c 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 393 c r e c u r s i v e B u t t e r w o r t h b a n d p a s s f i l t e r . I n o r d e r t o d e s i g n 394 c t h e f i l t e r a c a l l must be made t o BNDPAS a n d t h e n t h e 395 c f i l t e r may be a p p l i e d by c a l l s t o f i l t e r . The f i l t e r 396 c w i l l h a v e 8 p o l e s i n t h e s p l a n e a n d i s i n f o r w a r d a n d 397 c r e v e r s e d i r e c t i o n s s o a s t o h a v e z e r o p h a s e c u t o f f 398 c f r e q u e n c i e s . The a t t e n u a t i o n w i l l be -6DB a n d t h e r o l l o f f 399 c w i l l be a b o u t 96 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 400 c u s e d i n d e s i g n i n g t h e f i l t e r t o p r e v e n t a l i a s i n g p r o b l e m s . 401 c 402 COMPLEX P ( 4 ) , S ( 8 ) , Z 1 , Z 2 403 DIMENSION D ( 8 ) , X ( 2 0 4 8 ) , X C ( 3 ) , X D ( 3 ) , X E ( 3 ) 404 DATA ISW/0/,TWOPI/6.2831853/ 405 c 406 c T h i s s e c t i o n c a l c u l a t e s t h e f i l t e r a n d must be c a l l e d 173 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27, 1988 f o r C C i d = S E I S on G 407 C b e f o r e f i l t e r i s c a l l e d . 408 C 409 C F1 = Low f r e q u e n c y c u t o f f (6 DB Down) 410 C F 2 = H i g h f r e q u e n c y c u t o f f (6 DB down) 411 C DELT = S a m p l e i n t e r v a l i n m i l l i s e c o n d s 412 C D = W i l l c o n t a i n 8 Z d o m a i n c o e f f i c i e n t s o f r e c u r s i v e f i l t e r 4 13 C G = W i l l c o n t a i n t h e g a i n o f t h e f i l t e r 4 14 C 415 W R I T E ( 6 , 1 ) F 1 , F 2 , D E L T 416 1 FORMAT!'BANDPASS F I L T E R DESIGN FOR A BAND FROM',F8.3,'TO',F8 417 1.3,'HERTZ.',11' SAMPLE INTERVAL I S ' , F 5 . 2 , ' M I L L I S E C O N D S . ' ) 418 DT = D E L T / 1 0 0 0 . 0 419 TDT = 2.0/DT 420 FDT = 4.0/DT 421 ISW = 1 422 P ( 1 ) = C M P L X ! - . 3 8 2 6 8 3 4 , . 9 2 3 8 7 9 5 ) 423 P ( 2 ) = C M P L X ! - . 3 8 2 6 8 3 4 , - . 9 2 3 8 7 9 5 ) 424 P ( 3 ) = C M P L X ! - . 9 2 3 8 8 7 9 5 , . 3 8 2 6 8 3 4 ) 425 P ( 4 ) = C M P L X ( - . 9 9 2 3 8 7 9 5 , - . 3 8 2 6 8 4 ) 426 W1 = TWOPI'FI 427 W2 = TWOPI*F2 428 W1 = TDT * TAN(W1/TDT) 429 W2 = TDT*TAN(W2/TDT) 430 HWID = (W2-WU/2.0 431 WW = W1'W2 432 DO 19 1=1,4 433 Z1 = P ( I ) ' H W I D 434 Z2 = Z1'Z1-WW 435 Z2 = C S Q R T ( Z 2 ) 436 S ( I ) = Z1 + Z2 437 19 S ( I + 4 ) = Z1-Z2 438 C W R I T E ( 6 , 2 ) S 439 C 2 FORMAT('-S PLANE POLES ARE AT: ' ,/' ' , 8 ( / ' ' , E 1 2 . 6 , ' +I ' ,E 1 2 . 6 ) ) 440 G = .5/HWID 441 G = G"G 442 G = G" G 443 DO 29 1=1,7.2 444 B = - 2 . 0 * R E A L ( S ( I ) ) 445 Z1 = S ( I ) * S ( I + 1 ) 446 C = R E A L ( Z 1 ) 447 A = TDT+B+C/TDT 448 G = G * A 449 D ( I ) = ( C ' D T - F D T l / A 450 29 D( 1 + 1 ) = ( A - 2 . 0 ' B ) / A 451 G = G * G 452 C WRITEC6.3) 453 C 3 FORMAT!'-FILTER I S ( 1 - Z * ' 2 ) " * 4 / B 1 * B 2 * B 3 * B 4 ' ) 454 C W R I T E ( 6 , 4 ) D 455 C 4 F 0 R M A T ( 4 ( / ' B ( I ) = 1 + ' .E12.6, ' Z+' .E12.6, ' Z " 2 ' ) ) 456 C W R I T E ( 6 , 5 ) G 457 C 5 FORMAT!'-FILTER GAIN I S ' . E 1 2 . 6 ) 458 RETURN 459 C 460 C 461 ENTRY F I L T E R ( X , N , D , G , I G ) 462 C 463 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 f i l t e r e d 464 C 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 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27, 1988 f o r C C i d = S E I S on G 465 C G = F i I t e r G a i n 466 C IG = 1 Means t o r e m o v e t h e f i l t e r g a i n s o t h a t 467 C i s u n i t y . 468 C 469 I F ( I S W . E Q . 1 ) GO TO 31 470 W R I T E ( 6 , 6 ) 47 1 6 FORMATf'1BNDPAS MUST BE CALLED BEFORE F I L T E R ' ) 472 CALL E X I T 473 C 474 WRITE(6 , " ) X( 1 ) 475 C A p p l y f i l t e r i n f o r w a r d d i r e c t i o n 476 C 477 31 XM2=X(1) 478 XM1=X(2) 479 XM=X(3) 480 XC(1)=XM2 481 X C ( 2 ) = X M 1 - D ( 1 ) * X C ( 1 ) 482 X C ( 3 ) = X M - X M 2 - D ( 1 ) * X C ( 2 ) - D ( 2 ) * X C ( 1 ) 483 X D ( 1 ) = X C ( 1 ) 484 X D ( 2 ) = X C ( 2 ) - D ( 3 ) ' X D ( 1) 485 X D ( 3 ) = X C ( 3 ) - X C ( 1 ) - D ( 3 ) ' X D ( 2 ) - D ( 4 ) * X D ( 1 ) 486 X E ( 1 ) = X D ( 1 ) 487 X E ( 2 ) = X D ( 2 ) - D ( 5 ) * X E ( 1 ) 488 X E ( 3 ) = X D ( 3 ) - X D ( 1 ) - D ( 5 ) * X E ( 2 ) - D ( 6 ) ' X E ( 1 ) 489 X ( 1 ) = X E ( 1) 490 X ( 2 ) = X E ( 2 ) - D ( 7 ) * X ( 1) 49 1 X ( 3 ) = X E ( 3 ) - X E ( 1 ) - D ( 7 ) * X ( 2 ) - D ( 8 ) , X ( 1 ) 492 DO 39 1=4,N 493 XM2=XM1 494 XM1=XM 495 XM=X(I) 496 K = I - ( ( I - 1 )/3) *3 497 GO TO ( 3 4 , 3 5 , 3 6 ) ,K 498 34 M=1 499 M1 = 3 500 M2=2 501 GO TO 37 502 35 M=2 50 3 M1 = 1 504 M2 = 3 505 GO TO 37 506 36 M=3 507 M1 = 2 508 M2=1 509 37 X C ( M ) = X M - X M 2 - D ( 1 ) " X C ( M 1 ) - D ( 2 ) ' X C ( M 2 ) 510 X D ( M ) = X C ( M ) - X C I M 2 ) - D ( 3 ) * X D ( M 1 ) - D ( 4 ) * X D ( M 2 ) 51 1 X E ( M ) = X 0 ( M ) - X D ( M 2 ) - D ( 5 ) ' X E ( M 1 ) - D ( 6 ) " X E ( M 2 ) 512 39 X ( I ) = X E ( M ) - X E ( M 2 ) - D ( 7 ) ' X ( I - 1 ) - D ( 8 ) ' X ( I - 2 ) 513 C 514 c F i l t e r i n r e v e r s e d i r e c t i o n 515 C 516 XM2=X(N) 517 XM1=X(N-1) 518 XM=X(N-2) 519 XC( 1 )=XM2 520 X C ( 2 ) = X N 1 - D ( 1 ) " X C ( 1 ) 521 X C ( 3 ) = X M - X M 2 - D ( 1 ) ' X C ( 2 ) - D ( 2 ) ' XC( 1) 522 X D ( 1 ) = X C ( 1 ) 175 L i s t i n g 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 523 X D ( 2 ) = X C ( 2 ) - D ( 3 ) * X D ( 1 ) 524 X D ( 3 ) = X C ( 3 ) - X C ( 1 ) - D ( 3 ) ' X D ( 2 ) - D ( 4 ) " X D ( 1 ) 525 X E ( 1 ) = X D ( 1 ) 526 X E ( 2 ) = X D ( 2 ) - D ( 5 ) ' X E ( 1) 527 X E ( 3 ) = X D ( 3 ) - X D ( 1 ) - D ( 5 ) ' X E ( 2 ) - D ( 6 ) * X E ( 1 ) 528 X ( N ) = X E ( 1) 529 X ( N - 1 ) = X E ( 2 ) - D ( 7 ) 'X( 1 ) 530 X ( N - 2 ) = X E ( 3 ) - X E ( 1 ) - D ( 7 ) , X ( 2 ) - D ( 8 ) ' X ( 1 ) 531 DO 49 1=4,N 532 XM2=XM1 533 XM1=XM 534 J=N-I+1 535 XM=X(J) 536 K = I - ( ( I - 1 ) / 3 ) * 3 537 GO TO ( 4 4 , 4 5 , 4 6 ) ,K 538 44 M=1 539 M1 = 3 540 M2=2 54 1 GO TO 47 542 45 M-2 543 M1 = 1 544 M2=3 545 GO TO 47 546 46 M=3 547 M1=2 548 M2=1 549 47 X C ( M ) = X M - X M 2 - D ( 1 ) ' X C ( M 1 ) - D ( 2 ) * X C ( M 2 ) 550 X D ( M ) = X C ( M ) - X C ( M 2 ) - D ( 3 ) * X D ( M 1 ) - D ( 4 ) * X D ( M 2 ) 551 X E ( M ) = X D ( M ) - X D ( M 2 ) - D ( 5 ) " X E ( M 1 ) - D ( 6 ) * XE(M2) 552 49 X ( J ) = X E ( M ) - X E ( M 2 ) - D ( 7 ) * X ( J + 1 ) - D ( 8 ) * X ( J + 2 ) 553 I F ( I G . N E . I ) RETURN 554 DO 59 1 = 1 ,N 555 59 X ( I ) = X ( I ) / G 556 W R I T E ( 6 , " ) ( X ( I ) , 1=1,N) 557 RETURN 558 END 559 C 560 C S u b r o u t i n e DOT f i n d s t h e d o t p r o d u c t o f two a r r a y s . 56 1 C 562 SUBROUTINE DOT(L,X,Y,ANS) 563 C 564 I M P L I C I T REAL (A-H.O-Z) 565 REAL X ( 2 ) , Y ( 2 ) , A N S 566 ANS=0 .0 567 I F ( L ) 3 0 , 3 0 , 1 0 568 10 DO 20 1 = 1 ,L 569 20 ANS=ANS + X( I ) * Y( I ) 570 30 RETURN 571 END 572 C 573 C S u b r o u t i n e CROSS f i n d s t h e c r o s s p r o d u c t o f two a r r , 574 C 575 SUBROUTINE C R O S S ( L X , X , L Y , Y , L D , D ) 576 C 577 I M P L I C I T REAL (A-H.O-Z) 578 REAL X ( 2 ) , Y ( 2 ) , D ( 2 ) 579 DO 10 1 = 1 , LD 580 10 CALL D 0 T ( M I N 0 ( L Y + I - 1 , L X ) - I + 1 , X ( I ) , Y , D ( I ) ) 176 L i s t i n g o f CROSSCOR a t 16:44:17 on JUL 27, 1988 f o r C C i d = S E I S on G 581 RETURN 582 END 583 C 584 C S u b r o u t i n e MAXSN f i n d s t h e maximum v a l u e o f an a r r a y , 585 C t a k i n g i n t o a c c o u n t t h e a l g e r b r a i c s i g n s o f t h e e l e m e n t s . 586 C 587 SUBROUTINE MAXSN(LX,X,XMAX,INDEX) 588 C 589 REAL X ( 2 ) , X M A X 590 XMAX=X(1) 59 1 DO 10 1 = 1,LX 592 10 XMAX=AMAX1(XMAX,X(I)) 593 DO 20 J=1,LX 594 INDEX=J 595 I F ( X ( J ) - X M A X ) 2 0 , 3 0 , 2 0 596 20 CONTINUE 597 30 RETURN 598 END 177 APPENDIX (3 Program l i s t i n g of KALDAT, the d iscrete model of the second order system. 178 L is t ing of KALDAT at 18:25:29 on DEC 12. 1987 for CCId-SEIS 1 C This program generates synthetic geophone data Inorder 2 C to determine the performance of the Kalman f i l t e r . The Input 3 C noise, wk. is white, and has a gausslan d i s t r i b u t i o n with 4 C mean=0.0 and standard d e v l a t 1 o n « S I G M A 5 C 6 C 7 REAL A 11.A12.A21.A22.B11.821.WO.R0.PO.R01.P01 8 REAL TC.GAMA.HO.HI,H2,KI,WSIGMA,WMEAN,WS.ALPH,DELT 9 REAL DM,012,VS,AW, BW, OA ,SSS,SA.SS.BB11,BB1,ST,SH,SP 10 REAL Xl(40OO).X2(40OO).X3(4OOO).X4(40OO).X5(400O) 11 REAL X6(4000),X7(4000).X8(4000).X9(4000).2(4000) 12 REAL X81(40O0),X82(40O0).X91(4000).XS2(4000).V(4000) 13 INTEGER T 14 C 15 C 16 PIE=3.1415927 17 C 18 C length of desired time ser ies 19 T = 200 20 C 21 C sampling Interva1.de1ta 22 DELT=.01 23 C 24 C standard deviat ion of noise present 25 WSIGMA=.09 26 VSIGMA=.009 27 C 28 C natural frequency of Insrument 29 WO=8»2»PIE 30 C 31 C damping of Instrument 32 R0=.IB 33 C 34 C KI factor . 35 KI=1.0 36 C 37 c gain on measurement H 38 HO=.29 39 c 40 c time constant 41 T O . 10 42 c 43 DO IOO I » 1 . T 44 READ( 15.•) X1(I),X3(I ) 45 100 CONTINUE 46 c 47 WMEAN-0.0 48 GAMA- 1/TC 49 ALPH*ASIN(RO) 50 R01=R0»W0*DELT 51 P0«S0RT(1 -R0«R0) 52 POI-PO'WO'DELT 53 AA = EXP(-R01 )/P0 54 A11=AA« (C0S(P01-ALPH) ) 55 A !2-(AA*SIN(P01))/WO 56 A 2 1 » - A A » ( W 0 » S I N ( P 0 1 ) ) 57 A22"AA*(COS(PO1+ALPH)) 59.5 B 11-( 1 0 -A11) »K I/ (W0*W0) 179 L i s t i n g of KALOAT at 18:25:29 on DEC 12. 1987 for CCid=SEIS 59.7 B12=-A21'KI/(WO*WO) 60 H1«-W0*W0 6 1 H2=>-2»R0«WO 62 QA=GAMA*DELT 63 AW=£XP( -OA) 64 BW=1.0 -AW 65 X5(1)=0.0 66 X6(1)=0.0 67 X7(l)=0.0 68 X8(l)=0.0 70 X9(1)=0.0 72 2(1)=0.0 73 WRITE(6.') AW.BW 74 C 75 00 1 I - 1 . T 76 C 77 c At th i s point i t is necessary to generate the Gaussian noise 78 C as input input our f i r s t order transform funct1on(which represents 79 c a low pass f i l t e r ) which has X5 as our output. X5 is then feed 80 c into our second order model. 8 1 c 82 TK=I 83 ST=1 .0 «1 84 SH=2.0'I 85 SS=SQRT(D£LT) 85 . 5 SP=GAMA/2.0 86 c 87 CALL RNGAUS(WS.WMEAN,WSIGMA.TK,ST) 88 89 NN=I+1 90 c 91 X5(NN)=X5( I ) »AW+BW»WS 92 c 93 c 94 X6(NN ) = B 11 *X5(I) + A 11 *X6(I) + A12*X7(I) 95- c 96 c 97 X7(NN)=B12*X5(I ) + A2 1*X6(I ) + A 2 2 » X 7 ( I ) 98 c 99 BB 1t*Bl1*X1(I) 10O BB1=B1 1'X3(I ) 101 c 102 X8(NN)=A11*X8(I ) + A12*X9(I)+BB11+BB1 105 106 c 107 c 108 X 9 ( N N ) i A 2 1 * X 8 ( I ) + A 2 2 * X 9 ( I ) + B 1 2 * X l ( I ) + B 1 2 » X 3 ( I ) 1 1 1 c 112 c 113 CALL RNGAUS(VS.WMEAN.VSIGMA.TK.SH) 1 14 c 1 15 Z ( N N ) = H 0 « ( X 7 ( N N ) + X 9 ( N N ) ) + V S / S S 1 16 c 1 17 WRITE(2. • ) XS(I ) 1 18 c 1 19 1 CONTINUE 120 c 12 1 DO 2 J=1.T 180 L is t ing KALDAT at 18:25:29 on DEC 12. 1987 for CCId-SEIS 122 WRITE(1.•) Z(J) 123 2 CONTINUE 124 C 125 C 126 C DO 2 J - l . T 127 C W R I T E ( 2 . » ) X5(J)+1 .4 128 C 2 CONTINUE 129 C 130 STOP 131 END 132 C 133 SUBROUT INE RNGAUS ( W. W.ME AN. WS IGMA . TK . SS ) 134 C 135 REAL RNSEO,WMEAN,C1,C2.C3,WSIGMA,W,SS 136 INTEGER TK 137 DATA IRNSED/1/ 138 DATA TW0P1/6.283185307/ 139 C 140 C IF (TK.LT.10.0) THEN 141 C SS-TK*1000.0+1.0 142 C ELSE 143 C SS-TK»100.0+1 144 C END IF 145 C 146 C 147 GO TO (10.20.30) IRNSED 148 C 149 10 RNSED-ABS(RNSED) 150 CI-RAND(RNSED) 151 C 152 C F i r s t c a l l for a random number. 153 C 154 20 C1=S0RT(-2.*AL0G(RAN0(O.))) 155 C2 -TW0PI »RAND(0 . ) 156 C 3 ' C 1 » C 0 S ( C 2 ) 157 IRNSED-3 158 GO TO 40 159 C 160 C Second c a l l for random number. 161 C 162 30 C 3 - C 1 « S I N ( C 2 ) 163 IRNSED-2 164 C 165 40 W-WMEAN+WSIGMA-C3 166 C 167 RETURN 168 END 181 APPENDIX C Program l i s t i n g of the Kalman F i l t e r , KALMAN 4 182 L i s t i n g . KALMAN4 at 10:07:47 on DEC 13. 1987 for CCId-SEIS 1 C 2 C 3 £•»•*•*••*»•••*»»»•'****»***«*•*********•**••*•#««**•*••••**«••*••»#** 4 C* • 5 C» MAIN ' PROGRAM * 6 C* • 7 C* • 8 C* • 9 C« • 10 C* • 11 C« . • j2 £••••«••••»*•»•***••»•*•*****••*•*••»*•»•••*••••*•»**••»••»*•#«»•**•» 13 C 14 C 15 C 16 INTEGER SNP.SN.ST.SI.SP.NT.FT,KT.QT.T.IDEX4.IDEX2 16.5 INTEGER TE(4) 16.7 REAL X2.X22.X4.X44 17 C 18 REAL WO.RO.KI.BB.LL.RR.R.0ELT.A11.A21.A12.A22 19 REAL AA.ALPH.ESP.B1t,B12.0R.GAMA.C1.C2.C3.C4,P01.R01 20 REAL PIE.AN.RW.HH.TC,AW.QW.RA.PW.MV,HV(9) 21 REAL XMIN1.XMIN2,XMIN3.XMIN4.XMAX4.XMAX2.SNN.SPP 22 REAL PH1(500).PH2(500).PH3(500),PH4(500),XH1(500).XH2(SOO) 23 REAL XH3(500).XH4(50O).KK.QD.N1.NPP.NSS.SD.N3 24 INTEGER INOEX1.INDEX2.INDEX3.IN0EX4,N2.N4,NP(6).NS(6) 25 C 26 REAL Z(2000).WP(200).WS(200).BB!1.BB21.BB1.BB2 27 REAL KS(9).KH(9.9).BETA(10).QN(2.2).PS0(9.9) 28 REAL PS(9,9) .PS1(4.4) .PS2(4.5) .PS2T(5.4) .PS3(5.5) 29 REAL PU(9.9).PU1(4.4).PU2(4.5),PU2T(5.4),PU3(5,5) 30 REAL FS(9.9) .FS1(5.4) ,FS2(5.5) ,XU(9),XE(9).X0(9) 31 REAL FL(9.9).FLl(5.4) .FL2(5.5).H11(200).H12(200) 32 C 33 C 34 C Defining parameters: 35 C 36 C WO=natural frequency of Instrument R0=damplng of Instrument 37 C KI= coeff1cIent which takes into account instrument inductance. 38 C resistance and magnetic f i e l d . 39 C 0ELT=samp1tng Interval AN=white noise Input into second order 40 C system. 41 C ER=covariance on measurement noise QR=covarlance on state noise 42 C 43 C i n i t i a l estimates of the states 44 C 45 REAO(11.*)X0(1),X0(2),X0(3),X0(4).X0(5).X0(6).X0(7).X0(8) 46 1 .X0(9).FT 47 C 48 C 49 DO 76 1=1.9 50 DO 77 J=1.9 51 PS0(I.J)=0.0 52 77 CONTINUE 53 76 CONTINUE 54 C 55 READ(11.•)PS0(1.1).PS0(2.2).PSO(3.3).PS0(4,4).PSO(5.5). 56 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 for CC1d*SEIS 57 1PS0( 1 .8) .PS0(1.9) .PS0(2.B) .PS0(2 58 1PS0(3.8) .PS0(3.9) .PS0(4.8) .PS0(4 59 C eo C 61 PS0(8.1) -PS0(1.8) 62 PS0(9.1) -PS0(1.9) 63 P S 0 ( 8 . 3 ) » P S 0 ( 3 . 8 ) 64 PS0(9.3) -PS0(3.9) 65 P S 0 ( 8 . 4 ) » P S 0 ( 4 . 8 ) 66 PS0(9.4)=PS0(4.9) 67 PS0(8.2) -PS0(2.8) 68 PS0(9.2) -PS0(2.9) 69 P S 0 ( 7 . 6 ) « P S 0 ( 6 . 7 ) 70 PS0(9.8)=PS0(8.9) 7 1 C 72 DELT-.01 73 C 74 K I -1 .0 75 C 76 T = 200 77 C 78 PIE=3.1415927 79 C 80 W0=8»2*PIE 8t C B2 R0 =.18 83 C 84 C NP: length of Input P-wave 85 NS( 1) = 1 86 C NP(2)=2 87 C NP(3)=3 88 NSS=1 89 C 90 C NS: length of Input S-wava 91 DO 232 1 = 1. 10 92 NP(I)=I 93 C NS(3)=3 94 232 CONTINUE 95 NPP=9 96 c 97 c WP: c o e f f i c i e n t s de f in ing P-wave 98 WS( 1 )=0.0 99 c WP(2)=1.0 too c WP(3)=0.0 101 c WP(4)=0.0 102 c WP(5)-.6 103 c WP(6)=.8 104 c 105 c WS: c o e f f i c i e n t s de f in ing S-wave 106 00 323 1 = 1. 10 107 F1=50.0 108 QAS= .01 *(I ) 109 W P ( I ) = C 0 S ( F 1 » Q A S ) * E X P ( - Q A S ) 1 10 323 CONTINUE 1 1 1 c WS(2)=1.0 1 12 c W S ( 3 ) « 0 . 0 113 c WS(4)=0.0 1 14 . c WS(5)=.6 184 L i s t i n g di KALMAN4 at 10:07:47 on DEC 13. 1987 for CCId-SEIS 115 C WS(6)=8 116 C 1 17 C I 1 18 C MV: gain or ampl i f i cat ion on recorder 119 MV=.29 120 C 121 C FT: number of Iterations desi red 122 C FT = 4 123 C 124 00 78 1=1.T 125 R E A D ( 8 . « ) Z(I) 126 78 CONTINUE 127 C 128 C TC: time constant 129 TC=.10 130 C 131 C 132 ESP=1/TC 133 H1=-W0'W0*1.0 134 KI=1.0 135 H2=-2»RO»WO*1.0 136 P0=S0RT(1 -R0»R0) 137 R01«R0»WO*DELT 138 P01=P0'WO»DELT 139 GAMA=ATAN(RO/PO) 140 ALPH=ASIN(RO) 141 AA=EXP(-R01)/P0 142 A 1 2 « ( A A ' S I N ( P 0 1 ) ) / W 0 143 v A 2 1 = - A A « ( W 0 » S I N ( P 0 1 ) ) 144 A22=AA»(C0S(P01+ALPH)) 145 A11=AA« (C0S(P01 -ALPH) ) 146 B1 1 = (1 -A11) »K I/ (WO*WO) 147 B12=-A21«KI/(WO*WO) I 148 C 1 = ( ( V O * W O ) » E X P ( - R 0 1 ) ) / P 0 149 C2=SIN(P01)/W0 150 C3=C0S(P0t+ALPH) 151 C4=P0»EXP( -R01) 152 AW=EXP(-ESP'DELT) 153 BW=1.0-AW 154 0R=.O9 ! 155 0W=BW*BW»0R | 156 RW=.009 157 C 158 DO 9 1=1.9 159 XE(I)=X0(I) 160 9 CONTINUE 161 C 162 DO 1 111=1.FT 163 C 163.5 TEC 1 ) = IFiX(XE( 1 )) 163.7 TE(3) = IFIX(XE(3) ) 163.75 XE( 1 ) = FLOAT(TE( 1 ) ) 163.77 XE(3)=FL0AT(TE(3) ) 163.8 C 163.9 C 164 X2=XE(2)/OELT 165 TE(2)=IFIX(X2) 166 X22=TE(2)'DELT 185 L i s t i n g of KALMAN4 at 10:07:47 on DEC 13. 1987 for CCid-SEIS 167 XE(2)-X22 167.5 C 168 X4-XE(4)/DELT 169 TE(4)-IFIX(X4) 170 X44-TE(4)'DELT 171 X E ( 4 ) » X 4 4 171.5 C 171.6 C W R I T E ( 6 . » ) XE(1).XE(2 ) .XE(3).XE(4) 171.7 C 172 XE(5)=X0(5) 173 XE(6)=X0(6) 174 XE(7)=XO(7) 175 X E ( 8 ) « X 0 ( 8 ) 176 . X c ( 9 ) » X 0 ( 9 ) 177 C 178 DO 70 1 -1.9 179 DO 80 J - 1 . 9 180 PS(I .J )=PS0(I .J) 181 80 CONTINUE 182 70 CONTINUE 183 C 184 DO 30 1=1.4 185 DO 31 J - 1 . 4 186 PS1(I .J)=PS(I .J) 187 31 CONTINUE 188 30 CONTINUE 189 C 190 DO 32 1=5.9 191 DO 33 J - 5 . 9 192 D = d-4 193 D0=I-4 194 PS3(DD,D)=PS(I.J) 195 33 CONTINUE 196 32 CONTINUE 197 C 198 DO 34 1=1.4 199 DO 35 d - 5 . 9 200 D=J-4 201 PS2( I . D ) » P S ( I .'J) 202 35 CONTINUE 203 34 CONTINUE 204 C 205 C 206 DO 17 1=1.9 207 00 18 0 - 1 . 9 208 FS(I .J )=0.0 209 FL( I . J )=0.0 210 18 CONTINUE 211 17 CONTINUE 212 - C 213 C 214 FS(5.5)=AW 215 FS(6.5)=811 216 FS(6.6)=A11 217 FS(6.7)-A12 218 FS(7.5)=B12 219 FS(7.6)=A21 220 FS(7.7)=A22 186 L i s t i n g u. KALMAN4 at 10:07.47 on DEC 13. 1987 for CC1d=SEIS 221 FS(8.8)=A11 222 FS(8.9)=A12 223 FS(9.8)=A21 224 FS(9.9)=A22 225 C 226 DO 21 1*1.5 227 DO 22 J - 1 . 5 228 FL( I .U) -FS( I . J ) 229 22 CONTINUE 230 21 CONTINUE 231 C 232 DO 23 1-1.4 233 FS<I.I )- 1 .0 234 FL(I . I )=1.0 235 23 CONTINUE 236 C 237 SI » 0 238 SP = 0 239 C 240 C 241 C 242 DO 2 11 = 1 .T 243 C 244 KT=II 245 C 246 C 247 C Lines 328-356 define the error t r a n s i t i o n matrix. BS, which 248 C has nonnutl components defined by the d iscrete second order 249 . C system equat1ons(I.e. B1I.B12). This matrix Is a lso a component 250 C of the state t rans i t ion matrix, F, and is a function of the P and 251 C S-wave time of a r r i va1s ( tp -X2 . t s -X4 ) . 252 C 253 C 254 C 255 N1=NPP»DELT 256 C 257 0D=XE(2)+N1 258 KK=KT *DELT 259 IF (KK.GE.XE(2) OR.KK.GE.XE(4) ) THEN 260 F S ( 8 . a ) » A t 1 261 FL(8.8)=A11 262 FS(8.9)=A12 263 FL(8.9)=A12 264 FS(9.8)=A21 265 FL(9.8)-A21 266 Fi(9.9)=A22 267 FL(9.9)=A22 268 C 269 C 270 ELSE 271 GOTO 109 272 C 273 END IF 274 C 275 109 IF (KK.LT.XE(2).0R.KK.GT.0D) THEN 276 C 277 BB11=0.0 278 BB21=0.0 187 L i s t i n g of KALMAN4 at 10:07:47 on DEC 13, 1987 for CCIOSEIS 279 GOTO 173 280 C 281 ELSE 282 C 283 SI=SI+1 284 C 285 IF (SI.GT.NPP+1) THEN 286 GOTO 173 287 ELSE 288 GOTO 174 289 END IF 290 174 N2=NP(SI) 291 BB 1 1 *B 1 1*WP(N2) 292 BB21=B12*WP(N2) 293 C 294 173 END IF 295 C 296 C 297 N3=NSS*DELT 29B SD=XE(4)+N3 299 C 300 IF (KK.LT.XE(4).OR.KK.GT.SD ) THEN 301 C 302 BB1=0.0 303 BB2=0.0 304 C 305 ELSE 306 SP=SP+1 307 IF (SP.GT.NSS+1) THEN 308 GOTO 171 309 ELSE 310 GOTO 172 311 END IF 312 172 N4=NS(SP) 313 BB1=B11*WS(N4) 314 BB2»812«WS(N4) 315 316 C W R I T E ( 6 . » ) KK.XE(4).SD.WS(N4) 317 C • 318 17 1 END IF 319 C 320 C WRITE(6.•)BB11,BB21.BB1,BB2 321 C 322 FS(8. 1 )°BB 1 1 323 FS(9.1)=BB21 324 FS(8.3)=BB1 325 FS(9.3)=BB2 326 FL(8. 1)«BB1 1 327 FL(9.1)=BB21 328 FL(8.3)=BB1 329 FL(9.3)=8B2 330 C 331 C WRITE(6.*) FS(8. 1 ) , FS (9 .1 ) . FS (8 .3 ) , FS (9 .3 ) 332 C 333 DO 303 1=1.9 334 HV(I)=0.0 335 303 CONTINUE 336 C 188 L i s t i n g o. KALMAN4 at 10:07:47 on DEC 13. 1987 for CCId-SEIS 337 HV(7)=MV 338 HV(9)-MV 339 C 340 CALL XUPDAT(XU.XE.FS) 341 C 342 C 343 CALL EXTEND(FL.C1,C2.C3.C4.XU.KK.N1,N2.WS.WP.DELT.N3,N4.FT) 344 C 345 C 346 C 347 CALL EXTRAP(PU.PS.FL) 348 C 349 C 350 P J ( 5 . 5 ) » P U ( 5 . 5 ) + Q W 351 C 352 CALL GAIN(HV.PU.RW.KS) 353 C 354 CALL NWSTAT(XE.KS.XU.MV.KT.Z) 355 C 356 CALL PUPDAT(HV.KS.PU.PS) 357 C 358 C 359 P H 1 ( K T ) « P S ( 1 . 1) 360 PH2(KT)*PS(2,2) 361 PH3(KT)=PS(3.3) 362 PH4(KT)=PS(4.4) 363 XH1(KT)=XE(1) 364 XH2(KT)=XE(2) 365 XH3(KT)=XE(3) 366 XH4(KT)=XE(4) 367 C 368 WRITE(1.*) PHI(KT) 369 W R I T E ( 2 . « ) PH2(KT) 370 W3ITEO.*) PH3(KT) 371 W R I T E ( 4 . « ) PH4(KT) 372 W R I T E ( 5 , « ) XE(1) 373 W R I T E ( 7 , « ) XE(2) 374 W R I T E ( 1 0 . » ) XE(3) 375 W R I T E ( 9 . » ) XE(4) 376 C 377 2 CONTINUE 378 C 379 C 380 C 381 CALL MINSN(T.PHI.XMINI.INOEXI) 382 C W R I T E ( 6 . « ) XMIN1.XH1(INDEX 1).INDEX 1 383 C 384 XE( 1) »XH1( INDEX 1 ) 3B5 C 386 CALL MINSN(T.PH2.XMIN2.INDEX2) 387 XE(2)=XH2(INDEX2) 388 C 389 CALL MINSN(T.PH3.XMIN3.INDEX3) 390 C W R I T E ( 6 . « ) XMIN3,XH3(INDEX3).INDEX3 391 C 392 XE(3)* XH3(INDEX3) 393 C 394 CALL MINSN(T.PH4.XMIN4.IN0EX4) 189 L i s t i n g of KALMAN4 at 10:07:47 on DEC 13. 1987 for CCfd=SEIS 395 X E ( 4 ) » X H 4 ( I N D E X 4 ) 396 C W R I T E ( 6 . » ) XMIN4.XE(4).INDEX* 397 C 398 CALL MAXSN(T.PH2.XMAX2.I0EX2) 399 C WRITE(6.•) XMAX2.XH2(IDEX2+2).IDEX2 400 C 401 C 402 CALL MAXSN(T.PH4.XMAX4.I0EX4) 403 C W R I T E ( 6 . » ) XMAX4.XH4(IDEX4).IDEX4 404 C 405 C IF (INDEX2.E0.1) THEN 406 C 407 C X E ( 2 ) » X H 2 ( I D E X 2 ) 408 C 409 C ELSE 410 C 4 1 1 C XE(2)=(XH2(INDEX2+2)+XH2(INDEX2)+XH2(INDEX2+1))/3 412 C 4 13 C END IF 4 14 C 4 15 C 4 16 C 417 c IF (INDEX4.E0.1) THEN 418 c 419 c XE(4)=(XH4(IDEX4+2)+XH4(IDEX4)*XH4(IDEX4+1))/3 420 c XE(4)=XH4(IDEX4) 421 c ELSE 422 c 423 c X E ( 4 ) » ( X H 4 ( I N D E X 4 + 2 ) + X H 4 ( I N D E X 4 ) + X H 4 ( I N D E X 4 + 1 ) ) / 3 424 c 425 c END IF 426 c 427 c 428 c S N N » X E ( 2 ) / D E L T 429 c SN=IFIX(SNN) 430 c XE(2)=SN*0ELT 431 c W R I T E ( 6 . » ) XE(2) 432 c 433 c SPP=XE(4)/DELT 434 c SNP=IFIX(SPP) 435 c XE(4)=SNP*DELT 436 c WRITE(6.*) XE(4) 437 c 438 W R I T E ( 6 . » ) XMIN2.XE(2).INDEX2 439 c 440 W R I T E ( 6 . » ) XMIN1.XE(1).INDEX 1 44 1 c 442 1 CONTINUE 443 c 444 STOP 445 END 446 c 447 c 448 c This subroutine extrapolates the state matrix X(I ) . 449 c 450 SUBROUTINE XUPDAT(XU.XE.FS) 451 c 452 REAL XU(9) ,XE(9) ,FS(9.9) 190 L i s t i n g KALMAN4 at 10:07:47 on DEC 13. 1987 for CCid=>SEIS 453 INTEGER OR 454 C 455 OR • 9 456 C 457 CALL GMATV(FS.XE.XU.OR.OR.OR) 458 C 459 RETURN 460 ENO 461 C 462 C 463 C This subroutine propagates and updates the Kalman gain matrix K. 464 C which is a function of P and S-wave a r r i v a l s . 465 C 466 SUBROUTINE GAIN(HV.PU.RW.KS) 467 C • 468 REAL RW,PU(9,9),KS(9),HV(9).TMP1(9) 469 REAL D0T.TMP2(9) 470 INTEGER OR 47 1 C 472 OR = 9 473 C 474 CALL GVMAT(HV.PU.TMP1.OR.OR.OR) 475 C 476 DOT » GVV(TMP1.HV.OR) • RW 477 C W R I T E ( 6 . » ) DOT 478 C 479 CALL GMATV(PU.HV.TMP2.0R,0R.0R) 480 C 481 DO 1 1=1.9 482 K S ( I ) » T M P 2 ( I ) / D O T 483 1 CONTINUE 484 c 485 RETURN 486 END 487 c 488 c 489 c This subroutine gives us s tate estimate updates. 490 c 49 1 SUBROUTINE NWSTAT( XE,KS,XU.MV,KT,Z) 492 c 493 REAL DELTZ.Z(200O).XU(9).XE(9) 494 REAL KS(9).MV 495 INTEGER KT 496 497 c D E L T Z ( K T ) « Z ( K T ) - H ( 1 ) » X 1 ( K T ) - H ( 2 ) • X E ( 2 ) - H ( 3 ) » X 3 ( K T ) - H ( 4 ) » X E ( 4 ) 498 c 1 - H ( 5 ) * X 5 ( K T ) - H ( 6 ) * X 6 ( K T ) - H ( 7 ) » X 7 ( K T ) - H ( 8 ) « X 8 ( K T ) - H ( 9 ) » X 9 ( K T ) 499 c 500 DELTZ = Z(KT)-MV(XU(7)+XU(9) ) 501 c 502 c W R I T E ( 6 . » ) KS(1).KS(2).KS(3).KS(4).KT 503 c 504 XE( 1 )=XU(1)+KS(1)*DELTZ 505 X E ( 2 ) = X U ( 2 ) + K S ( 2 ) » D E L T Z 506 X E ( 3 ) = X U ( 3 ) + K S ( 3 ) » 0 E L T Z 507 X E ( 4 ) = X U ( 4 ) + K S ( 4 ) ' 0 £ L T Z 508 X £ ( 5 ) = X U ( 5 ) + K S ( 5 ) « D E L T Z 509 X E ( 6 ) = X U ( 6 ) + K S ( 6 ) « D E L T Z 510 XE(7)=XU(7)+KS(7)'DELTZ 191 I i L i s t i n g of KALMAN4 at 10:07:47 on DEC 13. 1987 for CCld=SEIS 511 X E ( 8 ) « X U ( 8 ) + K S ( 8 ) » D E L T Z 512 X E ( 9 ) - X U ( 9 ) + K S ( 9 ) » D E L T Z 513 C 514 515 RETURN 516 END 517 C ^ 518 C Since the t r a n s i t i o n matr ix , F, Is nonl inear with respect to 519 C the state parameters X2 and X4, i t Is necessary to l i n e a r i z e 520 C It. This Is accomplished In the fo l lowing subrout ine, by 521 C apply ing L e i b n i t z ' s ru le and the extended Kalman f i l t e r 522 C equations. 523 C 524 SUBROUTINE EXTEND(FL,CI,C2.C3,C4,XU.KK.N1,N2.WS.WP.OELT,N3.N4.FT 525 1) 526 C 527 INTEGER MS.SL.SH,N2,N4.FT 528 REAL MM,NN.C1.C2.C3.C4.DELT,N1,N3,KK 529 REAL XU(9).WP(200).WS(200).FL(9.9) 530 C 531 C 532 MM=XU(2)+N1 533 C WRITE(6.*)KK.XU(2).MM 534 IF (KK.LT.XU(2).OR.KK.GT.MM) THEN 535 C 536 FL(8.2)=0.0 537 F L ( 9 . 2 ) - 0 . 0 538 C 539 ELSE 540 C 54 1 SH-N2+1 542 C 543 F L ( 9 , 2 ) = C 1 * ( C 3 » W P ( N 2 ) - C 4 » W P ( S H ) ) * X U ( 1 ) 544 F L ( 8 . 2 ) = C 1 » C 2 * W P ( N 2 ) » X U ( 1 ) 545 C 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) 546 C 547 C 548 END IF 549 C 550 C 551 NN=XU(4)+N3 552 C 553 IF (KK.LT.XU(4).0R.KK.GT.NN) THEN 554 C 555 FL(8.4)=0.0 556 FL(9.4)=0.0 557 C 558 ELSE 559 C 560 SL=N4+1 561 C 562 F L ( 8 . 4 ) = C 1 » C 2 » W S ( N 4 ) * X U ( 3 ) 563 F L ( 9 . 4 ) » C 1 » ( C 3 ' W S ( N 4 ) - C 4 » W S ( S L ) ) » X U ( 3 ) 564 C 565 W R I T E ( 6 . « ) FL(8.4 ) .FL(9 .4) 566 104 END IF 567 C 568 RETURN 192 L is t ing r KALMAN4 at 10:07:47 on DEC 13. 1987 for CCId'SEIS 569 END 570 C 57 1 C This subroutine extrapolates the error covariance matrix 572 C 573 SUBROUTINE EXTRAP(PU.PS.FL) 574 REAL PS(9.9) .PU(9.9) .FL(9.9) .FLT(9.9) 575 REAL TMPI(9.9).TMP2(9.9) 576 INTEGER OR 577 C 578 C 579 OR = 9 580 CALL GTRAN(FL.FLT,OR,OR.OR.OR) 581 C 582 CALL GMULT(FL.PS.TMP1.OR.OR.OR.OR.OR.OR) 583 C 584 CALL GMULT(TMP 1 . FLT . PU, OR . OR . OR, OR .OR. OR) 585 C 586 587 C 588 RETURN 589 END 590 C This subroutine updates the error covir lance matrix P. 591 C 592 SUBROUTINE PUPDAT(HV.KS.PU.PS) 593 REAL KS(9),PU(9,9).HV(9),ID(9.9).TMP1(9.9) 594 REAL TMP2(9.9),PS(9.9) 595 INTEGER OR 596 C 597 OR = 9 598 C 599 CALL UNIT(ID.OR.OR) 600 C 601 DO 1 1=1.9 602 DO 2 d=1.9 603 TMP1(1.J)=KS(I)*HV(J) 604 2 CONTINUE 605 1 CONTINUE 606 C 607 CALL GSUBUD.TMP1.TMP2,OR.OR.OR.OR.OR) 608 C 609 CALL GMULT(TMP2,PU,PS.OR,OR,OR.OR,OR,OR) 610 C 6 1 1 RETURN 612 END 613 C 614 SUBROUTINE MlNSN(LX.X.XMIN.INDEX) 615 C 616 REAL X(2) 617 XMIN=X(1) 6 18 DO 10 I = 1.1.* 619 10 XMIN=AMIN1(XMIN.X(I)) 620 DO 20 J=1.LX 62 1 1NDEX=J 622 IF(X(d)-XMIN) 20.30.20 623 20 CONTINUE 624 30 RETURN 625 END 626 C 193 L i s t i n g of KALMAN4 at 10:07:47 on DEC 13. 198*7 for CCid-SEIS 627 SUBROUTINE MAXSN(LX.X.XMAX.INDEX) 628 C 629 REAL X(2) 630, XMAX-X(l) 631 \ DO 10 1-1.LX 632 10 XMAX-AMAXKXMAX.X(I)) 633 DO 20 J -1 .LX 634 1NDEX-J 635 IF(X{JJ-XMAX) 20.30.20 636 • 20 CONTINUE 637 30 RETURN 638 END 194
Cite
Citation Scheme:
Usage Statistics
Country | Views | Downloads |
---|---|---|
China | 6 | 8 |
Japan | 6 | 0 |
United States | 1 | 0 |
City | Views | Downloads |
---|---|---|
Beijing | 6 | 0 |
Tokyo | 6 | 0 |
Unknown | 1 | 0 |
{[{ mDataHeader[type] }]} | {[{ month[type] }]} | {[{ tData[type] }]} |
Share
Share to: