APPLICATION OF BAYESIAN RECURSIVE ESTIMATION FOR SEISMIC SIGNAL PROCESSING by ERICK BAZIW B . A . S c , The University o f British Columbia, 1986 M . A . S c . , The University of British Columbia, 1988 A THESIS S U B M I T T E D IN P A R T I A L F U L F I L M E N T OF T H E R E Q U I R E M E N T S F O R THE D E G R E E OF DOCTOR OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES (Geophysics) T H E U N I V E R S I T Y OF BRITISH C O L U M B I A A p r i l 2007 © Erick Baziw, 2007 Abstract Bayesian recursive estimation ( B R E ) requires that the posterior density function be estimated so that conditional mean estimates o f desired parameters or states can be obtained. B R E has been referred to as a complete solution to the estimation problem since the posterior density function embodies all available statistical information (i.e., prior, likelihood and evidence). Until recent advances i n B R E , most applications required that the system and measurement equations be linear, and that the process and measurement noise be Gaussian and white. A Kalman filter, K F , (closed form solution to the B R E ) could be applied to systems that met these conditions. Previous applications o f the K F to solve seismic signal processing problems (e.g., deconvolution) have had very limited success and acceptability in the geophysics signal processing community due to the restrictive nature o f the K F . The recently new B R E development o f sequential Monte Carlo ( S M C ) techniques for numerically solving non-stationary and non-linear problems has generated considerable interest and active research within the last decade. This thesis focuses upon the implementation o f S M C techniques (e.g., particle filtering) for solving seismic signal processing problems. A l l the associated filters o f B R E (hidden Markov model filter, K F , particle filter, Rao-Blackwellised particle filter, and jump Markov systems) and a new and highly robust and unique model o f the seismic source wavelet are implemented in two innovative algorithms for solving the important problems o f passive seismic event detection and blind seismic deconvolution. A ground-breaking concept i n blind seismic deconvolution referred to as principle phase decomposition (PPD) is outlined and evaluated i n this thesis. The P P D technique estimates and separates overlapping source wavelets instead o f estimating high bandwidth reflection coefficients. It is shown that one can then easily generate reflection coefficients from the separated source wavelets. In this thesis many advantages o f the P P D are outlined. Simulated seismogram data with low signal-to-noise ratios is blindly deconvolved where non-stationary, mixed-phase, and zero-phase source wavelets are present. I believe that there are currently no existing blind seismic deconvolution techniques which could obtain comparable performance results o f the P P D technique. The work in this thesis has resulted in three I E E E reviewed conference publication. n publications and one peer Table of Contents Abstract ii Table o f Contents iii List o f Tables v List o f Figures vi List o f Abbreviations xii Acknowledgments xiv Dedication xv I. Introduction 1 II. Bayseian Recursive Estimation and Seismic Modeling 4 III. Mathematical background 8 A. Bayesian Recursive Estimation 8 B. Kalman Filter and Jump Markov Linear Gaussian System 9 C. Hidden Markov M o d e l Filter 11 D. Particle Filter and Rao-Blackwellised Particle Filter 13 a) Bayesian Importance Sampling 13 b) Sampling Importance Re-sampling (SIR) 17 c) Choice of Importance Function 20 IV. Seismic Signal Processing 23 A. Passive Seismic Monitoring 23 1) P S - S E E D Filter Outline 25 a) State-Space Formulation 25 b) P S - S E E D R B P F and H M M Filter Implementation 26 c) Evaluating the P S - S E E D with Simulated Data 28 B. Seismic Deconvolution 32 1) Application o f Kalman Filter Smoother to Seismic deconvolution 33 a) Kalman Filter Smoother Governing Equations 34 b) Kalman Filter Seismic Deconvolution ( K F S D ) Algorithm Outline 36 c) A R M A Parameter Estimation by the Least Squares Method 38 d) Evaluating the K F S D Algorithm with Simulated Finite Difference Data 41 2) Application of Particle Filtering to B l i n d Seismic Deconvolution 45 a) Amplitude Modulated Sinusoid 47 b) P P D Algorithm Outlined 49 c) P P D Kalman Filter and Jump Markov Linear Gaussian System Formulation....52 d) P P D Parameter Specification and Algorithm Implementation 53 e) Evaluating the P P D Algorithm with Simulated Data 57 f) Proposed Methodology in Processing Seismograms with the P P D Algorithm...69 g) Estimation of the Dominant Frequency o f the Wavelet to be extracted within the P P D - W E Algorithm 89 h) Utilization o f the P P D - W E Algorithm within Standard Frequency Domain Deconvolution Techniques 104 i) Monte Carlo Simulations Ill j) Sensitivity Analysis 114 3) Implementation o f the P P D - W E Algorithm on Real Data 118 m V. Conclusions 134 References 137 Appendix 143 A. Derivation o f Chapman-Kolomogorov Equation 143 B. Derivation o f the B R E Update Equation 143 C. Derivation o f the Recursive Weights Update Equation 144 D. Derivation o f the "Likelihood" Formula 146 E. P S - S E E D F S M C D and H M M Filter Implementation 148 F. P P D F S M C D and H M M Filter Implementation 150 G. Implementation o f the Three State F S M C D within the P P D - W E Algorithm 152 H. P P D - W E Estimated Source and Residual Wavelets for S C P T Test Hole S C 64.... 154 1) Estimated Source Wavelet at 3m 154 2) Estimated Source Wavelet at 4m 155 3) Estimated Source Wavelet at 5m 156 4) Estimated Source Wavelet at 6m 157 5) Estimated Source Wavelet at 7m 158 6) Estimated Source Wavelet at 8m 159 7) Estimated Source Wavelet at 9m 160 8) Etimated Source Wavelet at 10m 161 9) Estimated Source Wavelet at 11 m: 162 10) Estimated Source Wavelet at 12m 163 11) Estimated Source Wavelet at 13m 164 12) Estimated Source Wavelet at 14m 165 13) Estimated Source Wavelet at 15m 166 14) Estimated Source Wavelet at 16m 167 15) Estimated Source Wavelet at 17m 168 iv List of Tables Table Table Table Table Table Table Table Table 3.1 K F Governing Equations for J M L G S 10 3.2 H M M Filtering Algorithm 12 3.3. Typical Particle Filtering Algorithm 16 4.1 PS-Seed Filter Formulation 27 4.2 K F Fixed-Interval Smoother Governing Equations 35 4.3 P P D Filter Formulation 55 4.4 P P D - W E Filter Formulation 72 4.5 Interval Velocities ( P P D - W E Extracted Source Wavelets) from the Crosscorrelation Technique for S C P T S C 64 128 Table 4.6 Interval Velocities (Processed Seismograms) from the Crosscorrelation Technique for S C P T S C 6 4 128 Table 4.7 Estimated Secondary Arriving Source Wavelet Time Offset for S C P T S C 64 .... 132 v List of Figures Figure 3.1. Schematic illustrating the bootstrap re-sampling methodology Figure 4.1. Source wavelet embedded with varying types o f Gauss-Markov background noise 19 30 Figure 4.2. P S - S E E D output results for estimating state x2 ^ after processing the data shown in Fig. 4.1 31 Figure 4.3. Estimated probability o f an event for the test case shown in Fig. 4.1(b) 32 Figure 4.4. Finite difference simulated source wavelet 41 Figure 4.5. Estimating source wavelet in Fig. 4.4 with 5 order A R M A model 42 Figure 4.6. The reflection coefficients utilized to test the performance o f the K F S D algorithm 422 Figure 4.7. the output after convolving A R M A wavelet with reflection coefficients shown i n Fig. 4.6 43 Figure 4.8. The K F S D estimated reflection coefficients 433 Figure 4.9. Seismic data o f Fig. 4.7 with Gauss-Markov noise o f T = 0.02 ms and s = 0.08 added 44 Figure 4.10. The K F S D estimated reflection coefficients for data shown in Fig. 4.9 45 Figure 4.11. A zero-phase Ricker wavelet with superimposed 50 H z sinusoid with phase o f 96° 48 Figure 4.12. Amplitude modulating term for estimating Ricker wavelet o f Fig. 4.14 48 Figure 4.13. Reconstructed Ricker wavelet superimposed upon actual Ricker wavelet o f Fig. 4.11 49 Figure 4.14. Illustration o f the convolution operation as the summation of A M S source wavelets 50 Figure 4.15. A mixed phased Berlage wavelet 57 Figure 4.16. Reflection coefficients convolved with Berlage wavelet shown in Fig. 4.15 58 Figure 4.17. The output after convolving Berlage wavelet with reflection coefficients illustrated in Fig. 4.16 and adding measurement noise 58 Figure 4.18. Estimating dominant frequency and initial phase component interactively 59 Figure 4.19. The addition o f a 50 H z sinusoid with phase o f 100° at time of 104 ms 60 Figure 4.20. The output o f P P D algorithm where the Berlage wavelets have been separate. .61 Figure 4.21 The error residual which is defined to be the difference between the filtered seismogram o f Fig. 4.17 and the estimated Berlage wavelets o f Fig. 4.20 61 Figure 4.22. The dipole reflection coefficients convolved with Berlage wavelet shown in Fig. 4.15 62 Figure 4.23. The output after convolving Berlage wavelet with reflection coefficients illustrated in Fig. 4.22 and adding Gauss-Markov measurement noise 62 Figure 4.24. The output o f PPD algorithm where the algorithm did an impressive job i n extracting the dipole Berlage source wavelets 63 Figure 4.25. The error residual which is defined to be the difference between the filtered seismogram o f Fig. 4.23 and the estimated Berlage wavelets o f Fig. 4.24 63 th c vi Figure 4.26. The output after superimposing Berlage wavelet with arrival time o f 99 ms and zero-phase Ricker wavelet o f Fig. 4.11 64 Figure 4.27. The output o f P P D algorithm where the algorithm did an impressive job in extracting the Berlage source wavelet from the overlapping zero-phase Ricker wavelet 65 Figure 4.28. The error residual which is defined to be the difference between the filtered seismogram o f Fig. 4.26 and the estimated source wavelets o f Fig. 4.27 65 Figure 4.29. M i x e d phased Berlage wavelet with f= 40 Hz, n = 2, a = 100, and (/) = 10°.. ..66 Figure 4.30. Output after superimposing the time-variant Berlage source wavelets with additive measurement noise (variance = 50 units and time constant = 0.001 ms). .67 Figure 4.31. The addition o f a 40 H z sinusoid with phase o f 100° at time of 104 ms to the sinusoid shown in Fig. 4.18 67 Figure 4.32. The output o f P P D algorithm where the algorithm did an impressive job in extracting the non-stationary Berlage source wavelets 68 Figure 4.33. The error residual which is defined to be the difference between the filtered seismogram o f Fig. 4.30 and the estimated source wavelets o f Fig. 4.32 68 Figure 4.34. The output after convolving Berlage wavelet o f F i g 4.15 with reflection coefficients illustrated in F i g . 4.16 and adding Gauss-Markov measurement noise with a variance of400 units and a time constant o f 0.01 ms 74 Figure 4.35. The output after applying a 100 H z , eighth order, zero-phase shift, low-pass Butterworth frequency filter to the synthetic seismogram shown in Fig. 4.34 75 Figure 4.36. The initial H M M / 3 k phase estimates for first pass o f P P D - W E algorithm on data shown in Fig. 4.35 75 Figure 4.37. The estimated and true Berlage source wavelets for the first reflection coefficient illustrated i n F i g 4.16 77 Figure 4.38. The estimated and true Berlage source wavelets for the second reflection coefficient illustrated in F i g 4.16 77 Figure 4.39. The A M T s for the two estimated A M S s shown in Figs. 4.37 and 4.38 78 Figure 4.40. The normalized cross-correlation function for estimated source wavelets illustrated in Figs. 4.37 and 4.38. Peak value occurs at 4.7 ms 78 Figure 4.41. The estimated and true Reflection coefficients for synthetic seismogram illustrated in Fig. 4.34 79 Figure 4.42. Three state F S M C D realization for time t = 40 ms 79 Figure 4.43. The reflection coefficients convolved with Berlage wavelet shown in Fig. 4.15.80 Figure 4.44. The output after convolving Berlage wavelet o f Fig 4.15 with reflection coefficients illustrated in Fig. 4.43 and adding Gauss-Markov measurement noise with a variance of400 units and a time constant o f 0.01 ms 80 Figure 4.45. The output after applying a 100 H z , eighth order, zero-phase shift, low-pass Butterworth frequency filter to the synthetic seismogram shown in Fig. 4.44 81 Figure 4.46. The initial H M M / 3 k phase estimates for first pass o f P P D - W E algorithm on data shown in Fig. 4.44 82 2 2 2 vn Figure 4.47. Estimated and true Berlage source wavelets for the first reflection coefficient illustrated in F i g 4.44 83 Figure 4.48. Estimated and true time series residual for the last two reflection coefficients illustrated in F i g 4.43 convolved with the Berlage source wavelet o f Fig. 4.15 83 Figure 4.49. The A M T s for the two estimated A M S s shown in Figs. 4.47 and 4.48 84 Figure 4.50. The estimated time series residual of Fig. 4.48 with A M S - E sinusoid with phase of 44° superimposed 84 Figure 4.51. The initial H M M / 3 phase estimates for first pass o f P P D - W E algorithm on data shown in Fig. 4.50 85 Figure 4.52. The estimated and true Berlage source wavelets for the first second coefficient illustrated i n F i g 4.44 86 Figure 4.53. The estimated and true Berlage source wavelets for the third reflection coefficient illustrated in F i g 4.44 86 Figure 4.54. The A M T s for the two estimated A M S s shown in Figs. 4.52 and 4.53 87 Figure 4.55. The normalized cross-correlation function for the first and estimated second arriving Berlage source wavelet shown in Fig. 4.52. Peak value occurs at 4.5 ms. .87 Figure 4.56. The normalized cross-correlation function for the first and estimated third arriving Berlage source wavelet shown in Fig. 4.54. Peak value occurs at 8.9 ms ..88 Figure 4.57. The estimated and true reflection coefficients for synthetic seismogram illustrated in Fig. 4.44 88 Figure 4.58. Illustration o f the filtered synthetic seismogram o f Fig. 4.44 with the A M S - E sinusoid and corresponding initial zero crossing o f 4.2 ms shown 90 Figure 4.59. Illustration o f the output o f the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.44 91 Figure 4.60. Illustration o f B S W 1 with parameters/= 55 Hz, n - 2, a = 170 and k (j) = 60° specified 92 Figure 4.61. Illustration o f B S W 2 with parameters/= 50 Hz, n = 2, a = 170 and <j) = 60° specified 93 Figure 4.62. Illustration o f B S W 3 with parameters/= 45 Hz, n = 2, a = 170 and tj) = 60° specified 93 Figure 4.63. Illustration o f B S W 4 with parameters/= 40 Hz, n = 2, a = 170 and (j) = 60° specified 94 Figure 4.64. Synthetic seismogram generated by summing time variant source wavelets B S W 1 , B S W 2 , B S W 3 , and B S W 4 94 Figure 4.65. The reflection coefficients utilized to generated synthetic seismogram illustrated in Fig. 4.64 95 Figure 4.66. The output after convolving B S W 1 , B S W 2 , B S W 3 , and B S W 4 with the first, second, third, and fourth reflection coefficients, respectively, shown in Fig. 4.65 ..95 Figure 4.67. The output after applying a 100 H z , eighth order, zero-phase shift, low-pass Butterworth frequency filter to the synthetic seismogram shown in Fig. 4.64 96 vin Figure 4.68. Illustration o f the output o f the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.64 97 Figure 4.69. The estimated and true BSW1 source wavelets 97 Figure 4.70. The estimated residual wavelet and the actual residual wavelet (i.e., BSW1+BSW3+BSW4) 98 Figure 4.71. The estimated time series residual o f Fig. 4.70 with the time parameters t* and f set to 11 ms and 8.96 ms, respectively 99 Figure 4.72. Illustration o f the output of the H M M - F E filter after processing the synthetic seismogram shown i n Fig. 4.71 99 Figure 4.73. The estimated and true B S W 2 source wavelets 100 Figure 4.74. The estimated residual wavelet and the actual residual wavelet (i.e., BSW3+BSW4) 100 Figure 4.75. The estimated time series residual o f Fig. 4.74 with the time parameters t* and t' set to 17 ms and 15 ms, respectively 101 Figure 4.76. Illustration o f the output o f the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.75 102 Figure 4.77. The estimated and true B S W 3 source wavelets 103 Figure 4.78. The estimated and true B S W 4 source wavelets 103 Figure 4.79. Typical synthetic seismogram 105 Figure 4.80. The amplitude spectrum o f seismogram illustrated in Fig. 4.79 105 Figure 4.81. Superposition o f filtered seismogram onto seismogram shown in Fig. 4.79 ....106 Figure 4.82. Illustration o f the output of the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.81 107 Figure 4.83. The estimated (solid black line) and true (dotted line) first arriving source wavelet 108 Figure 4.84. The estimated (solid black line) first arriving source wavelet o f Fig. 4.83 with the time series beyond 105 ms set to zero 109 Figure 4.85. The superposition o f the true reflection series (black line) onto the estimated reflection series (dotted line) 110 Figure 4.86. Ten P P D - W E estimates for noisy synthetic seismogram shown in Fig. 4.79. The estimate illustrated by the dotted line differs from the nine other approximations. The bold light grey line identifies the true source wavelet Ill Figure 4.87. Synthetic seismogram where the noise free seismogram o f Fig. 4.79 has additive Gauss-Markov noise with a time constant o f 0.001 ms and a variance o f 0.005 units 112 Figure 4.88. The P P D - W E source wavelet estimates for ten additive Gauss-Markov noise (time constant o f 0.001 ms and a variance o f 0.005 units ) realizations. The bold light grey line identifies the true source wavelet 113 Figure 4.89. The output after overlapping two Berlage source wavelets with additive measurement noise. These two source wavelets have a phase difference o f 170°. 114 Figure 4.90. The output of the P P D algorithm where the algorithm was able to separate the two source wavelets, but has some difficulty in estimating the modulating amplitude term 115 2 2 ix Figure 4.91. The error residual between the estimated P P D Beralge wavelet and the true Berlage wavelet with arrival time o f 99 ms 115 Figure 4.92. The error residual between the estimated P P D Berlage wavelet and the true Berlage wavelet with arrival time o f 108 ms 116 Figure 4.93. Illustration o f the null spaces with a 50 H z (T = 20 ms) source wavelet. Overlapping wavelets which arrive within the null spaces cannot be separated from the source wavelet to be extracted 117 Figure 4.94. Typical S C P T configuration. 120 Figure 4.95. R a w X axis V S P from S C P T hole S C 64 121 Figure 4.96. R a w Y axis V S P from S C P T hole S C 64 122 Figure 4.97. V S P for the processed seismic data captured at test site S C 64 for depths 2 m to 17 m. There is clear evidence o f overlapping source wavelets 123 Figure 4.98. The amplitude spectrum o f the seismogram captured at 2 m 125 Figure 4.99. Estimating t^ in for the seismogram acquired at 12 m 125 Figure 4.100. Extracted (utilizing P P D - W E technique) primary source wavelets V S P for depths 2 m to 17 m o f S C P T hole S C 64 127 Figure 4.101. Residual wavelets V S P for depths 3 m to 17 m o f S C P T hole SC 64 129 Figure 4.102. Estimated reflection coefficients for depth 1 6 m utilizing the W L T 130 Figure 4.103. Estimated reflection coefficients for depth 17 m utilizing the W L T 130 Figure 4.104. Superposition of estimated primary source wavelets for depths 4m to 17m... 131 Figure 4.105. Estimated source wavelet (black line) obtained b y averaging estimated P P D W E source wavelets shown in Fig. 4.104. The light grey lines signify the maximum and minimum bounds based upon the response outlined in Fig. 4.104 131 Figure H . 1. Thirty five P P D - W E primary source wavelet estimations at depth o f 3.0 m 154 Figure H.2. The Averaged Estimated First Arriving Source Wavelet at 3.0m 154 Figure H.3. The estimated first arriving source wavelet superimposed upon residual wavelet at 3.0 m 154 Figure H.4. Thirty five P P D - W E primary source wavelet estimations at depth o f 4.0 m 155 Figure H.5. The Averaged Estimated First Arriving Source Wavelet at 4.0m 155 Figure H.6. The estimated first arriving source wavelet superimposed upon residual wavelet. One can have high confidence in the results due to the fact that the Primary and Secondary source wavelets are very similar in form 155 Figure H.7. Thirty five P P D - W E primary source wavelet estimations at depth o f 5.0 m 156 Figure H.8. The Averaged Estimated First Arriving Source Wavelet at 5.0m 156 Figure H.9. The estimated first arriving source wavelet superimposed upon residual wavelet at5.0m 156 Figure H.10. Thirty five P P D - W E primary source wavelet estimations at depth o f 6.0 m ...157 Figure H . 11. The Averaged Estimated First Arriving Source Wavelet at 6.0m 157 Figure H.12. The estimated first arriving source wavelet superimposed upon residual wavelet at 6.0 m 157 Figure H.13. Thirty five P P D - W E primary source wavelet estimations at depth o f 7.0 m ...158 Figure H.14. The Averaged Estimated First Arriving Source Wavelet at 7.0m 158 x Figure H.15. The estimated first arriving source wavelet superimposed upon residual wavelet at 7.0 m 158 Figure H.16. Thirty five P P D - W E primary source wavelet estimations at depth o f 8.0 m ...159 Figure H.17. The Averaged Estimated First Arriving Source Wavelet at 8.0m 159 Figure H.18. The estimated first arriving source wavelet superimposed upon residual wavelet at 8.0 m 159 Figure H.19. Thirty five P P D - W E primary source wavelet estimations at depth o f 9.0m...160 Figure H.20. The Averaged Estimated First Arriving Source Wavelet at 9.0m 160 Figure H.21. The estimated first arriving source wavelet superimposed upon residual wavelet at 9.0 m 160 Figure H.22. Thirty five P P D - W E primary source wavelet estimations at depth o f 10.0 m .161 Figure H.23. The Averaged Estimated First Arriving Source Wavelet at 10.0m 161 Figure H.24. The estimated first arriving source wavelet superimposed upon residual wavelet at 10.0 m l 161 Figure H.25. Thirty five P P D - W E primary source wavelet estimations at depth o f 11.0 m . 162 Figure H.26. The Averaged Estimated First Arriving Source Wavelet at 11.0m 162 Figure H.27. The estimated first arriving source wavelet superimposed upon residual wavelet at 11.0m 162 Figure H.28. Thirty five P P D - W E primary source wavelet estimations at depth o f 12.0 m .163 Figure H.29. The Averaged Estimated First Arriving Source Wavelet at 12.0m 163 Figure H.30. The estimated first arriving source wavelet superimposed upon residual wavelet at 12.0 m 163 Figure H.31. Thirty five P P D - W E primary source wavelet estimations at depth o f 13.0 m . 164 Figure H.32. The Averaged Estimated First Arriving Source Wavelet at 13.0m 164 Figure H.3 3. The estimated first arriving source wavelet superimposed upon residual wavelet at 13.0 m 164 Figure H.34. Thirty five P P D - W E primary source wavelet estimations at depth o f 14.0 m .165 Figure H.35. The Averaged Estimated First Arriving Source Wavelet at 14.0m 165 Figure H.36. The estimated first arriving source wavelet superimposed upon residual wavelet at 14.0m 165 Figure H.37. Thirty five P P D - W E primary source wavelet estimations at depth o f 15.0 m .166 Figure H.38. The Averaged Estimated First Arriving Source Wavelet at 15.0m 166 Figure H.39. The estimated first arriving source wavelet superimposed upon residual wavelet at 15.0 m 166 Figure H.40. Thirty five P P D - W E primary source wavelet estimations at depth o f 16.0 m . 167 Figure H.41. The Averaged Estimated First Arriving Source Wavelet at 16.0m 167 Figure H.42. The estimated first arriving source wavelet superimposed upon residual wavelet at 16.0 m 167 Figure H.43. Thirty five P P D - W E primary source wavelet estimations at depth o f 17.0 m . 168 Figure H.44. The Averaged Estimated First Arriving Source Wavelet at 17.0m 168 Figure H.45. The estimated first arriving source wavelet superimposed upon residual wavelet at 17.0 m 168 xi List of Abbreviations Abbrv. Meaning A AMS AMT AMS-E ARMA amplitude modulated sinusoid amplitude modulating term A M S source wavelet to be extracted autoregressive moving average B BRE BSW BSD BIS Bayesian recursive estimation Berlage source wavelet blind seismic deconvolution Bayesian importance sampling D DF dominant frequency F FSMCD FMDSM finite state Markov chain distribution Forward Modeling / Downhill Simplex H HMM HMM-FE hidden Markov model H M M frequency estimation J JMLGS jump Markov linear Gaussian systems K KF KFSD Kalman filter Kalman filter seismic deconvolution xn p PF PDF PSM PS-SEED PPD PPD-WE PPD-WEMC particle filters probability distribution function passive seismic monitoring passive seismic signal enhancement and event detection principle phase decomposition P P D filter wavelet extraction P P D filter wavelet extraction Monte Carlo R RBPF Rao-Blackwellised particle filter S S/N SMC SIR STA/LTA SWE SC SCPT signal to noise ratio sequential Monte Carlo sampling importance re-sampling short term averaging / long term averaging source wavelet to be extracted seismic cone seismic cone penetration test V VSP Vp Vs vertical seismic profiling P-wave velocity S-waves velocities W WLT water level technique xiii Acknowledgments The work outlined in this thesis would not have been possible without the guidance o f Dr. Tadeusz Ulrych. Dr. Ulrych's openness to exploring new ideas in seismic signal processing, and his strong background in this area, facilitated the advancement o f m y research interests. Dr. Vikram Krishnamurthy, o f the University o f British Columbia's Electrical and Computer Engineering Department, is an expert in the design and implementation o f Bayesian recursive estimation filters, and I greatly appreciate his agreement to act as one o f m y committee members. Dr. Krishnamurthy is responsible for introducing me to the relatively new and exciting technique o f sequential Monte Carlo filtering or more commonly known as particle filtering. Dr. Michael Bostock's expertise i n seismology was fundamental for a thorough review o f my work. In addition, Dr. Bostock's assistance proved essential i n the completion o f m y thesis and its subsequent submission to the university's Department o f Graduate Studies. xiv Dedication This thesis is dedicated to my beautiful and loving wife, Colleen: without her support, this work would not have been possible. xv Application o f Bayesian Recursive Estimation for S e i s m i c S i g n a l P r o c e s s i n g I. Introduction There is considerable interest in the engineering community in the relatively new implementation o f sequential Monte Carlo ( S M C ) methods [3], [26], [27], [35] to solve optimal estimation problems where the physical system is non-stationary and non-linear. In addition, S M C techniques are ideally suited for real-time optimal estimation. S M C filters are a form o f Bayesian recursive estimation ( B R E ) and they are most commonly referred to as particle filters (PF). In S M C filtering, the required posterior density function is represented by a set o f random samples with associated weights where state estimates are computed based on these samples and weights. To the best of my knowledge, there currently is little published work with regard to the implementation of PF and its variants for solving seismic signal processing problems. Section III o f this thesis outlines the mathematical background o f B R E where the requirement o f modeling the physical problem in a state-space formulation is addressed and the fundamental Chapman-Komolgorov equation is outlined. Additional topics covered include the Kalman filter formulation, jump Markov linear Gaussian systems, hidden Markov model filter, particle filtering, and Rao-Blackwellised particle filtering. Crump [22] and Bayless and Brigham [6] were among the first researchers to fit geophysical problems into a B R E formulation. Mendel and his colleagues [41], [42] carried out work in implementing B R E techniques for solving the seismic deconvolution problem. In Mendel's work, the basic seismogram was modeled as an A R M A process. In the work o f Bayless, Brigham, Crump, and Mendel the Kalman filter (KF) set o f equations were utilized. Kalman filtering is a particular restrictive form o f B R E where certain conditions are required. These special conditions consist o f the case where the measurement and process noise are zero mean independent Gaussian white noise processes, the system equation is a linear function o f the state vector and process noise, the measurement is a linear function o f the state vector and measurement noise, and the initial estimate o f state vector has a Gaussian distribution. The focus o f m y research has been to structure geophysical optimal estimation problems 1 into B R E formulations, so that the relatively new and powerful P F technique and its variants could be utilized. Two specific seismic signal processing problems I address are real-time, non-linear, non-stationary passive seismic event detection and blind seismic deconvolution (BSD). A passive seismic monitoring ( P S M ) system is an assembly o f hardware and software components designed to acquire and analyze, in real time, the acoustic signals collected by an array o f appropriate seismic transducers. Section III(A) outlines a novel, unique and robust algorithm for identifying seismic events within low signal-to-noise ratio (S/N) passive seismic data in real-time. This algorithm is referred as the passive seismic signal enhancement and event detection (PS-SEED) filter. Since the event detection problem is a continuous, real-time process that has non-linear mathematical representations, a Rao- Blackwellised particle filter is utilized. In this algorithm, a jump Markov linear Gaussian system ( J M L G S ) is defined where changes (i.e., jumps) i n the state-space system and measurement equations are due to the occurrences and losses o f events within the measurement noise. The Rao-Blackwellised particle filter ( R B P F ) obtains optimal estimates of the possible seismic events by individually weighting and subsequently summing a bank of Kalman filters. These Kalman filters are specified and updated by samples drawn from a Markov chain distribution that defines the probabilities and transitional probabilities o f the individual dynamical systems which compose the J M L G S . In addition, a hidden Markov model filter is utilized within the R B P F filter formulation so that real-time estimates o f the phase o f the seismic event can be obtained. The P S - S E E D filter is demonstrated to provide up to an 80-fold improvement in the S/N when processing simulated seismic data with highly variable Gauss-Markov measurement noise. Seismic deconvolution is a very widely utilized signal enhancement technique in seismic signal processing. Ideally, by deconvolving the source wavelet from the recorded time series data, only the reflection coefficients remain. The reflection coefficients identify the impedance mismatches between different stratigraphic layers. B R E is ideally suited for seismic deconvolution, due to the fact that one easily accounts for non-stationarity o f the source wavelet and a wide variety o f additive measurement noise models can readily be incorporated into the model. There are two types o f seismic deconvolution problems. The first consist o f the typical case where the source wavelet is assumed known and the 2 reflection coefficients are unknown. The second problem is referred to as blind deconvolution. In this case, both the source wavelet and reflection coefficients are assumed unknown. Section IV(B)-1 outlines a K F algorithm for carrying out seismic deconvolution on simulated data where the source wavelet is known and the seismogram is modeled as an A R M A process (where the reflection coefficients are assumed Gaussian and white). Section IV(B)-2 outlines a innovative and powerful R B P F algorithm for solving the B S D problem. This algorithm is referred to as principle phase decomposition (PPD) and it incorporates a ground-breaking concept when carrying out B S D . The P P D algorithm is shown to have many advantages such as simple filter formulation with minimal parameter specification, conducive to blind seismic deconvolution, assumption o f minimum phase source wavelet is not required, ability to apply frequency filters so that high frequency measurement noise is readily removed, avoids problems associated with band-limited source wavelets when carrying out frequency domain deconvolution, easily handles non-stationary source wavelets, provides for time-variant estimations o f the source wavelet, relatively low associated computer computer-processing cost, reflection coefficients are not required to be represented by discrete state levels, and a whiteness assumption governing the reflection coefficient series is not required. In both the P S - S E E D and P P D formulations, an innovative model o f the source wavelet is outlined and utilized. This source wavelet model is referred to as the amplitude modulated sinusoid ( A M S ) . The A M S representation o f the source wavelet is a highly robust and accurate (with only two defining parameters) model which allows for time variance. The A M S is demonstrated to be a highly accurate representations o f seismic source waveform, mixed-phase wavelets approximation such the for many analytical exponentially decaying cyclic Berlage wavelet, zero-phase Ricker wavelet, and zero-phase Klauder wavelet. In addition, the A M S wavelet has proven very accurate in modeling seismic data acquired during passive seismic monitoring and vertical seismic profiling. In both the P S - S E E D and P P D formulations, all the associated filters of B R E (hidden Markov model filter, K F , particle filter, Rao-Blackwellised particle filter, and jump Markov systems) are implemented. 3 II. Bayseian Recursive Estimation and Seismic Modeling Bayesian recursive estimation is an optimal filtering technique that is based on state-space, time-domain formulations o f physical problems. Application o f this type o f filter requires that the dynamics o f the system and measurement model that relates the noisy measurements to the system state equations be describable i n a mathematical representation and probabilistic form which, with initial conditions, uniquely defines the system behaviour. The potentially non-linear discrete stochastic equation describing the system dynamics is defined as follows: k= x fk-\( k-\> k-\) x p( k\ k-l) u x x (1) In (1), the vector ft is a non-linear function o f the state vector Xk and the process or system noise «£. (1) describes a Markov process o f order one. The sampled, potentially non-linear, measurement equation is given as follows: Zk= k( k> k) h x <-> v P(zk\ k) x (2) In (2), hk depends on the index k, the state Xk, and the measurement noise v* at each sampling time. The probabilistic state-space formulation described by (1) and the requirement for updating the state vector estimate based upon the newly available measurements described by (2) are ideally suited for the Bayesian approach to optimal estimation. The main advantages o f utilizing a state-space formulation in describing physical problems is three-fold. 1) Time variance o f the system and measurement dynamics and statistics can readily be accounted for. 2) Complicated time variant measurement noise models can easily be incorporated into the measurement equations. 3) The ability to utilize process or system noise to compensate for errors i n the mathematical model o f the system dynamics. I f the seismic signal processing problem can be cast into state-space form, then the previously described benefits can be used to obtain an optimal solution. 4 The most important seismic model is, i n general, written as [51] z(t ) = S(t)* ju(t) + v(t) (3) where z(t): is the measured seismogram. S(t) : is the seismic wavelet which is a superposition o f earth and instrument responses. p(t): is the reflectivity o f the earth which consists o f all primary reflections as well as all surface and internal multiples. v(t): is the additive noise, generally taken to be white with a Gaussian pdf. *: denotes the convolution operation. A fundamental task for the seismologist is to estimate the impedance at depth from the recorded seismogram. A commonly adopted simple model, referred to as the Goupillaud layered medium [51], in applied seismology is that o f a horizontally layered, one dimensional earth. The impedance o f the Ath layer of the pancake earth model is defined as (4) k = Pk k s where andF^ v are, respectively, the density and velocity in the M i layer. The relationship between and p^, the reflection coefficient (assuming only primary reflection) o f the kth layer, is k+\~ k s f*k = s k+\ s 5 (5) +k s Rearranging (5) gives k k+\ £ = *i n ~k € 7=1 1-M Therefore, it is required to estimate the reflection series, (6) J J in order to determine s^. In extracting the reflection series, the source wavelet must first be estimated and then deconvolved from the recorded seismogram. If the measurement noise term in (3) is ignored, then the Z transform o f (3) is given as z(t) where ^(z) = S(t)*fi(t)<^Z(z) = S(z) V(z) () x y denotes the Z transform o f the reflection series ju(t). Ulrych and Sacchi [51] eloquently illustrate that the seismogram given in (7) defines a autoregressive moving average ( A R M A ) process i f it is assumed that the reflection series are Gaussian and white. A process is A R i f it is assumed that the present observation is a linear combination o f past observations plus a Gaussian random variable. For stability, the A R process is required to be minimum phase (i.e., zeros outside the unit circle (geophysical convention)). A time series is defined to be a M A process i f it is assumed to be generated by a finite linear combination of past and present inputs only where the input is defined to be white. A s expected, an A R M A time series is generated by the combination o f an A R and M A process. A s outlined in Section IV(B), Mendel [42] fitted the A R M A representation o f the basic seismogram into a linear state-space formulation. Mendel then applied a linear Kalman filter smoother (a particularly restrictive form o f B R E ) to solve the deconvolution problem. A s pointed out in Section IV(B), the main disadvantages o f utilizing a K F S is that it is very difficult to model seismograms by low order A R M A representations when the source wavelet is mixed or zero phase. The A R M A model is highly susceptible to sampling rates where higher rates typically require a higher order for the A R M A model. Most importantly, 6 the K F S algorithm is poorly suited for deconvolving seismograms with non-stationary source wavelets and for addressing the blind deconvolution problem. Another important problem in seismic signal processing is the real-time detection o f acoustic events during passive seismic monitoring. The solution to this problem also requires a non-linear approach. In order to address the issues o f blind seismic deconvolution and passive seismic event detection within the state-space and B R E contexts, the relatively new techniques o f sequential Monte Carlo ( S M C ) filtering are both reviewed and adapted i n this thesis. In general terms, S M C filters are a form o f Bayesian recursive estimation ( B R E ) where a set o f particles define the possible values o f a variable (in the one dimensional case) that one is interested in. These particles are propagated according to the system dynamics, (1), and they are subsequently weighted in a maximum-likelihood sense by using the measurement equations, (2). The set o f particles and their corresponding weights allows one to numerically calculate the posterior density function and subsequently derive a conditional mean estimate o f the variable o f interest. 7 III. Mathematical background A. Bayesian Recursive Estimation In the Bayesian approach to optimal estimation, it is attempted to construct the posterior estimate o f the state given all available measurements [3]. In general terms, it is desired to obtain estimates o f the discretized system equation states x k based on all available measurements up to time k (denoted as zi k) by constructing the posterior p(x : k \zi.-k) and having the initial (prior) pdf o f the state p(xo) specified as an initial condition. The posterior pdf allows one to calculate the conditional mean estimate of the state (E[xk\zi.kJ)B R E is a two-step process consisting o f prediction and update [3]. The prediction step involves using the system equation defined by (1) to obtain the prior pdf o f the state at time k via the Chapman-Kolmogorov equation that is given as follows: P( k x I z i . * - l )=\p( k\ k-\ x I z i . * - i )<bck-\ )P( k-\ x x (8) The Chapman-Kolmogorov equation is derived based upon the transitional densities o f a Markov sequence. The derivation o f the Chapman-Kolmogorov equation is outlined in Appendix A . The prediction step generally deforms, translates, and spreads the state pdf due to the system equations and process noise. The update step computes the posterior pdf from the predicted prior pdf and a newly available measurement. The posterior p d f is updated via Bayes' rule as follows: n/v \r P( k I *i.-jfc) x )-P( k\ k)P(Xk\Z\±-\) Z , (9) X ;—: m ; where the normalizing constant is given as follows: p(z k I zi-k-i) = \P(z k Ix k 8 )p(x k | z\ -\ :k )dx k Equation (4) can also be represented as follows: posterior = likelihood • prior evidence where the prior is given by the system equation, the likelihood is given by the measurement equation, and the evidence is the normalizing constant in the denominator. The update step usually concentrates the state pdf by the action of combining the likelihood of the current measurement with the predicted state [43]. Appendix B outlines the derivation of the BRE update equation and the derivation of the "likelihood" formula is given in Appendix D. The recurrence equations defined by (8) and (9) form the basis for the optimal Bayesian solution. The BRE of the posterior density can be solved optimally (exact solution) when the state-space equations fit into a Kalman filter formulation or a Hidden Markov Model [3]. Otherwise, BRE requires an asymptotically optimal numerical estimation approach when deriving the posterior pdf. B. Kalman Filter and Jump Markov Linear Gaussian System As previously stated, the standard set of KF equations can be implemented as an optimal solution to the BRE when certain conditions are met. These special conditions consist of the case where Uk and v* are zero mean independent, Gaussian white noise processes, fk is a linear function of the state vector and process noise, hk is a linear function of the state vector and measurement noise, and the initial estimate of x has a Gaussian distribution [33]. 0 Similar to the KF, a JMLGS is also defined as a linear Gaussian system, but in this case the system and/or measurement equations (fk and hk) evolve with time according to afinitestate Markov chain [24], [28], [29]. Table 3.1 outlines the KF governing equations for a JMLGS [24]. The index, i, denoted in Table 3.1, facilitates the implementation of a bank of Kalman filters when implementing a RBPF (subsequently outlined). 9 TABLE 3.1 K E GOVERNING EQUATIONS F O R J M L G S Eq. Description Mathematical Representation 10 Finite state Markov chain. 11 System equation. 4~4ii^J = (y kK-i (y kK-i xi 12 Measurement equation. 13 State estimate extrapolation. 14 Error covariance extrapolation. F i +G i k 4=^4>4 4 + k\k-\ = (y k ) 'k-\\k-\ x F l x i\k-i= (y ) U -i (y ) p F i p F i k k T+ k (y k)QUk-i (y ) G i G i T k 15 Measurement extrapolation. 4 = (y k 16 Innovation. 17 Variance o f innovation. 4=4-4 H l k= (y k) \k-x (y k) i si H i pi H i T +R k 18 Kalman gain matrix. 19 State estimate update. 20 Error covariance update. x + K k\k- k\k-l kk xl x x l +A i A' a In (11) and (12) and Uk are mutually uncorrected i.i.d Gaussian zero mean white noise processes with variances o f Qk and Rk, respectively (i.e., v ~ N(0,R ) and Matrices F and G define the k k recursive nature o f the system and process noise, respectively. Matrix F J defines the relationship between the system states and the available measurements. k t u ~N(0,Q )). 10 In terms o f the B R E , the Kalman Filter can be viewed as the following recursive relationship: P( k-\ I Z\.-k-\) = ( k-l x N >' k-l\k-l > k-\\k-\) x x I Z\:k-\ ) = { k: P( k x N P( k\Z\:k) x k\k-\ • k\k-\) x x p = { k: k\k> k\k) N x x ( ) p 21a (21b) (21c) P The K F has the sufficient finite dimensional statistics o f mean and error covariance. C. Hidden Markov Model Filter The H M M filter (also termed a grid-based filter) has a discrete state-space representation and has a finite number o f states. In the H M M filter the posterior pdf is represented by the delta function approximation as follows: p( k-\ x I )= I ^i-M-i^v^it-i - k-l j x ( ) 22 i=\ where and tv k-l\k-l< l = 1 >•••>%> represent the fixed discrete states and associated conditional probabilities, respectively, at time index k-l. The governing equations for the H M M filter are derived b y substituting (22) into the Chapman-Kolmogorov equation (8) and the posterior p d f update equation (9). This substitution results in the H M M prediction and update equations which are outlined in Table 3.2 [3]. 11 T A B L E 3.2 H M M Step 1 2 3 Description FILTERING ALGORITHM Mathematical Representation Initialization (k=0) initialize particle weights. e.g., w[ ~\/N ,i = l, ...,N - Prediction - predict the weights. N ™k\k-\ = Z s S t . s Update-update the weights. I wj^ (z \ xj) Wk\k = -jr, T — ( \ lP Obtain optimal minimum variance estimate o f the state vector and corresponding error k i k\k-A k\ k) 7=1 N . k ~ S ^Arl/t-^/t & w 4 . \ l-l) x z x s x i = 1 p _ ' / . / x ~ L*>k\k< k x k ^ wi kA k x x \T k) x covanance. 5 Let k = k+1 & iterate to step 2. In the above equations it is required that the likelihood pdf p(ik \ k ) x and the transitional probabilities p(x' \ x[_ ) be known and specified. k x 12 D. Particle Filter and Rao-BIackwellised Particle Filter 1) Particle Filter: A s stated previously, the recurrence equations defined by (8) and (9) form the basis for the optimal Bayesian solution, and, except for the ICF and H M M exact solutions, the B R E requires an asymptotically optimal numerical estimation approach. To solve the B R E numerically, a new family o f filters that rely upon sequential Monte Carlo ( S M C ) methods have been made popular within the last decade. This family o f new filters are most commonly referred to as particle filters. Similar to the H M M filter, the particle filter (PF) represents the posterior pdf by the delta function approximation, but in this case a randomized grid is utilized for the estimation o f the posterior pdf. If there was perfect Monte Carlo sampling ( x' ~ p(x | z, ) ) then: k lim N s -> oo E I TV [ k x I Z\:k k :k ] according to the strong law o f large numbers. For non-linear, asymptotically optimal estimation problems, it is nearly impossible to sample i\omp(x k \z . ); x k therefore, the weights in (22) are obtained using Bayesian importance sampling (BIS) [3]. Bayesian Importance Sampling The basic idea in BIS is to represent the posterior density by a set o f samples and corresponding weights utilizing an importance density. Suppose p(x) is the desired probability density function which is to be estimated, where the weighted grid estimate o f p(x) is defined as follows: (23) 13 In (23), it is assumed that the pdf p(x) is related to pdf p*(x) via constant of proportionality Y. A proposal importance density distribution where samples x can easily be obtained l 1 then defined as follows: 1 s q(x) = —I,S(x-x ) N ,=i (24) N l K ) s The numerator and denominator of (23) are multiplied by q(x) to give the following: Y q(x) Y q(x) N s i = x Based upon the results of (25) it is clear that the constant 7 must be given as shown: 2 y.. I'P(X') 1 Substituting (26) into (25) gives the following: ^ N f s [ x - x ' i— i= !j q(x t ^P w (x ) N 1=1 q(x l i ) Comparing (27) with (23) we see that the normalized weights for estimating pdf p(x) are defined as shown: 1 2 Thus q(x) can be estimated numerically via realizations of x'. The normalization requirement for the probability density function. 14 is P = -*L w = { where iv' = ll^l (28) (=1 i=\ q(x ) l If it desired to obtain samples from the posterior densityp(x\z)xp(z\x)p(x), one could assign the prior as the proposal (i.e., q(x) - p(x)) to give the following weight: ~i p(z\x )p(x ) l = (z\x l )p( x ) _ l = ) l P q(x ) x i (29 p(x ) l l For the P F set o f equations, it can be shown as outlined i n Appendix C that the following recursive relationship holds for the importance weights as shown: i k= k-l w :w P(z \x )p(x \x _ ) l k q(x l k l k k x (30) \x _ z ) k k v k The posterior density function is then approximated as follows: I yv\8(x P( k x I Zl.-jfc) = — -x ) l k (31) k 1 z=l A typical particle filtering algorithm is outlined in Table 3.3 [3]. In the algorithm outlined in Table 3.3, the Bayesian importance density is defined to be the prior (p(x' )). k particle degeneracy check is made as indicated in Step 4. 15 In addition, a T A B L E 3.3 TYPICAL PARTICLE FILTERING ALGORITHM Step ~ l 2 Description Mathematical Representation Initialization (k=0) generate sample particles and set corresponding weights. [ Px ~\/N W K Update the weights by the likelihood function and then normalize. Obtain asymptotically optimal minimum variance estimate o f the state vector and corresponding error ~ * [ x 0 ,i s = 1,...,N . S K N . ™\ = k-\P(z J x' ) & w = w / I w\ S W l k l k k N k ~ X Wk k & S x x i = 1 „ i, i . . - \, i « K - J covanance. Sampling Importance Re-sampling (SIR). Re-sample i f N EJF A Uk 7 r eff ~ ~N S <N £ ( i=l T Prediction: take u\ ~ p - N & propagate x l k + l the state particles. Let k = k+1 & iterate to step 2. 16 =Axl ~ W K < ? N ) +Bu\,i = 1,...,N . S b) Sampling Importance Re-sampling (SIR) A common problem with the SIS approach is that after a few iterations, most particles have a negligible weight (the weight is concentrated on a few particles only). This phenomenon is referred to as the "degeneracy problem", and it is due to the fact that the variance o f the importance weights increase over time. A simple statistic to monitor that gives an indication of the degeneracy is the effective sample size N f. A n approximation to N jf is given in Step eJ e 4 o f Table 3.3, where N defines the number o f particles utilized and NT is a user specified s threshold (e.g., N = (0.6 to 0.8) N ). A small value of N jf indicates severe degeneracy. The T s e standard technique to counter the degeneracy problem is to re-sample the particles utilizing a Bayesian bootstrap technique [3] i f the effective number o f particles is less than Nj. The SIR methodology implemented is a Bayesian bootstrap technique. The Bayesian bootstrap technique maps the weighted random measure , w j on to the equally weighted l k random measure jx^, iVJ*j by sampling uniformly with replacement from JJ4 with j \N probabilities jw>£ t . F i g . 3.1 [43] outlines the general idea behind SIR. F i g . 3.1 illustrates t s = 1 the initial particles with uniform weight. After a pass o f the SIS algorithm a degeneracy check is made. If N jf < NT, then a re-sampling type algorithm is implemented. In the e qualitative example shown in F i g 3.1, N particles are generated (with weightNJ^) from the c weighted random measure jx^, w j . N is approximately equal to the number o f uniform k c ~yv jN~ and N > l weights contained within w l k (i.e., N c l k c 1). The newly created and uniformly distributed particles are then propagated through the system equations and the SIR cycle is repeated. 17 A SIR algorithm similar to that described in F i g 3.1 is outlined as follows [3]: SIR Algorithm: • Initialize the C D F : c =w\ • Fori = 2 : N - Construct C D F : c, = c,_, + w' x s k End For • Start at the bottom o f the C D F : i=l • Draw a starting point from the uniform distribution: w, ~ U(0,N] • Forj = l : N ) s - M o v e along the C D F : « =u +N; (j-l) l y l - While Uj > C; - i = i+l - End While - Assign sample: x{ = x' - Assign weight: w{ = J V ; k 1 End For In summary, the SIR technique generates JV children such that ; and var[Nj] = N w[(\ - w[) . s 18 = N withE[N ] = N w[ s t s t If- t t ! ? T O New weights Figure 3.1 Schematic illustrating the bootstrap re-sampling methodology (after, Muhlich [43]). 19 c) Choice of Importance Function The degeneracy problem can be minimized b y selecting an importance density that minimizes the variance o f the importance weights so that N gr is maximized. Consider the eJ following points on importance sampling: • If q(x) =p*(x), (27) reduces to simple Monte Carlo estimation. • Importance sampling is effective when q(x) approximates p(x) over most o f the domain. • There are no absolute requirements for how well q(x) approximates p(x), except that q(x) must not be zero anywhere p(x) is non-zero. In general, all importance functions give correct answers but some require larger TV, values. • It fails when q(x) misses high probability regions of p(x) and systematically yields samples with small weights. • T o overcome this problem, it is critical to obtain data points from importance regions o f p(x) by explicitly searching for significant regions i n the target distribution. The optimal choice o f the importance density that minimizes the variance o f the importance weights is defined to be as follows [3]: q( x |x _, z ) l = (x\\ l k k x k op P x l k l , z ) (32) k Substitution o f (32) into (30) gives the following result: ™k =™k-\P( k z I k-\) x = ™k-liP( k z \ )p(x x / k \x _ )dx' 1 k k x k (33) Equation (33) is the optimal importance density because for a givenx^. ,w takes the same l k value independent from the sample drawn 20 fromq(x k \ x _\, z ) . k k op Therefore, conditional on x _ , l k x the variance o f different ^ r e s u l t i n g from different sampled x'^is zero (i.e., var^hv j = 0). The optimal importance density is advantageous because one can k compute the importance weights n ^ a n d sample - ^ i n parallel, since w ^ i s independent of x . Two major drawbacks of the optimal importance density function are as follows: k • The need to be able to sample from p( x' • The need to compute the integral defined by (33) in closed form. k \ x' ,z kl )• k The two cases i n which the optimal importance can be used are firstly when x k is a member o f a finite set. In this case the integral in (33) becomes a sum and it is subsequently possible to sample from p( x \ x ' , z ) . The second case in which the optimal l k kl k importance can be used is when there are non-linear system dynamics, linear measurement equations, and white Gaussian process and measurement noise [3]. In this case the density defined by (32) is Gaussian [25]. For problems which do not meet the previously defined two conditions, a popular choice o f the importance density is defined as (i.e., the prior) [3],[43]: 9(x \x _ ,z ) l = p(x \x _ ) l k k l l k l k k l (34) Substitution o f (34) into (30) gives: *k = KiP(*k\* ) i k w The prior importance density function is easy to implement, and it tends to move samples towards the region o f high likelihood. T w o drawbacks to the prior importance density function are that it does not take measurements into account and that it is sensitive to outliers. 21 The Rao-Blackwellised particle filter ( R B P F ) allows for the reduction in the number o f particles when implementing B R E on non-linear systems. This is highly advantageous because, even though there is a theoretical independence o f accuracy on the particle dimension, practice has shown that the number o f particles needs to be quite high for high dimensional systems [35]. The R B P F combines a bank o f K F s with a P F . In this case, the K F s are utilized for generating a set o f particles, where the weights o f the particles are calculated with a P F [24]. In essence, the posterior pdf is approximated with a recursive, stochastic mixture o f Gaussians [4], [24], [26]. This type o f particle filtering algorithm is referred to as Rao-Blackwellisation because it is related to the Rao-Blackwell theorem. The Rao-Blackwell theorem is named after Calyampudi Radhakrishna Rao and David Blackwell and it describes a technique that transform a crude estimator into an estimator that is optimal by typically the mean-squared-error criterion [21]. In the R B P F s implemented for the subsequently outlined algorithms, a set o f particles are generated by firstly computing the finite state Markov chain distribution, which is denoted Secondly, based upon the samples drawn from a bank o f Kalman filters (as outlined in Table 3.1) is utilized to compute a set o f particles. The posterior pdf o f the state vector is then calculated and subsequent asymptotically optimal estimates are obtained [24]. 22 IV. Seismic Signal Processing A. Passive Seismic Monitoring There is considerable interest in the engineering community in the real-time identification o f "events" within time series data with low S/N. This is especially true for acoustic emission analysis that is utilized for monitoring and inspecting the integrity and safety o f many structures, such as metal and concrete bridges, gas and oil pipe lines, large storage tanks and aerospace vehicles. A particularly important application o f acoustics emission analysis is within the field o f passive seismology. A passive seismic monitoring (PSM) system is designed to acquire and analyze, in real time, the acoustic signals collected by an array o f appropriate seismic transducers. Seismic activity is often observed in the vicinity o f underground excavations, deep open pits and quarries, around and below large reservoirs where fluids are being injected into, or removed from, permeable subsurface formations, and adjacent to the sites o f large underground explosions [34]. Extreme examples o f energy release can cause violent "rockbursts", resulting in fatalities and injuries among underground personnel and damage to mine structures (e.g., drifts, stopes, shafts, etc.). Passive seismic systems are capable o f detecting rock failures in the vicinity o f underground excavations caused by the sudden release o f strain energy resulting from the redistribution o f stresses around openings. Various hydrocarbon production sites also benefit from seismic monitoring systems during certain phases o f production. Primary or secondary extraction or the injection o f material into the reservoir to enhance production can cause significant stress changes. These stress changes can result in failures o f the overlying strata and the migration o f hydrocarbons to aquifers or to the ground surface. Thus, P S M can be used to satisfy environmental concerns, meet regulatory requirements and assess the development o f induced fracturing within the reservoir. In addition, the passive seismic monitoring systems have been successful in identifying and locating casing failures due to steam stimulation in oil sands [49]. 23 During filling o f hydroelectric or large irrigation reservoirs, changes i n regional loading and pore pressures cause significant stress variations within the surrounding rock mass. These can induce a wide range o f micro and macro-seismic events, some o f which are capable o f causing damage to adjacent structures or to the dam itself. P S M can locate and characterize these potentially hazardous, induced events. Irrespective o f the cause o f seismic events, their reliable detection and placement on a common time base is o f critical importance. This is because the arrival times and phase identification o f the source wavelets (P- and S-wave) at various detector packages, within a three-dimensional array, provide the basis for the calculation o f the location o f the seismic event. Imprecision or uncertainty i n arrival time and phase determination reduces the precision o f the source location operation. For example, in many passive seismic monitoring situations where there is interest i n the behaviour o f specific geological features, or where events must be related to specific structures in a mining or hydrocarbon extraction environment, the accurate determination o f event arrival times is the primary rationale for the installation o f the system. In regions where the level o f induced seismicity is high, and it is accompanied by significant ambient noise, it is essential that the passive seismic monitoring systems possess the capability o f automatically identifying the P- and S-waves generated by seismic events within the noise contaminated seismic time series. Reliable automated identification allows for the timely analysis o f a large volume o f data, and the delivery o f results to the end user in a useful manner. The ability to locate seismic events accurately is directly dependent upon the ability to identify the P-wave, SV-wave and SH-wave responses (phase association) and to determine subsequent arrival times (phase picking) [32]. In addition, the derivation o f important source parameters, such as attenuation, seismic moment, source radius, static stress drop, peak particle velocity, seismic energy and failure mechanism (eg., shear, tensile and casing failure) is directly related to the ability to carry out phase association and identification [50]. Since the monitoring o f seismic acoustic emissions is a continuous, real-time process that typically runs twenty-four hours a day, seven days a week, it requires the incorporation o f real-time, passive seismic signal enhancement and event detection ( P S - S E E D ) digital filters. These digital filters should also provide for non-stationary and non-linear optimal estimation 24 capabilities. To meet these requirements, an algorithm based upon Bayesian recursive estimation ( B R E ) is implemented. The innovative P S - S E E D filter outlined in this section substantially builds upon previous original designs [10], [11] in fitting the P S M event detection filter into a Kalman and particle filter formulation. I also believe that the algorithms outlined in [10], [11] are the first published works that attempt to solve the P S M event detection problem utilizing B R E techniques. In the P S - S E E D algorithm, a novel and unique model o f the seismic event is implemented. This robust source wavelet model is an amplitude modulated sinusoid ( A M S ) . A s outlined in this section, the major improvements to the P S - S E E D algorithm consist o f utilizing a R B P F , which individually weights and subsequently sums a bank o f linear Kalman filters with J M L G S characteristics. The J M L G S describes linear system and measurement equations that can change ("jump") with time. In the P S - S E E D filter, a H M M filter is also applied, which allows for the determination o f time variant phase estimates i f a seismic event is present. 1 ) PS-SEED Filter Outline The P S - S E E D filter is based upon the standard short term averaging / long term averaging ( S T A / L T A ) technique, where an event is declared within the filtered time series, when the S T A / L T A ratio exceeds a predefined threshold [2]. In the P S - S E E D digital filter, the seismic event is modeled as an A M S , which is contained within measurement noise [7], [8], [10], [11], [12]. The rationale for modeling the seismic event as an A M S is due to the difficulty in characterizing time variant, passive seismic source wavelets (e.g., A R M A model [25]) and for added robustness, so that the filter can be applied in varied site conditions with minimal modifications. a) State-Space Formulation A s outlined in [10], [11], [12] the measurement noise is modeled as a Gauss-Markov process with defining parameters o f variance and time constant. This does not preclude any other types o f ambient noise models, but it is required that the measurement noise be specified 25 within a state-space formulation. A s previously stated, the seismic event is modeled as an A M S as follows (continuous form): x\(t) = X2(t)sin[cot + <p(t)] where xj(t) is an approximation to the P-wave or S-wave seismic wavelet (the frequency anomaly), X}(t) is the seismic wavelet's amplitude modulating term, ? is the dominant frequency of the wavelet, a n d / (t) is the corresponding phase. In [10], [11], state X2(t) was approximated as a random walk process, ? was assumed constant, and/(i) was derived i n an ad hoc manner. To allow for greater flexibility, the P S S E E D models state x (t) as a Gauss-Markov process similar to that outlined for the ambient 2 noise model. More sophisticated amplitude models can be implemented, such as the formulation o f a Taylor series on the amplitude dynamics carried out to a third term. The third term (acceleration) is then modeled as a Gauss-Markov process. In the P S - S E E D formulation outlined in this section, it is again assumed that ? is constant (i.e., the investigator is looking for a dominant frequency that represents the P-wave or S-wave), w h i l e / (t) is derived by implementing a H M M filter. b) PS-SEED R B P F and H M M Filter Implementation In general terms, the P S - S E E D obtains estimates o f the possible seismic events by utilizing a R B P F , which individually weights and subsequently sums a bank o f Kalman filters that describe the physics o f the ambient noise and seismic event. These Kalman filters are specified and updated by samples drawn from a Markov chain distribution that defines the probability o f the individual dynamical systems that compose the J M L G S [8], [24]. In addition, a H M M filter is applied, which determines optimal /(t) values i f a seismic event is present. Table 4.1 outlines the P S - S E E D filter formulation where Parameters a(l-2) k b(l-2) k and define the Gauss-Markov processes and A is the sampling rate. Ambient noise parameters ah and bh are adaptively derived from the autocorrelation o f the noise portion of the recorded time series. Parameters a2* and b2 are set b y specifying 1/3 o f the expected k maximum possible variance o f the amplitude o f the seismic event and the corresponding time constant (i.e., correlation between samples) [11]. Appendix E provides for more detail on the implementation o f the P S - S E E D algorithm. 26 T A B L E 4.1 PS-SEED F I L T E R F O R M U L A T I O N Step 1 Description Mathematical Representation Specify and initialize J M L G S system equations and N . System Dynamics a k = x Note: The measurement noise (v ) variance Ru is set to alk (sk) , where (s^ is the variance o f the Gauss-Markov noise. This allows the filter to put more weight on the prior when the GaussMarkov noise becomes more correlated from sample to sample. 0 h-\ s _x2 _ 0 k X\ _ k k-l x k-l a2 x2 k x 2 2 bh-l + h-\ 0 0 u bl _ k x k-i u2 Measurement Equations Case 1 (y\ - ambient noise): k 4 =h +k Case 2 (y2 - ambient noise + event): x v k = x\ + x2 sin(a>kA + <Pk) + k v k k Initialize the prior and transitional pdf for the JMLGS. p(y\),p(yl),&p(y Initialize the prior and transitional pdf for the fixed-grid phase. p( p )&p( p \( J \y _ ), i i J k i ( ( P k ),ij x = l,..,,Np. N = number of fixed-grid values. p 4 Draw samples for finite state Markov chain. 5 Propagate sinusoids based current time index (t =k?) and phase estimates. Sinusoid(t) = sin((Ot+ cp' ) Utilizing (13)-(17) outlined in Table 2.1 and updated sinusoids o f Step 5, propagate the system and measurement equations, calculate importance weights for particles and then update and normalize the weights. Obtain asymptotically optimal estimate o f the state vector. w[=W _,N(z \z ,S )„ Sampling Importance Re-sampling (SIR). Re- y k~ k\yU p i k i k k s • i k N k i = l,...,N . s • w =w /zK k k k~z A k x x i=l sample i f N jr < N e T i=i Use x and H M M filter equations (Table 2.2) to k estimate 10 <p l k if y = l k y2 . k Utilize (18)-(20) to update the bank o f K F s 11 Let k = k+1 & iterate to step 4. 1) de Freitas [23] demonstrates that the importance weights for y\ are given by the predictive density p(z \ z,. _ ,y\. ) = N(z ;z\,S' k k x k k k ), where Ndenotes a Gaussian distribution. 27 c) Evaluating the PS-SEED with Simulated Data This section outlines the performance o f the P S - S E E D for the specific case o f a periodic exponentially decaying source wavelet embedded within Gauss-Markov noise, with widely varying time constants and high variances. The seismic event is modeled with defining parameters o f initial amplitude, damping factor and dominant angular frequency [11]. The simulated seismic event had a frequency o f 200 H z , initial amplitude o f 160 mm/s and 2 damping factor o f 79/s specified. The sampling rate was set at 0.05ms, and a total sampling time o f 300 ms was specified. F i g . 4.1(a) illustrates the source wavelet generated with the previously specified parameters. Figs. 4.2(b)-(g) show the simulated source wavelet with additive Gauss-Markov noise, with variances and time constants specified as (1000 mm /s , 0.0001 ms), (1000 mm /s , 0.0001 2 4 2 4 ms), (4000 mm /s , 0.0001 ms), (1000 mm /s , 0.1 ms), (1000 mm /s , 1.0 ms) and (2000 2 4 2 4 2 4 mm /s , 10 ms), respectively. The source wavelet had an arrival time o f 150 ms (f = 0° i n 2 4 (31)) specified for the simulated traces illustrated in Figs. 4.1(b), (d) and (f). In Figs. 4.1(c), (e) and (g), the source wavelet had arrival times o f 133 ms (f = 140°), 138.7 ms if = 90°) and 164.4 ms (f = 45°) specified, respectively. In Fig. 4.1, the units o f the amplitudes o f the time series data are not displayed in order to reduce clutter, and due to the fact that the S T A / L T A event detection technique is only dependent upon relative amplitudes. The initialization o f the J M L G S system equations (Case 1 and 2 o f Step 1 in Table 4.1) was carried out similarly to that outlined in [11]. The initialization o f the finite state Markov chain probabilities are based upon the likelihood o f an event occurring, and the transitional probability o f moving from a non-event to an event. For the simulations presented here, the 2 1 pdfs and the transitional pdfs were set to p(yo) P(yl I y[-\) = 0-8, p(yl I yl~\) = 0.9, p(yo) = 0- > and p(yl 2 | y\_ ) x 1 1 = 0.1, p(y k \ y _\) k = 0.8, = 0.2, respectively. In general terms, the probability o f an event is about 20 percent. The initialization o f the pdf o f the time variant phase was set to the uniform distribution, while the fixed-grid transitional pdf P(<Pk I Pk-l) w a s s e t high ( -§-> 0-996) and the remaining values were set to have e a uniform distribution. There were 90 (N ) phases specified to reflect the possible phase p 28 shifts o f 0° to 180° in 2° increments. The number of particles (N ) specified for the P S - S E E D S filter was 100. Degeneracy threshold parameter NT was set to 0.8 N . S The P S - S E E D filter was implemented on the simulated data illustrated in Fig. 4.1, with the output for state x2 k (amplitude component of A M S model) illustrated in Figs. 4.2(b)-(g). Fig. 4.2(a) shows the actual modulating amplitude component o f the simulated wavelets. A s is illustrated in Figs. 4.2(b)-(g), there is a very impressive real-time S/N improvement after implementation o f the P S - S E E D algorithm. In these test cases, the S / N is defined as the ratio o f the average maximum amplitude o f the signal to the maximum average amplitude o f the noise. The P S - S E E D was able to increase the S/N by approximately 80-fold, 80-fold, 30fold, 10-fold, 9-fold and 11-fold for the Gauss-Markov additive noise illustrated in Figs. 4.1(b) - (g), respectively. The lower S / N is due to the noise characteristics more closely matching those o f the source wavelet, making real-time noise/signal separation more difficult. The H M M filter portion o f the P S - S E E D algorithm also did a good job o f estimating the phase shifts for the source wavelets illustrated in Fig. 4.1 as indicated by the responses shown i n F i g . 4.2. F i g . 4.3 illustrates the estimated probability o f an event (P{y2 k | z £ )) for the test case shown in F i g . 4.1(b). A s expected, the probability o f an 1 ; event increases along with the seismic amplitude estimate illustrated in Fig. 4.2(b). 29 0 50 100 150 200 250 Time [ms] Figure 4.1. Source wavelet embedded with background noise. 30 varying types of Gauss-Markov 300 0 50 100 150 200 250 Time [ms] Figure 4.2. PS-SEED output results for estimating state x2 k data. 31 after processing Fig. 4.1 300 1 0.8 0.6 OL 0.4 0.2 50 100 200 150 250 300 Time [ms] Figure 4.3. Estimated probability of an event for the test case shown in F i g . 4.1(b). B . Seismic Deconvolution Seismic deconvolution is one o f the most widely researched and implemented seismic signal processing tools. The primary goal o f seismic deconvolution is to remove the characteristics of the source wavelet from the recorded seismic time series, so that one is ideally left with only the reflection coefficients. The reflection coefficients identify and quantify the impedance mismatches between different geological layers that are o f great interest to the geophysicist. In seismology, the recorded time series, z(t), is defined to be the linear convolution o f the source wavelet, S(t), with the earth's reflection coefficients, fi(t), with additive measurement noise, v(t). The mathematically representation o f this relationship is given as shown: z(0 = l M(T)S(t-T)dT + v(t) t o ( 3 ? ) The discrete representation o f (38) is given as follows: k z(k)=Y.MO)S(k-(i-l)) i=l + v(k), k = l,2 — N (38) There are many techniques o f seismic deconvolution that can be implemented so that an optimal estimate is made o f the earth model. A majority o f the standard 32 seismic deconvolution methods utilize the steady state Wiener digital filter that assumes a minimum phase source wavelet [48]. Other techniques implement inverse theory, minimum entropy deconvolution, adaptive deconvolution, complex cepstrum deconvolution and independent component analysis. M a n y o f these deconvolution techniques are ad hoc in nature [42], they are effected by the band-limited nature o f the source wavelet, they are highly susceptible to additive measurement noise, and they assume that the source wavelet is stationary. However, it is readily known that the higher frequencies are attenuated more rapidly than lower frequencies, resulting i n significant variation in the source signal as it travels through the earth. A more effective deconvolution technique would be one that potentially allows for time variation o f the source wavelet and it should be able to handle diverse additive measurement noise [9]. In addition, the case where one attempts to deconvolve an unknown source wavelet from an unknown reflection sequence (blind deconvolution) is also o f great interest. P F is an ideal tool to assess with regards to carrying out blind seismic deconvolution, due to the fact that the solution is non-linear. 1) Application of K a l m a n Filter Smoother to Seismic deconvolution Mendel [42] is one o f the original researchers to carry out work i n fitting the seismic convolution model into a state-space representation for the purpose o f minimum variance seismic deconvolution. In Mendel's approach to minimum variance state-space seismic deconvolution the basic seismogram is represented as an autoregressive moving average process ( A R M A ) [9], [42]. The A R M A model is a combination o f both an autoregressive (AR) process and a moving average ( M A ) process [19]. A n A R time series process is generated by a linear combination o f past observations plus a Gaussian random variable. The M A process is generated by a finite linear combination o f past and present inputs only. In the Kalman filter seismic deconvolution ( K F S D ) algorithm, the recorded seismogram is modeled as an A R M A process, which is driven by a forcing function defined as the in-situ reflection coefficients. 33 a) Kalman Filter Smoother Governing Equations Smoothing is an off-line data processing procedure that uses all measurements between 0 and T to estimate a state at a time t where 0 = t = T [33]. The smoothed estimate of x(t), based on all the measurements between 0 and T, is identified as x(t \ T). There are three types of Kalman Filter smoothers, which are identified as follows: 1. Fixed-interval smoothing: the time interval of measurements (i.e., the data span) is fixed, and we seek optimal estimates at some, or perhaps all, interior points. 2. Fixed-point smoothing: an estimate at a single fixed point in time is obtained, and the data span time, T, is assumed to increase indefinitely. 4. Fixed-lag smoothing: it is again assumed that T increases indefinitely, but in this case we are interested in an optimal estimate of the state at a fixed length of time in the past (i.e.,x(t-5\T), withdheld fixed). In the KFSD algorithm, a Fixed-Interval Smoother is applied. The discrete optimal fixedinterval smoothed estimate x(k \ N) (where N = T/A and A is the sampling rate) is defined in Table 4.2 [33]. The computational sequence for the discretefixed-intervalKalman Filter smoother can be thought of as a two pass process. In the first pass, optimal real-time state estimates are obtained by implementation of equations (9)-(20) outlined in Table 3.1. In the second pass, results from thefirst-passestimates (i.e., x and P ) are reprocessed starting from time index N and utilizing equations (39) to (43). Gelb [33] defines a state to be smoothable (Smoothability Condition) if an optimal smoother provides a state estimate with lower error covariance values compared to the real-timefilteredestimate. Only those states, such as randomly time-varying states, which are controllable by the noise driving the system state vector are smoothable. Therefore, constant states are not smoothable. 34 T A B L E 4.2 K F FIXED-INTERVAL S M O O T H E R GOVERNING Eq. Description 39 State smoother estimate. Backward recursive equation. 40 41 Smoothing state transition matrix. 42 Smoothing error covariance matrix 43 Covariance matrix of r . lW In In In In (39) (40) (42) (43) EQUATIONS Mathematical Representation xtm = k\k-l k\N = x + Pi k\k-\ k\N x +r "JW J\ = ~ j r r 1S p p = k + + Hj [HjPjj MUf ^ [ F r . 7 K H k k R j \ l [ Zj - H _} jXj{j \ l P \N = k\k-\ - k\k-\ k\N k\k-\ p p k J\N = S F J + i j+\\N j i S p + s p + Hj k = N - l , N - 2 , •••, 1 and n x 1 vector r is a n x 1 vector. j = N , N - l , •••, 1 and r(N+l \N) = 0. k = N - l , N-2, , 1. j = N , N - l , •••, 1 andS =0. N+1 + {N 35 1-1 H ) J\}-\ 'j P H + j R H x b) Kalman Filter Seismic Deconvolution (KFSD) Algorithm Outline The first step i n the K F S D algorithm is for the user to specify the order o f the A R M A process. For example, the Z transform for a fourth order A R M A model is given as follows: bz +bz A V(z) = z 4 +bz 3 x +aiz __X(z) l 3 +a z 3 +bz 2 2 4 + a z^ + a 2 2 3 (44) U(z) 4 where the parameters, bi, bj, b$, and b , define the M A process coefficients, while the 4 parameters, ai, a , 05, and a define the A R process coefficients. V(z) is the Z transform o f 2 4 the source wavelet, X(z) is the Z transform o f the seismogram, and U(z) is the Z transform o f the in-situ reflection coefficients. In (44), the output Xk+ is estimated on the basis of x +3, x +2, x +i, x , p-k+4, Mk+3, jUk+2, and fik+i 4 K k K k according to the following equation: k + 4 + l k+3 x + 2 k+2 a x a + 3 k +l + 4 k = x lMk+4 a x a (45) x + 2Mk+3 + 3Mk + 2 + 4^k + \ b b b b For the 4 order case , the K F S D state-space formulation is derived b y introducing variable t h 3 Yl(z) into (44) [42] outlined as follows: bz 4 x V(z) = z 4 +b z +b z 3 2 2 +aiz 3 +a z 2 3 2 +b z Y\(z) l 4 +az X 3 +a X) Y z 4 _X(z) (46) U(z) Equating numerator and denominator terms i n (46) results i n the following two expressions: k = \y^k+4 x 3 b + 2y^k+3 +hy k+2 b l C a n be generalized to an A R M A model o f any order. 36 + 4y^k+\ b (47a) Mk =y k+4 + \y k+3 l a + 2ylk+2 1 + ^y k+l a + a y\ l 4 k B y choosing xl = yl , x2 = yh+i, x3 = yl +2, and x4 = yl +3, (47a) and (47b) can be written as into a state-space formulation as follows: k k k k k+l xl x2 k+l x3 /t+l k k k 0 1 0 0 0 0 1 0 x2 0 0 0 1 x 3 ~3 ~2 k+l x4 a k -ai a "0" k x l + k 0 *k_ 1 x where the impedance structure /i is defined as E\ju ju \ = Q d(k-r). k k T 0 k (48) ju is a Gaussian k white noise process with mean zero and variance Q . k The discrete measurement equation is given as follows: xl k k z x2 k = [b b A h] 2 (49) x3 k *k x The computational sequence for the discrete K F S D algorithm is outlined as follows: 1. Specify the order o f the A R M A process and derive coefficients (e.g., ay, a , a , a , bi, b , b , and 64) (subsequently addressed). 2 3 4 2 3 2. Define the state-space matrix equation (48) and measurement equation (49). 4. Obtain the filtered estimates (e.g., x \ , P \ , K ) b y implementation o f equations (9)k k k k k (20) outlined in Table 3.1. 4. Obtain smoothed estimates (e.g., x \ , k N r ^, P \jq, F^ ) by utilizing the values derived i n p k k step 3 and subsequently implementing equations (39) to (43) outlined in Table 4.2. 5. Derive reflection coefficients or in-situ impedance structure by implementing [42] 37 Mk\N (50) =Qk k k+\\N G r A s previously stated, the first step i n the K F S D algorithm is for the user to determine the order o f the A R M A process and subsequently derive the necessary model parameters. This portion o f the K F S D analysis is referred to as system identification. The maximum- likelihood o f approximating the true source wavelet with an A R M A model increases monotonically with increasing system order, while the computational cost o f increasing the system order is proportional to n where n is the A R M A model order [42]. In general terms, 3 the process o f determining the A R M A model order is a trial and error approach. In this analysis, the investigator chooses a model that has the smallest number o f parameters while meeting a performance index requirement that measures how well the A R M A model fits the actual in-situ model. The technique that I utilized to derive the A R M A model parameters is based upon the work o f Ogata [44]. In this approach to system identification, a least squares cost function that is defined as the difference between the A R M A model response, and the corresponding experimental response, is minimized. c) A R M A Parameter Estimation by the Least Squares Method The derivation o f the A R M A model parameters b y utilizing a least squares method is demonstrated by again considering the 4 order (i.e., n = 4) A R M A model given in (44). If th the numerator and denominator o f (44) are multiplied by z , we obtain the following: 4 X(z) U( ) z _ b +b z l \ + a\Z +b z 1 2 1 +b z 2 3 +a z 2 2 + a 3 3 z In (51), the output x is estimated on the basis o f x^-; x . k k according to the following equation: 38 (51) 4 2> + « 4 Z -4 x s, x . , ju , p, .i, fi -2, and /u ,3 k k 4 k k k k *k = - \ a X k - l - 2 k-2 a X ~ i a X k - 3 ~ a 4 k - 4 X V * + 2Mk-l b + + hMk-2 + b Mk-3 4 (52) where x is the estimated value of x . k k The error between the estimated output value and actual output value is defined as follows: = k - k k x e (53) x Since x depends on past data up to n sampling periods earlier, the error e is defined only for k k k=n. B y substituting (52) and k = n, n+1,... Ninto (53) and combining the resultingN-n+1 equations into a vector-matrix equation, we obtain the following: x where x N = [x x 4 5 ... x ], N = C N q + e N (54) N q = f-aj, -a , -a , -a , b;, b , b , b ], e = [e e ...e J, and CN is N N 2 3 4 2 3 4 N 4 5 N defined as x x C x 2 x 3 x 0 x x 2 3 4 l x l M4 M3 M2 Ml M5 ju M3 M2 4 -1 N X/V-l N-2 x N-3 N-4 x MN x MN-l MN-2 MN-3 The least squares performance index is defined as follows: J = ~ X s l X N k = 4 2 k =\s 2 / N s N (56) The least squares method involves minimizing (56), such that the A R M A parameter values w i l l best fit the observed data. In Ogata's [36] formulation, it is assumed that the input sequence {p. j is such that for N>4, k c' C is nonsingular. Ogata shows that the optimal N N 39 estimate o f q is given as shown: N IN N C C c yN N \ (57) N In (57) it is required that {ju J is sufficiently time-varying, so that d C k N is non-singular. N Equation (57) is a first best estimate (in a least squares sense) o f the A R M A parameters. Ogata presents a subsequent recursive formulation for the estimate o f A R M A parameters, utilizing (57) as an initial estimate. The recursive least square estimation is defined as follows: 9 N+l ~9N LvjV+1 + N+l K (58) - N+1 N] C x where CC\ N K N+l c'N+l N (59) 1 + c N+l and N+\ c = LViV : yN-l ' yN-2 : '• yN-3 ' MN+l ' MN '• MN-1 : : ' f^N-2 ] : (60) The method that I utilize for determining the A R M A parameters for the source wavelet is to firstly convolve the isolated wavelet with a highly variable and known noise process with mean zero and unit variance. This insures that {ju } is sufficiently time-varying so that k CC f N N is non-singular. Initial estimates o f the A R M A parameters are obtained by implementation o f (57), known /i , and the convolved output sequence y . The initial k estimates o f q N k are then fed into the recursive estimation equation, defined by (58), until the performance index (56) reaches a predefined minimum. 40 d) Evaluating the K F S D A l g o r i t h m with Simulated Finite Difference Data This section presents test test-bed simulation results when implementing the KFSD algorithm previously outlined. The first step in the simulation was to define a seismic source wavelet. A m i n i and Howie utilized the finite difference program F L A C to model seismic source wavelets [9], [36]. F i g . 4.4 illustrates the simulated source wavelet generated by A m i n i and Howie and obtained by personal communication [36]. The wavelet shown in Fig. 4.4 was generated by assuming a uniform half-space, with an in-situ shear wave velocity o f 180 m/s, and a sampling rate o f 0.02 ms. I' •i• 01 • i • • i ••I 2 3 4 5 6 7 8 | i i| ) M ,) | i , , i .,)n, J | i •j, i | . i | • • , i . , n , M | u 9 10 11 12 13 14 15 16 17 18 19 20 21 22 Time (ms) Figure 4.4. T h e finite difference simulated source wavelet. The K F S D algorithm was then implemented on the data illustrated in F i g . 4.4. In deriving the necessary A R M A parameters, a 5 th order A R and M A was utilized with a data compression ratio o f 16 to 1, so that the sampling rate became 0.32 ms. In addition, it was found that the A R M A model estimation algorithm worked best when the source wavelet was time reversed. The only impact that time reversing the source wavelet has on the K F S D algorithm is that the recorded seismic time series must be time reversed when processing. Fig. 4.5 illustrates the results o f the A R M A estimation o f the time reversed source wavelet shown in Fig. 4.4 41 1 The estimated A R M A model is then convolved with the impedance structure illustrated in Fig. 4.6 to give the output shown in Fig. 4.7. If the K F S D algorithm is applied to the seismic data shown in Fig. 4.7, the reflection coefficients illustrated in F i g . 4.6 are recovered exactly. 0.8-j 0.6 | : 0.4•g 0.2^ I °: I -0.2: -0.4-0.6-0.84 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 Time (ms) Figure 4.6. T h e reflection coefficients utilized to test the performance of the K F S D algorithm. 42 300 Figure 4.7. The output after convolving the A R M A wavelet with the reflection coefficients shown in Fig. 4.6. Gauss-Markov noise is then added to the synthetic seismogram shown in Fig. 4.7, where a time constant o f T = 0.02 ms and standard deviation o f s = 0.002 was specified. The K F S D c algorithm derived the output shown in F i g . 4.8. A s is illustrated in Fig. 4.8, the arrival time locations o f the reflection coefficients is recovered exactly, but there is some degradation in the estimation o f the amplitudes. The K F S D algorithm gives accurate estimates o f the relative amplitudes o f the reflection coefficients. 0.8 0.60.4at 0.2- Q. E < -0.2-0.4- -0.6 • 0 i • • i 20 40 1 1 i 60 • 1 i 80 • • i 100 ' • i 120 i 140 160 . i 180 Time (ms) Figure 4.8. The K F S D estimated reflection coefficients. 43 i , , . 200 , | 220 , . , . 240 . | 260 . i , i 280 i 300 The K F S D algorithm was next tested for its ability to derive the reflection coefficients within a high noise environment. F i g . 4.9 shows the seismic data o f Fig. 4.7 with GaussMarkov noise o f T = 0.02 ms and s =0.08 added. The K F S D algorithm derived the c estimates illustrated in Fig. 4.10. A s is shown in Fig. 4.10, the arrival time locations o f the reflection coefficients are recovered exactly, but the algorithm fails to determine the amplitudes. The K F S D algorithm gives accurate estimates o f the relative amplitudes o f the reflection coefficients which facilitates the investigator to recover the true amplitudes. Mendel [42] outlines a fixed—point smoother which can be utilized to estimate the reflection coefficient amplitudes once the corresponding time locations are known. i . . . . . . . . . . . 0 20 40 i . . . . . . . . . . . 60 80 100 i . . . 120 i . . . i . . . 140 160 Time (ms) i . . . 180 i . . . . . . . 200 220 i > , . , . . . 240 260 i . . . i 280 300 Figure 4.9. Seismic data of Fig. 4.7 with Gauss-Markov noise of T = 0.02 ms and s = C 0.08 added. 44 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 Time (ms) Figure 4.10. The K F S D estimated reflection coefficients for data shown in Fig. 4.9. Fig. 4.11 shows the K F S D output for the data derived when the true time reversed wavelet of Fig. 4.4 is convolved with the reflection coefficients in Fig. 4.6 and Gauss-Markov noise of T = 0.02 ms and s =0.08 is added. The K F S D algorithm provided very similar results to c those illustrated in Fig. 4.10. The main disadvantage o f utilizing the previously described K F S D algorithm is that it is very difficult to model seismograms as A R M A models when the source wavelet is mixed or zero phased [37]. In addition, the A R M A model is highly susceptible to sampling rates where higher sampling rates typically require a higher order for the A R M A model, the K F S D makes the assumption that the reflection coefficients are Gaussian and white, and the K F S D algorithm is poorly suited for deconvolving seismograms with non-stationary source wavelets and for addressing the blind deconvolution problem 2) Application of Particle Filtering to Blind Seismic Deconvolution L i u and Chen [40] were one o f the first researchers to address the blind deconvolution problem utilizing S M C techniques. Their work predominantly focused upon deconvolving signals within the digital communications field. In this case, the source wavelet, (S(k)) in 45 (38) is referred to as the blurring mechanism, and the reflection coefficients are the unobserved input signal, which is assumed to have discrete values with known state levels l'"'' m s (similar to a finite state Markov chain). Fortunately for digital communication s problems, the blurring function (i.e., source wavelet) has very few coefficients and there are typically a very limited number o f state levels describing the input signal. For example, i n the simulated results given b y L i u and Chen [40], the number o f blurring coefficients was three, with 16 possible state levels for the input signal. This system configuration results in a total o f 16 = 4,096 possible combinations. Unfortunately, a typical seismic source wavelet 3 has a minimum o f 30 coefficients, and the reflections coefficients would require a minimum of 40 state levels . This system configuration results in a total o f 4 0 4 30 = 10 48 combinations. In addition, the input data communication signals are characteristically highly variable, unlike the sparse reflection coefficients typically encountered in seismology. The number o f combinations is indicative o f the required particles for solving the blind deconvolution problem, and 1 0 48 is not a feasible number. For these reason, a new approach in carrying out blind seismic deconvolution, which relies upon the unique A M S model o f the source wavelet, is outlined in this section. Similar to the P S - S E E D algorithm, the source wavelet is modeled as an amplitude modulated sinusoid and the blind deconvolution is carried out by initially determining the seismogram's principle phases. Once the principle phases are determined, a R B P F algorithm is utilized to separate the corresponding overlapping source wavelets. The subsequently outlined techniques are referred to as the principle phase decomposition (PPD) and principle phase decomposition wavelet extraction ( P P D - W E ) algorithms. The P P D and P P D - W E techniques makes use o f the fact that the discrete convolution operation outlined by (38) can also be represented as the summation o f several source wavelets o f differing arrival times. 4 Range o f values of+0.2 to -0.2 in steps o f 0.01. 46 a) Amplitude Modulated Sinusoid Similar to the P S - S E E D algorithm, a fundamental and unique component o f the P P D and P P D - W E algorithms is that the seismic source wavelet is modeled as an amplitude modulated sinusoid. The A M S representation o f the source wavelet has proven to be highly robust and accurate i n the case of passive event detection, where it is desired to identify nonstationary source wavelets within complicated measurement noise in real-time [10], [11]. In addition, a source wavelet generated from vertical seismic profiling can be modeled as an A M S [9]. The A M S has also proven to be a good approximation for many analytical representations o f seismic source wavelets, such as the exponentially decaying cyclic waveform [47], mixed-phase Berlage wavelet [1], zero-phase Klauder wavelet [38], and zero-phase Ricker wavelet [47]. The following example demonstrates the robustness o f the A M S model. The Ricker wavelet is mathematical represented in the time-domain as [47]: /^^^(l-S^viO-^o) )^-^^ ^, 2 0 <>t 0 (61) where Ao = wavelet maximum amplitude (centered between two flanking lobes), VM = dominant or peak frequency o f the Ricker wavelet, and t = wavelet arrival time o f 0 maximum amplitude. Fig. 4.11 illustrates a zero-phase Ricker wavelet, with the following parameters specified, Ao = 160 units, VM = 50 H z , and t = 120. Superimposed upon the Ricker wavelet is a scaled 50 0 H z sinusoid with a phase o f 96°. I f the sinusoid shown i n Fig. 4.11 (unsealed) is multiplied by the P S - S E E D numerically estimated amplitude modulating term shown i n Fig. 4.12, then the source wavelet illustrated i n F i g . 4.13 is obtained (superimposed upon the true Ricker wavelet o f Fig. 4.11). Although the Ricker wavelet has a peak frequency, it doesn't have a specific sinusoidal term, and one is still able to reconstruct the desired wavelet b y applying an appropriate amplitude modulating term. 47 60 70 80 90 100 110 120 130 140 150 160 170 180 Time (ms) Figure 4.11. A zero-phase Ricker wavelet with superimposed 50 Hz sinusoid with a phase of 96°. 80 90 100 110 120 130 140 150 160 170 Time (ms) Figure 4.12. Amplitude modulating term for estimating the Ricker wavelet of Fig. 4.11. 48 180 Time (ms) Figure 4.13. Reconstructed Ricker wavelet superimposed upon actual Ricker wavelet of Fig. 4.11. b) P P D Algorithm Outlined. The formulation o f the P P D algorithm has similarities to the passive seismic signal enhancement and event detection ( P S - S E E D ) fdter previously outlined, with the added complexity o f modeling and reconstructing several overlapping source wavelets, as opposed to identifying a single source wavelet ("event") i n real-time. In the P P D algorithm, the only prior requirement is that the source wavelet be modeled as an A M S . A s opposed to standard seismic deconvolution techniques that attempt to derive reflection coefficients, the P P D algorithm decomposes the seismogram into its associated overlapping source wavelets. This mitigates the problem associated with obtaining high bandwidth reflection coefficients from band-limited source wavelets. If desired, reflection coefficients can be readily constructed from the derived overlapping source wavelets. The P P D algorithm also provides for timevariant estimations o f the source wavelet. This information could be utilized within nonblind (known source wavelet) deconvolution techniques. The concept o f decomposing a seismogram into overlapping A M S source wavelets, based upon inherent phase differences, is illustrated i n Fig. 4.14. F i g . 4.14 shows three sinusoidal wavelets that have an identical dominant frequency and varying phase components o f 49 0°, 110°, and 250°. Three source wavelets are generated by modulating the amplitudes o f the sinusoidal wavelets (outlined i n bold), as is shown in F i g 4.14. The seismogram can be generated and decomposed b y superimposing the A M S source wavelets and separating the A M S wavelets, respectively. A Y A A A A /A AA \/A A v/ A ^ / u f VA \ n y /7 SAA<\ A l A A A7 \ / V'/- \ / V Y: AA Figure 4.14. A n illustration of the convolution operation as the summation of A M S source wavelets. Similar to the P S - S E E D filter, the P P D algorithm implements a R B K F that individually weights and subsequently sums a bank o f linear Kalman filters with J M L G S . The K F s define the system and measurement dynamics for event and non-event conditions. The event condition is associated with the case o f a source wavelet or overlapping source wavelets being present within the recorded time series at time index, k. The non-event condition represents the case when only measurement noise is present. A s shown in Table 3.1, the Kalman filters are specified and updated by samples drawn from a finite state Markov chain distribution ( F S M C D ) . The F S M C D defines the probability and transitional probabilities o f the event and non-event conditions. In the P P D algorithm, there is the added complexity that more than one source wavelet can be present at time index, k. The number o f Kalman filters required and the associated definition o f the system and measurement 50 equations must reflect that possibility o f overlapping source wavelets. In shallow and deep reflection seismology, it is reasonable to assume there w i l l be a maximum o f three to four overlapping source wavelets recorded within the seismogram. For example, Sheriff and Geldart [47] demonstrate that the expected vertical resolution o f a noise free reflection seismic survey is ?/4, where ? defines the wavelength o f the source wavelet. For shallow investigations (20 m to 100 m), it is feasible to assume a dominant compression source wave frequency (/) o f approximately 100 H z (T = l/f= 10 ms) and a corresponding velocity (V ) o f 1500 m/s. Based upon these parameters, p the calculated vertical resolution is 4.75 m (note: X = V/f). The associated two way travel time (At) o f the source wavelet for a vertical resolution o f 4.75 m is 5 ms (i.e., At = 2*4.75/V). If the period (10 ms) o f the source wavelet is divided by the two way travel time, a ratio o f two is obtained. This implies that for this particular configuration, one would not expect more than three overlapping source wavelets. This assumes that the duration o f the source wavelet is equal to its period. A similar calculation can be made for a deep reflection (500 m to 5 K m ) survey, where V p « 3000 m/s to 5000 m/s and / ~ 50 H z . If it is assumed that there is a maximum o f three overlapping source wavelets for an event condition, then there are a total o f seven overlapping source wavelet combinations. These combinations are outlined as follows: • Only one source wavelet present resulting in three possible configurations ((1,0,0), (0,1,0), and (0,0,1)); • Two overlapping source wavelets present, resulting in three possible configurations (1,1,0), (1,0,1), and (0,1,1); • Three overlapping source wavelets present (1,1,1). If one assigns a total o f 200 event/no-event K F s (particles) for each time increment, and it is assumed that there is a 40 percent chance o f an event, then the total number o f required K f s (particles) is 200*0.4*7 (event) + 200*0.6 (no-event) = 560+120 = 680. This is a dramatically lower number than the previously outlined value o f 1 0 requirement o f modeling the reflection coefficients as discrete state levels. 51 48 without the In addition to utilizing a R B P F , and similar to the P S - S E E D algorithm, the P P D algorithm implements a set o f H M M filters in order to refine the initially specified phase components o f the A M S source wavelets. For example, i f an approximate phase o f 100° is specified with a fixed-grid phase interval resolution o f ±15°, then the H M M filter w i l l adjust the 100° phase within the -85° to 115° phase interval window to give an optimal estimate. c) P P D K a l m a n Filter and J u m p M a r k o v L i n e a r Gaussian System Formulation The P P D algorithm models the amplitude modulating term of the A M S source wavelet as a slowly varying (high time constant) Gauss-Markov process. When specifying the K F system equation, the user must take into account the maximum number o f possible overlapping source wavelets. For example, (62) illustrates the required K F system equation for the case o f a maximum o f three overlapping source wavelets. The state vector represents the amplitude component o f the three source wavelets and the additive background noise (state x4k) that is modeled as a Gauss-Markov process with defining parameters o f variance and time constant. A s previously stated for the P S - S E E D algorithm, this does not preclude any other types o f ambient-noise models, but it is required that the measurement noise be specified within a state-space formulation. The P P D measurement equation incorporates the J M L G S dynamics o f the S B D problem, where the measurement equation can change ("jump") based upon the event and no-event conditions. A n event condition signifies that source wavelets are present. The P P D algorithm must take into account the possibility that only one source wavelet or several overlapping source wavelets are present when an event condition occurs. Equation (63) outlines the P P D measurement equation for the B S D problem when it is assumed that a maximum o f three source wavelets are overlapped at time index, k, for the event condition. 52 \ 0 0 0 xl 0 a,2 0 0 x2 0 x3 + 0 a k+\ x2 0 0 k+\_ 0 0 ^k+i x x4 k z = a 0 ^k sin{cokA + (p\ ) + x2 x k 3 k 0 a_ n 0 k k k k_ x4 h 0 0 sin(a>kA + (p2^ ) + x3 0 k 0 0 " u\ 0 0 u2 h 0 0 u3 k k (62) k b _ _u4 n k sinicokA + (p3 ) + x4 k k (63) In (63), co represents the dominant frequency o f the source wavelet and f (1..3)k are the H M M filter phase estimates at time index, k, for source wavelets 1 to 3, respectively. A s previously outlined, for each event condition there are seven possible combinations ((1,0,0), (0,1,0), (0,0,1), (1,1,0), (1,0,1), (0,1,1) and (1,1,1)) for the case o f three source wavelets. The P P D algorithm also synthesizes measurements for the no-event condition. In these cases, synthetic measurements o f very low value are assigned to the respective source wavelets. This removes the residual influence o f source wavelets that do not contribute to the estimated measurement equation described by (68). d) PPD Parameter Specification and Algorithm Implementation The first step i n formulating the P P D algorithm consists o f specifying the maximum number of possible overlapping source wavelets and outlining the corresponding K F system and measurement equations. Subsequent to this, the dominant frequency and approximate values of the A M S source wavelets' initial phase estimates and corresponding fixed-grid phase interval resolutions are specified. This information is readily available from the seismogram under investigation. Section IV.B.2.g outlines an algorithm utilized for estimating the source wavelet's dominant frequency. After the dominant frequency and approximate phase of the first source wavelet is estimated, approximate values o f the remaining A M S source wavelet phases and corresponding phase resolutions are specified. This can be carried out 53 utilizing a parameter estimation program, such as the simplex method. It is also possible to implement the P P D filter without the requirement o f specifying all the initial phase estimates. This P P D algorithm configuration is outlined in Section IV.B.2.f. In this case, the overlapping source wavelets are extracted sequentially and chronologically from the seismogram. The only phase which is required to be specified is that o f the current source wavelet to be extracted. The final important parameter to specify within the P P D algorithm is the number o f Kalman filters or particles required. A s previously outlined, in the R B P F implemented for the P P D algorithm, a set o f particles (N ) ; s that define the event/no-event 5 condition are initially generated by computing a F S M C D . If an event has occurred, then a set o f additional K F s or particles are generated, equal to the total number o f overlapping source wavelet combinations, which is calculated as follows: N -N ^OSC - A • 1 i \STF+ + 1 2, — i=2 In (64), Nose = (64) SW'- N * — i'-(N -i)! sw the number o f overlapping source wavelet combinations and Nsw = the number o f source wavelets. The total number o f particles (Ns) and corresponding Kalman filters for the P P D algorithm is calculated as follows: N =N xPExN / / s 0SC s , / \ T J Z 7 W . ar/ +(\-PE)xN (65) / s In (65), P E = the probability o f an event. Table 4.3 outlines the P P D filter formulation. Appendix F provides for more detail on the utilization of the F S M C D and the corresponding implementation o f the H M M filter within the P P D algorithm. Source wavelet present/not-present condition for each time increment. 54 T A B L E 4.3 PPD Step 1 2 3 4 FILTER FORMULATION Description Mathematical Representation Specify the number o f source wavelets (Nsw) and their corresponding phases and phase intervals, respectively. Specify the measurement-noise variance and time constant. Based upon these parameter formulate the J M L G S system and measurement equations. System Dynamics Calculate Ns and initialize the prior and transitional pdf for the event/no-event JMLGS. For each phase interval, initialize the prior and transitional pdf for the fixed-grid phases. A t each time index, draw samples for finite state Markov chain for the event/no-event conditions. If event is present (i.e., See (57) for the case of 3 source wavelets Measurement Equations See (58) for the case of 3 overlapping source wavelets. p(y ).p(yZ).&p(y \y - ), l 0 k J k l ij=l,2 p( (p])&p( <p\ | q>{ ), i = 1, ...,N w, and j S = l,...,Np], whereNpj = number offixedgrid value (i.e., phase interval resolution). k i _pli ^ J k y' = yl ) generate additional particles (i.e., Kalman filters) based upon the calculated combinations. k 5 k Propagate sinusoids based current time index (t =k?) and phase estimates. 55 Sinusoid(1..3) = sin(cot+f (1..3) ) k T A B L E 4.3 PPD Step (CONTINUED) FILTER FORMULATION Description Mathematical Representation w =w _,N(z \Z ,S ),, i i k i k k i k k i 1,...,N . S 7 8 Obtain asymptotically optimal estimate o f the state vector (i.e., estimate source wavelet responses for the three specified overlapping wavelets). Sampling Importance Re-sampling (SIR). Re-sample i f N EJR 9 <N i=\ N. 1 T Use x and H M M filter equations (Table k 3.2)toestimate/f/..j>) t. / 10 Utilize (18)-(20) to update the bank o f K F s 11 Let k = k+1 & iterate to step 4. 1) de Freitas [23] demonstrates that the importance weights for y\ are given by the predictive density p(z k \ z, _ ,y\. ) = N(z ;z' ,S (t 1 k k 56 k k ), whereNdenotes a Gaussian e) Evaluating the P P D A l g o r i t h m with Simulated Data The deconvolution o f minimum phase wavelets is a standard procedure, and the Weiner filter has been used extensively for this purpose [48]. The real challenge for seismic signal processors is to deconvolve mixed phased wavelets. For this reason, the performance o f the P P D was assessed by deconvolving seismograms that contained mixed phased Berlage wavelets [1]. The Berlage wavelet is defined as follows: w(t) = AH(t)tn e -at n (66) cos( 2/ft + (/)), where H(t) is the Heaviside unit step function [H(t) = 0 for t = 0 and H(t) = 1 for t >0]. The amplitude modulation component is controlled by two factors: the exponential decay term, a, and the time exponent, n. These parameters are considered to be non-negative, real constants. Figure 4.15 illustrates a mixed phased Berlage wavelet, with f = 50 Hz, n = 2, a = 100 and </> = 10°. The parameter, A, was set to give a maximum amplitude o f approximately 80 units, and the sampling rate was 0.05 ms. The Berlage wavelet o f Fig. 4.15 was then convolved with the closely spaced reflection coefficients shown i n F i g . 4.16 to give the output illustrated in F i g . 4.17. The synthetic seismogram shown i n F i g . 4.17 and all subsequently illustrated synthetic seismograms i n Section IV.B.2.e have additive GaussMarkov noise with a variance o f 50 units and a time constant o f 0.001 ms. 2 100; 60; ST 0 10 20 30 40 50 60 70 Time (ms) Figure 4.15. A mixed phased Berlage wavelet. 57 80 90 100 110 120 130 140 100 80 120 140 160 180 200 Time (ms) 220 Figure 4.16. Reflection coefficients convolved with Berlage wavelet shown in Fig. 4.15. 0 20 40 60 80 100 120 140 160 180 200 Time (ms) Figure 4.17. The output after convolving Berlage wavelet with reflection coefficients illustrated in Fig. 4.16 and adding measurement noise. In all subsequent post analysis, a 200 H z , eighth order, zero-phase shift, low-pass Butterworth frequency filter [31] is applied so that the signal-to-noise ratio o f the synthetic seismogram is significantly increased. This is due to the fact the P P D algorithm is separating overlapping source wavelets and not estimating high bandwidth reflection coefficients as in the case o f standard seismic deconvolution algorithms. For example, a major concern in frequency domain deconvolution is that one is working with bandlimited source wavelets and any high frequency measurement noise w i l l significantly degrade the ability to estimate the reflection coefficients. In this case, one cannot simply apply a low-pass frequency filter because this w i l l remove information associated with the high bandwidth 58 220 reflection coefficients. In the K F S D algorithm, it would be required to incorporate the applied low-pass frequency filter within the state-space formulation [42]. This would result i n a relatively complicated filter formulation becoming more complicated and unstable. In order to demonstrate the general concept o f the P P D algorithm, the dominant frequency and associated phases were determined interactively as opposed to utilizing a parameter estimation algorithm. For example, Fig. 4.18 illustrates the seismogram o f Fig. 4.17 after the 200 H z low-pass frequency filter was applied. Superimposed on this seismogram is a scaled sinusoid with dominant frequency o f 50 H z and corresponding phase o f 190°. A s is evident in Fig. 4.18, the sinusoid nearly overlaps the first portion o f the filtered synthetic time series o f Fig. 4.17; therefore, when implementing the P P D algorithm on the filtered seismogram o f Fig. 4.17, a dominant frequency o f 50 H z is specified along with an initial phase component of 190°. Time (ms) Figure 4.18. Estimating the dominant frequency and initial interactively. 59 phase component In Fig. 4.18, it is also apparent that another source wavelet is present around 104 ms, due to the fact that the seismogram no longer follows the sinusoid with a frequency o f 50 H z and a corresponding phase o f 190°. For these reasons, an additional scaled 50 H z sinusoid with a phase o f 100° is interactively introduced and summed with the sinusoid o f 190° at 104 ms. A s is shown in F i g . 4.19, the two sinusoids add such that they nearly have identical oscillations as the synthetic seismogram; therefore, the P P D algorithm is initialized with a dominant frequency o f 50 H z and two phases o f 190° and 100° and corresponding phase interval resolution o f ±3° (i.e., 187° to 193° and 97° to 103°, respectively). Figure 4.20 illustrates the output o f the P P D algorithm, where the two Berlage wavelets have been separated exactly. This is demonstrated by the error residual illustrated in F i g . 4.21. The error residual in Fig. 4.21 is defined to be the difference between the filtered seismogram o f Fig. 4.17 and the estimated Berlage wavelets o f Fig. 4.20. Time (ms) Figure 4.19. The addition of a 50 Hz sinusoid with a phase of 100° at a time of 104 ms. 60 100 120 140 160 180 200 220 Time (ms) Figure 4.20. The output of the PPD algorithm, where the Berlage wavelets have been separated. 80 100 120 140 160 180 200 Time (ms) Figure 4.21 The error residual, which is defined to be the difference between the filtered seismogram of Fig. 4.17 and the estimated Berlage wavelets of Fig. 4.20. 61 220 The P P D algorithm was next tested for its ability to deconvolve a closely spaced dipole. For this test, the Berlage wavelet o f Fig. 4.15 was convolved with the dipole shown i n Fig. 4.22 to give the output (with additive Gauss-Markov noise) as illustrated i n Fig. 4.24. 1| o> : 1 0.5- -a = 0- h E "* -0.5 : 0 20 40 60 80 100 120 140 160 180 200 220 Time (ms) Figure 4.22. The dipole reflection coefficients convolved with the Berlage wavelet shown in Fig. 4.15. 0 20 40 60 80 100 120 140 160 180 200 220 Time (ms) Figure 4.23. The output after convolving the Berlage wavelet of Fig. 4.15 with reflection coefficients illustrated in Fig. 4.22 measurement noise. 62 and adding the Gauss-Markov Figure 4.24 illustrates the output o f the P P D algorithm, where the algorithm did an impressive job in extracting the dipole Berlage source wavelets. Figure 4.25 shows the error residual, which is defined to be the difference between the filtered seismogram o f Fig. 4.23 and the estimated Berlage wavelets o f Fig. 4.24. u . . . . < u 0 20 i 40 , . , 60 80 100 120 140 160 180 , , 200 220 Time (ms) Figure 4.24. The output of the PPD algorithm, where the algorithm did an impressive job in extracting the dipole Berlage source wavelets. 4; 3i •3 0 20 40 60 80 100 120 140 160 180 200 Time (ms) Figure 4.25. The error residual, which is defined to be the difference between the filtered seismogram of Fig. 4.23 and the estimated Berlage wavelets of Fig. 4.24. 63 220 The next two sets o f synthetic seismograms were designed to test the P P D algorithm's ability to carry out blind seismic deconvolution on non-stationary source wavelets. In the first test, the zero-phase Ricker wavelet o f Fig. 4.11 is superimposed upon the Berlage wavelet illustrated in F i g . 4.15 (with an arrival time o f 99 ms) to give the output (with additive measurement noise) shown in F i g . 4.26. Figure 4.27 illustrates the output o f the P P D algorithm, where the algorithm did an impressive job in extracting the Berlage source wavelet from the overlapping zero-phase Ricker wavelet. I believe that this is the first example where a Ricker wavelet has been separated from a mixed-phase Berlage wavelet, utilizing a blind deconvolution algorithm. Figure 4.28 shows the error residual, which is defined to be the difference between the filtered seismogram o f Fig. 4.26 and the estimated wavelets o f Fig. 4.27. 100 120 140 160 180 200 220 Time (ms) Figure 4.26. T h e output after superimposing a Berlage wavelet with an arrival time of 99 ms and the zero-phase Ricker wavelet of F i g . 4.11. 64 180 | 0 . . I I . | . I 20 I : i | I 40 I I . I | I 60 I I I I | I 80 I I | I ! I I I | 100 120 Time (ms) I , i , . 140 I , , , I I 160 | i 1 | , | 180 , I I | | | | , , | 200 220 Figure 4.27. The output of the PPD algorithm, where the algorithm did an impressive job in extracting the Berlage source wavelet from the overlapping zero-phase Ricker wavelet. 4: 3J 0 20 40 60 80 100 120 140 160 180 200 Time (ms) Figure 4.28. The error residual, which is defined to be the difference between the filtered seismogram of Fig. 4.26 and the estimated source wavelets of Fig. 4.27. 65 The next non-stationary source wavelet test demonstrates the ability o f the P P D algorithms to deconvolve two Berlage wavelets, where the first wavelet has a frequency o f 50 H z and the second wavelet has a frequency o f 40 H z . This test illustrates the effects o f the Q value on the source wavelet. It is readily known that the higher frequencies are attenuated more rapidly than lower frequencies as the source wavelet travels through the earth. This test attempts to replicate this attenuation. Figure 4.29 shows a mixed phased Berlage wavelet with f= 40 Hz, n = 2, a = 100, <j> = 10°, and an arrival time o f 104 ms. This source wavelet was then superimposed upon the Berlage wavelet o f Fig. 4.15 with arrival time o f 99 ms. Figure 4.30 illustrates the result o f superimposing the time-variant Berlage source wavelets with additive measurement noise. The dominant frequencies and associated phases were again determined interactively, as was illustrated in Figs. 4.18 and 4.19. A s outlined in Fig. 4.18, the first Berlage source wavelet was determined to have a dominant frequency o f 50 H z and corresponding phase o f 190°. Figure 4.31 illustrates an additional scaled 40 H z sinusoid with a phase o f 100°, which is interactively introduced and summed with the 50 H z sinusoid o f phase 190° at 104 ms. A s is shown in Fig. 4.31, the two sinusoids add, such that they nearly have identical oscillations as the synthetic seismogram; therefore, the P P D algorithm is initialized with two dominant frequencies o f 50 H z and 40 H z , and two phases o f 190° and 100°, respectively. Figure 4.32 illustrates the output o f the P P D algorithm, where the two Berlage wavelets have been separated. The error residual is shown in Fig. 4.33. Time (ms) Figure 4.29. A mixed-phased Berlage wavelet with / = 40 Hz, n = 2, a = 100, and ^ = 10°. 66 0 20 40 60 80 100 120 140 160 180 200 220 Time (ms) Figure 4.30. The output after superimposing the time-variant Berlage source wavelets with additive measurement noise (variance = 50 units and time constant = 0.001 ms). 2 Time (ms) Figure 4.31. The addition of a 40 Hz sinusoid with phase of 100° at time of 104 ms to the sinusoid shown in Fig. 4.18. 67 0 20 40 60 80 100 120 140 160 180 200 22 Time (ms) Figure 4.32. The output of the P P D algorithm, where the algorithm did an impressive job in extracting the non-stationary Berlage source wavelets. 80 100 120 140 160 180 220 Time (ms) Figure 4.33 Error residual that is defined to be the difference between the filtered seismogram of Fig. 4.30 and the estimated source wavelets of Fig. 4.32. 68 f) Proposed Methodology i n Processing Seismograms with the P P D Algorithm A s previously outlined, the main challenge in implementing the P P D algorithm is the specification o f the principle phase components o f the seismogram. For this reason, I recommend that the overlapping source wavelets be sequentially and chronologically extracted from the seismogram under analysis. This mitigates the requirement o f prespecifying the principle phase components. The only phase which is required to be specified is that o f the current (in chronological order) source wavelet which is to be extracted from the seismogram. This information is readily available from the seismogram understudy. In this P P D filter wavelet extraction ( P P D - W E ) configuration there are only two possible overlapping source wavelets: the source wavelet to be extracted and the remaining seismogram time series data. Furthermore, for each event condition there are only two permutations. These permutations reflect the situation where only the source wavelet understudy is present at time index k (i.e., (1,0)) and the case where the source wavelet understudy is overlapped with other time series data at time index k (i.e., (1,1)). The two event condition allows the investigator to implement a significantly simplified P P D where only a 3 state F S M C D (i.e., noise, event condition (1,0) + noise, and event condition (1,1) + noise) is required. This avoids the requirement o f calculating the total number o f overlapping source wavelet combinations (i.e., (59)). Appendix G provides for more detail on the utilization o f the 3 state F S M C D within the P P D - W E . In the P P D - W E filter formulation the state equation defined by (62) is slightly modified as follows: £+l 1 A 0 0 0 x\ k 0 0 0 0 0" k+l 0 a 0 0 0 x2 0 b 2 0 0 0 u2 2k+l 0 0 1 A 0 x3 k + 0 0 0 0 0 0 k+\ 0 0 0 0 x4 k 0 0 0 V 0 u4 k+l_ 0 0 0 0 0 x l x2 x x4 x5 2 0 k a _ _x5 _ n k 69 0 0 k (67) k un k The P P D - W E algorithm models the amplitude modulating terms ( A M T ) (states xl and x3 ) k k o f the A M S source wavelet with a first order Taylor series approximation . The rate o f 6 change terms (states x2 and x4 ) are then approximated as Gauss-Markov processes. This k k allows for considerable flexibility and controllability when assigning prior to the amplitude modulating terms (e.g., smoothness). The time constant terms o f states x2 k and x4 are k important and fundamental parameters within the P P D - W E algorithm. It is desired that states x2 and x4 result in a smooth trajectory o f the amplitude modulation terms o f the k k A M S while at the same time allow for sufficient maneuverability so that the A M S source wavelet to be extracted ( A M S - E ) follows the oscillations o f the user specified sinusoid. The P P D - W E applies a linear range of possible time constant terms equal to the number o f specified particles (i.e., Kalman Filters) for states x2 and x4 . In the subsequent test bed k simulations, the range o f time constant 0.92, k terms varies from 0.6 ms (i.e., ai = a 2 highly maneuverable) to 6.22 ms (i.e., aj = a 2 = = 0.992, sluggish). The SIR (i.e., degeneracy check) portion o f the P P D - W E estimates the optimal value o f the time constant for states x2 and x4 within the range specified. k k Another important P P D - W E parameter to specify is the time index t* where it is assumed that the overlapping time series data has not yet arrived (i.e., only the A M S - E exist). This value is an approximation and is readily estimated as the end time o f the A M S - E initial phase estimate. Time index t* allows the P P D - W E to initially lock onto the A M S - E . Equation (68) outlines the P P D - W E measurement equation for the B S D problem when it is assumed that a maximum o f two source wavelets are overlapped at time index, k, for the event condition. k z = x ^k sin{eokA + qA ) + x3^ sin(a>kA + <j?3 ^ ) + x5 ^ (68) k In (68), co again represents the dominant frequency o f the source wavelet and /l k is the H M M filter phase estimate at time index, k, for the source wavelet to be extracted. The phase window resolution for the source wavelet is set quite small (e.g., ±3°) because the Please note that the amplitude modulating terms defined by (62) and (67) are forced to be positive. 6 70 phase o f the source wavelet to be extracted is easily determined from the seismogram with high accuracy. The H M M filter phase estimate at time index, k, for time series data overlapping the source wavelet to be extracted is defined by parameter /3k. The phase window resolution for f 3k is initially set at 1° to 360°. This range o f phases is generally refined based upon the initial values estimated from the first pass o f data. For example, i f after the first pass o f processing the seismogram it is found that the H M M estimates o f /3k range between 100° to 150°, then the P P D - W E is re-applied to the data set with f3k restricted to this range of phase. Major advantages o f the P P D - W E algorithm are the simplicity o f implementation, minimal parameter specification and there is no theoretical limit on the number o f overlapping source wavelets. Table 4.4 outlines the P P D - W E filter formulation. A disadvantage o f the P P D - W E algorithm is that any errors generate during the wavelet extraction process w i l l propagate as the seismogram is sequentially and chronologically processed. 71 T A B L E 4.4 PPD-WE Step 1 FILTER FORMULATION Description Mathematical Representation Identify the source wavelet to be extracted by specifying its corresponding dominant frequency and corresponding phase and phase resolution window. Specify the measurement-noise variance and time constant and the variance and time constant range for states x2 and x4 - Based upon these parameter formulate the J M L G S system and measurement equations. System Dynamics k k Specify Ns and initialize the prior and transitional pdf for the 3 state F S M C D as outlined in Appendix F (e.g., (A20). For the two phase intervals, initialize the prior and transitional pdfs for the fixed-grid phases. See (62). Measurement Equations See (63) for the case o f 3 overlapping source wavelets. P(y\). P(vl), P(y\h & P(A I y{-\)> P(<p\ )&p(<p\ \ <p{ ),i = l,2 andy = P ...,Npi, where N = number of fixed-grid value (i.e., phase interval resolution). PJ A t each time index, draw Ns samples for 3 y'k state F S M C D . Update Kalman filter measurement equations based upon y\ as outlined in Appendix F ((A21) and (A22). I f y\ =y\, then states xl and x3 are jittered k k [25] (i.e., add on U[-1.5,+1.5]) to facilitate a diversity o f particles. Propagate sinusoids based current time index (t =k?) and phase estimates. 72 Sinusoid(l ,3) = sin(cot+f (l,3) ) k T A B L E 4.4 (CONTINUED) PPD-WE FILTER FORMULATION Step Description Mathematical Representation Utilizing (13)-(17) outlined in Table 2.1 and updated sinusoids o f Step 5, propagate the system and measurement equations, calculate importance weights for particles and then update and normalize the weights. _ N(z \z ,S ),, j ^ '"" ' w\ = w | / f w | 1 Obtain asymptotically optimal estimate o f the state vector (i.e., estimate source wavelet to be extracted). Sampling Importance Re-sampling (SIR). Re-sample i f N <N eff T i = wk i ] k i k k i s N i=l N . . k ~ S ^k \ s x x 1 = 1 ^ ff e 1 ~TV i=i 9 Use x and H M M filter equations (Table k 2.2) to estimate/(l,3) . k 10 Utilize (18)-(20) to update the bank o f K F s 11 Let k = k+1 & iterate to step 4. 1) de Freitas [23] demonstrates that the importance weights for y\ are given b y the predictive density p(z k \ z _ ,y\. ) = N(z ;z' ,S' ), X:k x k k 73 k k whereNdenotes a Gaussian The implementation o f the P P D - W E algorithm is outlined in more detail by considering the following analysis o f synthetic seismograms. In the first test case, the mixed phased Berlage wavelet illustrated in Fig. 4.15 is convolved with the reflection coefficients outlined i n F i g 4.16 to give the output illustrated in Fig. 4.34. The synthetic seismogram shown in Fig. 4.34 is similar to that illustrated i n Fig. 4.17, but in this case the additive Gauss-Markov noise has a significantly higher variance of 400 units and a lower time constant o f 0.01 ms (i.e., 2 significantly higher correlation in measurement noise) which results i n a substantially lower signal-to-noise ratio. Time [ms] Figure 4.34. The output after convolving Berlage wavelet of Fig 4.15 with reflection coefficients illustrated in Fig. 4.16 and adding Gauss-Markov measurement noise with a variance of 400 units and a time constant of 0.01 ms. 2 A 100 H z , eighth order, zero-phase shift, low-pass Butterworth frequency filter is then applied to the synthetic seismogram shown in F i g 4.34 to give the output illustrated in F i g 4.35. Figure 4.35 shows the filtered time series data superimposed upon the raw synthetic seismogram illustrated i n F i g 4.34. In F i g 4.34 the initial 20 ms o f the synthetic seismogram time series data shown in Fig. 4.34 has been ignored due to the fact that there is no source wavelet information present within this time window. Figure 4.35 also illustrates the initial phase estimate o f 105° (±3°) o f the A M S - E source wavelet and a t* estimate of 7 ms. 74 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 Time [ms] • Original Data -AMS-ESinusoid(105°Phase) \ - Filtered Data Figure 4.35. The output after applying a 100 Hz, eighth order, zero-phase shift, lowpass Butterworth frequencyfilterto the synthetic seismogram shown in Fig. 4.34. The P P D - W E algorithm was then implemented on the synthetic seismogram illustrated in Fig. 4.34. The variance specified for states x2 and x4 was 90000 (units/s) . There were a 2 k total o f N s k = 500 particles (Kalman filters) specified. The first pass o f the P P D - W E algorithm implemented on the synthetic seismogram shown in Fig. 4.34 resulted in the initial/3 k phase estimates illustrated in Fig. 4.36. A s is evident from Fig. 4.36, the phases 30° to 50° are dominant for f 3 ; therefore, the P P D - W E was re-executed with this phase k restriction specified f o r / 3 . k 30" | | 0 10 i i • I | I I ' 1 20 . | I . 1 30 I • | 40 I . . i I | i 1 . \ 50 I | 60 I I I I I j 70 . i ) , I 80 , I , , 1 90 , j I I I I I | 100 Time [ms] Figure 4.36. The initial HMM/3 phase estimates forfirstpass of PPD-WE algorithm k on data shown in Fig. 4.35. 75 The second pass o f the P P D - W E algorithm provided for the A M S wavelet estimates illustrated in Figs. 4.37 and 4.38. Figure 4.37 shows the estimated and true Berlage source wavelets for the first reflection coefficient illustrated i n F i g 4.16. Figure 4.38 shows the estimated and true Berlage source wavelets for the second reflection coefficient illustrated i n F i g 4.16. A s is evident from Figs. 4.37 and 4.38, there is very close agreement between the estimated and true mixed-phased Berlage source wavelets. F i g 4.39 illustrates the A M T s for the two estimated A M S s shown in Figs. 4.37 and 4.38. The A M T s in Fig. 4.39 have smoothed trajectories while at the same they allow for significant maneuverability o f the corresponding A M S s . The P P D - W E algorithm obtained an optimal estimate time constant estimate o f 0.96 ms (i.e., a/ = a = 0.949) for states x2/ and x4 . 2 ( k One can readily synthesize reflection coefficients from the P P D - W E estimated source wavelets. A simplistic approached consist o f determining and averaging peak amplitude offsets and corresponding amplitude ratios from the estimated A M S s . If the source wavelet is stationary then the investigator can cross-correlate the estimated source wavelets and obtain relative arrival time estimates o f the reflection coefficients. For example, F i g . 4.40 illustrates the normalized cross-correlation function for estimated source wavelets shown i n Figs. 4.37 and 4.38. The peak value o f the cross-correlation function in F i g . 4.40 occurs at 4.7 ms; therefore the second reflection coefficient is estimated to arrive 4.7 ms after the previous (sequentially) reflection coefficient. The corresponding amplitude o f the reflection coefficient is estimated by averaging the dominant amplitude ratios o f the two estimated A M S s . F i g 4.41 illustrates the estimated and true reflection coefficients for the synthetic seismogram illustrated in Fig. 4.34 and the P P D - W E output shown i n Figs. 4.37 and 4.38. 76 100 I 0 j . 10 I I I • | i 20 I I i I | i 30 I I i I | 1 I I ! I 40 | ! I I 50 i I , I 60 I , , , j : 70 I i 1 I [ 80 I 90 I ! I I ,' 100 Time [ms] Figure 4.37. The estimated and true Berlage source wavelets for the first reflection coefficient illustrated in Fig 4.16. Time [ms] Figure 4.38. The estimated and true Berlage source wavelets for the second reflection coefficient illustrated in Fig 4.16. 77 i 0 i 5 . 10 , r , , , . 15 , i , , , 20 , , , , 25 i , 30 35 / 40 45 Time [ms] Figure 4.40. T h e normalized cross-correlation function for estimated source wavelets illustrated in Figs. 4 . 3 7 and 4 . 3 8 . Peak value occurs at 4.7 ms. 78 50 1 0.8 0.6 E 0.4 < 0.2A 0 10 ! I ! i i l i l | ! i ! i l j I i I ! ! | 12 14 16 1 18 MTTTTTTTTTi i ! 20 22 24 26 28 30 32 l iil ii 34 l 36 38 40 Time [ms] Estimated Reflection Coefficients True Reflection Coefficients Figure 4.41. Estimated and true Reflection coefficients for synthetic seismogram illustrated in Fig. 4.34. Fig. 4.42 shows a realization o f the three state F S M C D for time t = 40 ms. F S M C D values of 0, 1 and 2 correspond to the cases y ' = y \ , y ' = y \ and y k k Appendix G . 79 l k = y \ , respectively, as outlined i n In the second test bed analysis o f the P P D - W E algorithm, the mixed phased Berlage wavelet illustrated in Fig. 4.15 is convolved with the reflection coefficients outlined i n F i g 4.43 to give the output illustrated in Fig. 4.44. The synthetic seismogram shown i n Fig. 4.44 has additive Gauss-Markov noise with a variance of400 units and a time constant o f 0.01 ms. 2 1 0.8-I 0.6- 0.4- ^ 0.2: — 0 Q. E -0.2 < -0.4: -0.6-0.8-1- I 0 10 I I I 20 I | 30 i I ! I i | i 40 i I I 50 60 70 i : i 1 i 80 i i 90 100 Time [ms] Figure 4.43 Reflection coefficients convolved with Berlage wavelet shown in Fig. 4.15. 0 10 20 30 40 50 60 70 80 90 100 110 120 T i m e (ms) Figure 4.44. The output after convolving Berlage wavelet of Fig 4.15 with reflection coefficients illustrated in Fig. 4.43 and adding Gauss-Markov measurement noise with a variance of 400 units and a time constant of 0.01 ms. 2 80 A 100 H z , eighth order, zero-phase shift, low-pass Butterworth frequency filter is then applied to the synthetic seismogram shown in F i g 4.44 to give the output illustrated in F i g 4.45. Figure 4.45 shows the filtered time series data superimposed upon the raw synthetic seismogram illustrated in F i g 4.44. In F i g 4.45 the initial 20 ms o f the synthetic seismogram time series data shown in Fig. 4.44 has again been ignored due to the fact that there is no source wavelet information present within this time window. Figure 4.44 also illustrates the initial phase estimate o f 105° (±3°) o f the sinusoid o f the source wavelet to be extracted. Time index t* was again set to 7 ms. 200 Time [ms] | — Original Data Filtered Data AMS-E Sinusoid (105° Phase) | Figure 4.45. The output after applying a 1 0 0 Hz, eighth order, zero-phase shift, lowpass Butterworth frequency filter to the synthetic seismogram shown in Fig. 4.44. The P P D - W E algorithm was then implemented on the synthetic seismogram illustrated in Fig. 4.44. There were a total of N = 600 particles (Kalman filters) specified. The first pass S o f the P P D - W E algorithm implemented on the synthetic seismogram shown i n F i g . 4.44 resulted in the initial /3 k phase estimates illustrated in F i g . 4.46. A s is evident from Fig. 4.46, the phases 50° to 80° are dominant for /3 with this phase resolution specified f o r / 3 . k 81 k ; therefore, the P P D - W E was re-executed 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 Time [ms] Figure 4.46. The initial H M M / 3 k phase estimates for first pass of P P D - W E algorithm on data shown in Fig. 4.44. The second pass o f the P P D - W E algorithm provided for the A M S wavelet estimates illustrated i n Figs. 4.47 and 4.48. Figure 4.47 shows the estimated and true Berlage source wavelets for the first reflection coefficient illustrated i n F i g 4.43 Figure 4.48 shows the estimated residual wavelet and the actual residual wavelet generate by convolving the Berlage wavelet illustrated in Fig. 4.15 with the last two reflection coefficients illustrated in F i g 4.44. A s is evident from Figs. 4.47 and 4.48, there is very close agreement between the estimated and true time series responses for time interval 0 ms to 27 ms. After time 27 ms the estimated wavelets start to slightly diverge from the true wavelets. F i g 4.49 illustrates the A M T s for the two estimated A M S s shown in Figs. 4.47 and 4.48. The P P D - W E algorithm obtained an optimal estimate time constant estimate o f 1.14 ms (i.e., aj = a = 2 0.957) for states x2 and x4 . k k Subsequent to the extraction o f the first arriving Berlage source wavelet illustrated in F i g . 4.47, the P P D - W E algorithm is then applied to the estimated residual wavelet shown i n F i g . 4.48. Figure 4.50 illustrate the estimated residual wavelet shown in F i g . 4.48 with an initial phase estimate o f 44° (±3°) for the sinusoid o f the source wavelet to be extracted and a time index t* o f 9 ms. The first pass o f the P P D - W E algorithm implemented on the estimated residual wavelet shown in F i g . 4.50 resulted in the i n i t i a l / 3 k phase estimates illustrated in Fig. 4.51. A s is evident from F i g . 4.51, the phases 110° to 140° are dominant for f3 k therefore, the P P D - W E was re-executed with this phase resolution specified for 82 ; f3 . k 100 0 10 20 30 40 50 60 70 80 90 100 Time [ms] ^JEstimatedB^ Figure 4.47. The estimated and true Berlage source wavelets for the first reflection coefficient illustrated in Fig 4.43 100 0 10 20 30 40 50 60 70 80 90 100 Time [ms] ^^stimatedlRes^^ Figure 4.48. The estimated and true time series residual for the last two reflection coefficients illustrated in Fig 4.43 convolved with the Berlage source wavelet of Fig. 4.15. 83 Time [ms] — A M T of AMS-E AMT of Overlapping AMS "j Figure 4.49. The A M T s for the two estimated AMSs shown in Figs. 4.47 and 4.48. : / \ 4o; -o o_ 20 Vs. \ £ -20: < -40' /J '^ V— N t* = 9 ms 10 20 30 40 50 70 80 90 100 Time [ms] • Estimated Residual • AMS-E Sinusoid ^4£Riase)J Figure 4.50. The estimated time series residual of Fig. 4.48 with A M S - E sinusoid with phase of 44° superimposed. 84 < > I j' 10 15 20 25 30 35 40 45 50 55 60 65 70 75 1 1 1 1 80 1 1 1 1 85 ' 90 95 100 Time [ms] Figure 4.51. T h e initial H M M / 3 k phase estimates for first pass of P P D - W E algorithm on data shown in F i g . 4.50. The second pass o f the P P D - W E algorithm provided for the A M S wavelet estimates illustrated in Figs. 4.52 and 4.54. Figure 4.52 shows the estimated and true Berlage source wavelets for the second reflection coefficient illustrated in F i g 4.43 A s is evident from Fig. 4.52, there is very close agreement between the estimated and true Berlage source wavelets. Figure 4.53 shows the estimated (with 100 H z low-pass filter applied) and true Berlage source wavelets for the third reflection coefficient illustrated i n F i g 4.43 The P P D - W E algorithm did a good job i n obtaining accurate amplitude oscillations but it had difficulty in obtaining the correct amplitude scaling. This test bed outlines the current limitation o f the P P D - W E algorithm i n that any errors generate during the wavelet extraction process w i l l propagate as the seismogram is sequentially and chronologically processed. Figure 4.54 illustrates the A M T s for the two estimated A M S s shown i n Figs. 4.52 and 4.54. The P P D - W E algorithm obtained an optimal estimate time constant estimate o f 1.14 ms (i.e., ai = a = 0.957) for states x2 and x4 . 2 85 k k 0 10 20 30 40 50 60 70 80 90 100 Time [ms] Figure 4.52. The estimated and true Berlage source wavelets for the first second coefficient illustrated in Fig 4.43 80] 40 : 0 10 20 30 40 50 60 70 80 90 100 Time [ms] Estimated Berlage Wavelet — • True Berlage Waxelet | Figure 4.53. The estimated and true Berlage source wavelets for the third reflection coefficient illustrated in Fig 4.43. 86 10 20 30 40 50 60 70 80 90 100 Time [ms] -AMT of AMS-E • AMT of Overlapping AMS Figure 4.54. The A M T s for the two estimated AMSs shown in Figs. 4.52 and 4.53. Figure 4.55 illustrates the normalized cross-correlation function for the first and estimated second arriving Berlage source wavelets shown in F i g . 4.52. The peak value o f the crosscorrelation function in Fig. 4.55 occurs at 4.5 ms; therefore the second reflection coefficient is estimated to arrive 4.5 ms after the previous (sequentially) reflection coefficient. Figure 4.56 illustrates the normalized cross-correlation function for the first and estimated third arriving Berlage source wavelets shown in Fig. 4.53. The peak value o f the cross-correlation function in F i g . 4.56 occurs at 8.9 ms; therefore the third reflection coefficient is estimated to arrive 8.9 ms after the first reflection coefficient. 1- 0 25 30 35 40 45 50 Time [ms] Figure 4.55. The normalized cross-correlation function for the first and estimated second arriving Berlage source wavelet shown in Fig. 4.52. Peak value occurs at 4.5 ms. 87 0 5 10 15 20 25 30 35 40 45 Time [ms] Figure 4.56. The normalized cross-correlation function for the first and estimated third arriving Berlage source wavelet shown in Fig. 4.53. Peak value occurs at 8.9 ms. The corresponding amplitudes o f the reflection coefficients are estimated by averaging the four dominant amplitude ratios o f the first arriving Berlage source wavelet and the estimated overlapping source wavelets arriving at 24 ms and 25 ms, respectively. F i g 4.57 illustrates the estimated and true reflection coefficients for the synthetic seismogram illustrated in F i g . 4.44 and the P P D - W E output shown in Figs. 4.47, 4.52 and 4.53. 1 0.8 0.6 0.4 _=3 0.2-I o . 0E -0.2 ~| < -0.4' -0.6 -0.8 18 : : 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 Time [ms] -~ Estimated Reflection Coefficients - True Reflection Coefficients Figure 4.57. The estimated and true reflection coefficients for synthetic seismogram illustrated in Fig. 4.44. 88 34 g) Estimation of the Dominant Frequency of the Wavelet to be extracted within the PPD-WE Algorithm In the previously outlined P P D - W E algorithm the dominant frequency (DF) o f the source wavelet to be extracted ( S W E ) was assumed known. In order to make the P P D - W E truly blind, it is necessary to estimate the S W E ' s D F . The approach taken is to implement a H M M filter similar to the H M M filters utilized to estimate the A M S phases. In this P P D - W E filter formulation, a range o f possible D F s (e.g., 30 H z to 70 H z in 1 H z increments) and corresponding phases are pre-specified within a H M M filter formulation. A s the data is processed, the H M M frequency estimation ( H M M - F E ) filter obtains an optimal estimate o f the dominant frequency. The P P D - W E filter is re-executed with the estimate o f the H M M F E D F specified. The P P D - W E measurement equation for the case when there is only the S W E present is given as k z ^k =x sin{o)t + (pl k ), In (69), co represents the D F o f the S W E and / h where t = kA (69) is the corresponding H M M - F E D F estimate at time index k. A s previously outlined, knowing the D F one can easily obtain an initial estimate for the A M S - E sinusoid. For example, in F i g 4.35 the D F o f the A M S - E sinusoid was known (i.e., 50 Hz) and the initial phase was estimated to be 105° with a corresponding initial sinusoidal zero crossing of /' = 4.2 ms as illustrated in Fig. 4.58. 89 Time [ms] Filtered Data — — AMS-E Sinusoid (105° Phase) \ Figure 4.58. Illustration of the filtered synthetic seismogram of Fig. 4.34 with the A M S E sinusoid and corresponding initial zero crossing of 4.2 ms shown. A t the first zero crossing o f the A M S - E sinusoid, we have the following relationship (for negative amplitude first break): a>t' + <pl' = lS0° (70) From (70), it is clear that there is a linear relationship between the D F and the corresponding phase at the first zero crossing denoted by t' which is given as p l ' = 180°-fi>r' ( ) 71 B y utilizing (71) one can easily calculate a corresponding phase for any D F specified. For example, i n the case outlined i n F i g . 4.58 a D F o f 40 H z would result i n a corresponding phase o f 119.5° and a D F o f 70 H z result in a corresponding phase o f 74.2°. 90 The implementation o f the H M M - F E filter within the P P D - W E is outlined by considering the synthetic seismogram illustrated i n Fig. 4.34. In this case, the D F is assumed unknown and an initial frequency range o f 30 H z to 70 H z is specified with a increment resolution o f 1 H z and (71) is utilized to calculate the corresponding phases. Figure 4.39 shows the output of the H M M - F E filter after the first pass o f the H M M - F E filter and an initial value o f co = 35 Hz is specified. 25 Time [ms] Figure 4.59. Illustration of the output of the H M M - F E filter after processing the synthetic seismogram shown in F i g . 4.44. 91 A s is evident from F i g 4.59, the H M M - F E filter did an impressive job on locking onto the correct D F value o f 50 H z within a relatively short time frame. The D F value o f 50 H z is then specified within the P P D - W E algorithm and the synthetic seismogram illustrated i n Fig. 4.34 is then reprocessed to give the output shown in Figs. 4.37 and 4.38. In the next test bed analysis a synthetic seismogram is generated were there is time variance o f the A M T s and dominant frequency components o f the source wavelet. In this analysis there are four overlapping source wavelets. The four overlapping source wavelets B S W 1 , B S W 2 , B S W 3 , and B S W 4 are Berlage wavelets with dominant frequencies o f f = 55 Hz, f = 50 Hz, f = 45 Hz and / = 40 H, respectively, and parameters n = 2, a = 170, (j) = 60° specified. Figures 4.60, 4.61, 4.62, and 4.63 illustrate Berglage source wavelets B S W 1 , B S W 2 , B S W 3 , and B S W 4 , respectively. 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 Time [ms] Figure 4.60. Illustration of B S W 1 with parameters f = 55 Hz, n = 2, a = 170 and (j) = 6 0 ° specified. 92 I I i I i 0 I 1 ' 5 ' I I I I I I 10 I I I I 15 ' I P I 20 P I I I i I i 25 1I 30 I • I i ! | I • i I | • i i i 35 40 45 Time [ms] I . I | I 50 I I | . i . I | i 55 60 . . ; 65 i . ! , • . . i • 70 75 . • , Figure 4.61. Illustration of BSW2 with parameters f = 50 Hz, n = 2, a = 170 and <f> = 60° specified. 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 Time [ms] Figure 4.62. Illustration of BSW3 with parameters / = 45 Hz, n = 2, a =170 <t> = 60° specified. 93 | 80 and -60- TTTTTi 0 5 10 15 25 25 20 30 ' i i •I > M T r r r 35 40 45 35 40 Time [ms] 50 50 55 60 60 Figure 4.63. Illustration of BSW4 with parameters f=40 Hz, n=2,a= 65 H i 70 75 80 170 and <j> = 60° specified. The synthetic seismogram illustrated i n F i g 4.64 was generated b y convolving B S W 1 , B S W 2 , B S W 3 , and B S W 4 with the first, second, third, and fourth reflection coefficients, respectively, shown in F i g . 4.65 and summing the results. The time series data illustrated in Fig. 4.64 has additive Gauss-Markov measurement noise with variance o f 100 units and a 2 time constant o f 0.01 ms. Figure 4.66 illustrates B S W 1 , B S W 2 , B S W 3 , and B S W 4 after implementing the previously described convolution process. i 10 i i i ! 20 40 50 70 Time (ms) 80 i I ' •• > • i I I 100 110 120 i i I I I 130 140 Figure 4.64. Synthetic seismogram generated by summing time variant source wavelets BSW1, BSW2, BSW3, and BSW4. 94 0 10 20 30 40 50 60 70 80 90 100 Time [ms] Figure 4.65. The reflection coefficients utilized to generated synthetic seismogram illustrated in Fig. 4.64. 50- 25=3 = 0- -25- -50- " "I 10 15 I | I I I I | I I I I | i I I I 20 25 30 35 40 I I ! i 45 ! j 50 I M ; | I 55 ! II | I II I 60 J I I I 65 I 70 j I I I ! | ! f 75 Ii 80 i! i j 85 ! | I11111111111 90 95 100 Time [ms] •BSW1 BSW2 BSW3 BSW4 Figure 4.66. The output after convolving BSW1, BSW2, BSW3, and BSW4 with the first, second, third, and fourth reflection coefficients, respectively, shown in Fig. 4.65. 95 A 100 H z , eighth order, zero-phase shift, low-pass Butterworth frequency filter is then applied to the synthetic seismogram shown in F i g 4.64 to give the output illustrated i n F i g 4.67. Figure 4.67 shows the filtered time series data superimposed upon the raw synthetic seismogram illustrated in F i g 4.64 (with the initial 20 ms o f the synthetic seismogram time series data being ignored). The time parameters t* and t' were estimated to be 6 ms and 4.6 ms, respectively. In the implementation o f the H M M - F E an initial frequency range o f 30 H z to 58 H z was specified with a increment resolution o f 0.1 H z . In the specification o f the H M M - F E frequency range it is mandatory that the investigator does not specify frequency values which incorporate the seismogram's overall frequency components. For example, the synthetic seismogram i n F i g . 4.67 has overall frequency components ranging from 60 H z (e.g., peak 1 to peak 3 and peak 2 to peak 4) to 70 H z (e.g., peak 1 to peak 2). If these frequencies are incorporated into the H M M - F E , then the P P D - W E w i l l just track the seismogram response and not the principle phase components which comprise the seismogram. Figure 4.68 shows the output of the H M M - F E filter after the first pass o f the H M M - F E filter and an initial value o f m = 35 Hz is specified. A s is evident from F i g 4.68, the H M M - F E filter did an impressive job on locking onto the correct D F value o f 55 H z within a relatively short time frame. Time [ms] Figure 4.67. The output after applying a 100 Hz, eighth order, zero-phase shift, lowpass Butterworth frequency filter to the synthetic seismogram shown in Fig. 4.64. 96 Figure 4.68. Illustration of the output of the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.64. The P P D - W E algorithm provided for the A M S wavelet estimates illustrated i n Figs. 4.69 and 4.71 . Figure 4.69 shows the estimated and true B S W 1 source wavelet. Figure 4.70 7 illustrates the estimated residual wavelet and the actual residual wavelet (i.e., B S W 1 + B S W 3 + B S W 4 ) . A s is evident from Figs. 4.69 and 4.70, there is very close agreement between the estimated and true time series responses. Time [ms] — Estimated BSW1 True BSW1 Figure 4.69. The estimated and true B S W 1 source wavelets. Please note that the P P D - W E output has a 100 H z , eighth order, zero-phase shift, lowpass Butterworth frequency filter applied. 7 97 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 Time [ms] Estimated Residual (BSW2+BSW3+BSW4) — True Residual (BSW2+BSW3+BSW4) | Figure 4.70. T h e estimated residual wavelet and the actual residual wavelet (i.e., B S W 1 + B S W 3 + BSW4). Subsequent to the extraction o f the first arriving Berlage source wavelet illustrated i n F i g . 4.64, the P P D - W E algorithm is then applied to the estimated residual wavelet shown i n F i g . 4.70. Figure 4.71 illustrates the estimated residual wavelet shown i n Fig. 4.70 with the time parameters t* and t' set to 11 ms and 8.96 ms, respectively. In the implementation o f the H M M - F E an initial frequency range o f 30 H z to 55 H z (we know that the D F cannot be greater than the previous D F o f 55 H z ) was specified with an increment resolution o f 0.1 H z . Figure 4.72 shows the output o f the H M M - F E filter after the first pass o f the H M M - F E filter and an initial value o f co = 40 Hz specified. A s is evident from F i g 4.72, the H M M - F E filter locked onto the correct D F value o f 50 H z within 23 ms. 98 Bayesian Recursive Estimation 60^ 401 <u -a 20^: = 0= 3 a. E -20<£ •40| f = 8.9 ms -60; -80 0 t* = 11 ms ^ : 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 Time [ms] Figure 4.71. The estimated time series residual of Fig. 4.70 with the time parameters t* and t' set to 11 ms and 8.96 ms, respectively. 55-| 0 5 10 15 20 25 30 35 40 45 Time [ms] Figure 4.72. Illustration of the output of the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.71. Implementation o f the P P D - W E algorithm on the data shown in F i g . 4.71 resulted in the A M S wavelet estimates illustrated in Figs. 4.73 and 4.74. Figure 4.73 shows the estimated and true B S W 2 source wavelets. Figure 4.74 illustrates the estimated residual wavelet and the actual residual wavelet (i.e., B S W 3 + B S W 4 ) . A s is evident from Figs. 4.73 and 4.74, the estimated responses deviate from the true responses between approximately 17 ms to 27 ms. 99 80 Figure 4.73. The estimated and true BSW2 source wavelets. |II 0 l i i | I I I 5 10 •I| 15 I 1 • I I | 20 : •• i 25 ( I | I I : I I I . . I I 30 I j I !i I 1 I I I I I > I I: I I| i I : I I : . i I I I i i 35 40 45 50 55 60 65 70 Time [ms] Estimated Residual (BSW3+BSW4) True Residiual (BSW3+BSW4) | 1 75 Figure 4.74. Estimated residual wavelet and the actual residual wavelet (i.e., BSW3+BSW4). 100 1 I . . | I I , I i | 80 | Subsequent to the extraction o f Berlage source wavelets B S W 1 and B S W 2 , the P P D - W E algorithm is then applied to the estimated residual wavelet shown in F i g . 4.74. Figure 4.75 illustrates the estimated residual wavelet shown in Fig. 4.74 with the time parameters t* and t' set to 17 ms and 15 ms, respectively. 101 Bayesian Recursive Estimation A n initial frequency range o f 30 H z to 47 H z was specified with an increment resolution o f 0.1 H z within the H M M - F E filter. Trie 50 H z frequency component was avoided due to the fact that the overall seismogram had a dominant frequency o f approximately 50 H z (peak 1 to peak 3). Figure 4.76 shows the output o f the H M M - F E filter after the first pass o f the H M M - F E filter and an initial value o f co = 40 Hz specified. A s is evident from F i g 4.77, the H M M - F E filter locked onto a D F value o f approximately 44.5 H z within 30 ms. | 0 5 10 15 . 20 ! I i I | I i I i 25 . j | 30 35 i I I . I ; I > j I 40 . | I 45 I I i I | I M 50 I I | I II I 55 60 Time [ms] Figure 4.76. Illustration of the output of the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.75. Implementation o f the P P D - W E algorithm on the data shown in F i g . 4.76 resulted in the A M S wavelet estimates illustrated in Figs. 4.77 and 4.78. Figure 4.77 shows the estimated and true B S W 3 source wavelets. Figure 4.74 illustrates the estimated and true B S W 4 source wavelets. A s is evident from Figs. 4.77 and 4.78, the P P D - W E algorithm but had some difficulty in estimating the A M T s . This is due to the fact that any errors generate during the wavelet extraction process will propagate as the seismogram is sequentially and chronologically processed. To the best o f m y knowledge, the previously processed synthetic seismogram is the first example o f the blind extraction o f time variant overlapping source wavelets with closely spaced dipole reflection coefficients. The low signal-to-noise ratio synthetic seismogram contained two dipole reflection coefficients offset within 10 ms and time variant source 102 wavelets with an approximate duration o f 40 ms. These results have significant implication for site characterization investigations (SCI) which utilize vertical seismic profiling (VSP). In many SCI there exists significant in-situ impedance contrast due to remediation efforts (e.g., insertion o f stone columns utilizing vibro-compaction). This makes the post analysis o f the V S P data appreciably difficult due to the closely spaced reflection coefficients and the time variant source wavelets. The P P D - W E algorithm could significantly mitigate the difficulty i n separating the closely spaced overlapping source wavelets. II 30 35 Iii Ii 40 I | H 45 I , , | I I , , , | • I , I I j , , I , 50 55 60 65 70 . 75 80 Time [ms] — Estimated BSW3 • True BSW3 Figure 4.77. T h e estimated and true B S W 3 source wavelets. 6040; T3 =3 "a. E < 20; _ / / fl- V \\ -20; -40-6010 15 20 25 30 35 40 45 i " 1 50 Time [ms] • Estimated BSW4 — True BSW4 Figure 4.78. T h e estimated and true B S W 4 source wavelets. 103 55 60 65 70 75 80 h) Utilization of the P P D - W E A l g o r i t h m within Standard Frequency D o m a i n Deconvolution Techniques A standard frequency domain methodology i n estimating the reflection series, ju , is the k water level technique ( W L T ) [51]. The W L T is outlined by considering (7), representing the Z transform by the Fourier transform (i.e., z = e " ) and rearranging terms as follows: 10 *r«; =^ (72) S(co) where ¥(co) denotes the Fourier transform o f the reflection series /J(t). Theoretically, for a x known source wavelet, one could simply implement (72) and calculate ¥(co). x series, n , k is then estimated The reflection by taking the inverse Fourier transform of ^(co). Unfortunately, due to inaccuracies in the specification o f the source wavelet, the bandlimited nature o f the source wavelet and additive measurement noise, the implementation o f (72) is highly unstable and inaccurate. To mitigate the previously outlined limitations, (72) is modified b y firstly multiplying the numerator and denominator by the complex conjugate o f S(co) (denoted as S*(co)) and secondly, b y introducing an additive scalar value to the denominator which is referred to as the water level (A) [51]. Implementation o f these two modifications to (72) gives Z(co)S*(co) S(co)S*(co) + A ^Z(co)S*(co) ( ? 3 ) P (co) + A s where P (co) denotes the power spectrum o f the source wavelet (i.e., the Fourier transform s o f the autocorrelation o f S(t)). In general terms, the setting o f the water level is a trial and error approach. A s A —> 0 , the resulting estimated reflection coefficients approach Dirac delta functions. When A » P(co) the resulting estimated reflection coefficients become significantly bandlimited and the result converges to the Fourier transform o f the crosscorrelation between the recorded seismogram Z(co)S*(co)). 104 and the source wavelet (i.e., The implementation o f the W L T in conjunction with the P P D - W E technique is illustrated by considering the synthetic time series illustrated i n F i g . 4.79. The simulated time series shown in this figure is a typical seismogram which one may encounter and it was generated without m y prior knowledge o f the source wavelet(s) or reflection series. 0 100 50 150 200 250 300 350 400 450 500 Time [ms] Figure 4.79. A typical synthetic seismogram. The amplitude spectrum o f the seismogram shown in Fig. 4.79 is illustrated in Fig. 4.80. A s is evident from Fig. 4.80, the dominant seismic bandwidth resides between 40 H z to 60 H z . 0.081 — Frequency[Hz] Figure 4.80. T h e amplitude spectrum of the seismogram illustrated in F i g . 4.79. 105 A 100 H z low-pass Butterworth frequency filter was then applied to the seismogram o f Fig. 4.79. Figure 4.81 illustrates the filtered seismogram (black line) superimposed upon the noisy seismogram (light grey plot) o f Fig. 4.79. 0.8 o.6; 0.4; • 3 0 2 0 §- -0.2; •* -0.4; -0.6! -0.8 -1-1 50 100 150 200 250 300 350 400 450 Time [ms] Figure 4.81. The superposition of filtered seismogram onto seismogram shown in Fig. 4.79. The P P D - W E algorithm was then implemented on the filtered seismogram so that the first arriving source wavelet could be extracted and estimated. The first 45 ms o f the filtered time series was ignored due to the fact no signal was present. A s a result o f the estimated frequency bandwidth o f the seismogram (i.e., F i g . 4.80), the H M M - F E component o f the P P D - W E algorithm had a frequency range o f 40 H z to 60 H z specified. Figure 4.82 illustrates the estimated dominant frequency output o f the H M M - F E . A s is evident from F i g . 4.82, the H M M - F E locks onto a dominant frequency range o f approximately 49 H z to 51 H z within 5 ms from the onset o f the first arriving source wavelet. 106 500 45 : . . •. . )I .,,,|,., .,, 0 5 10 ,,,,,, 15 20 |, , , I , , : I , ,, 25 30 35 4 Time [ms] Figure 4.82. Illustration of the output of the H M M - F E filter after processing the filtered seismogram shown in F i g . 4.81. Figure 4.83 shows the P P D - W E estimated first arriving source wavelet superimposed upon the true source wavelet. A s is evident from Fig. 4.83, the P P D - W E algorithm estimated the correct dominant frequency o f approximately 50 H z and also did a very impressive job in estimating and extracting the first arriving source wavelet. The estimated responses after 105 ms are a P P D - W E filter residual where the algorithm is estimating overlapping source wavelet responses with a similar phase and dominant frequency to the currently estimated source wavelet. For this reason and due to the fact that the source wavelet w i l l not likely, from a physics point o f view, decay to zero and start up again, the estimated time series beyond 105 ms is set to zero. F i g . 4.84 illustrates the estimated source wavelet o f Fig. 4.83 with the time series beyond 100 ms set to zero. 107 108 109 A s opposed to going through the lengthy process o f separating all possible overlapping source wavelets, the first deconvolution attempt was performed using the estimated source wavelet o f Fig. 4.84 and the W L T . If there is minimal source wavelet variation within the time series, the estimated reflection coefficients w i l l be very similar in shape. Figure 4.85 illustrates the superimposition o f the true reflection series onto the estimated reflection series for the synthetic seismogram shown in Fig. 4.79. For these results, an 8 order Butterworth th low pass filter was applied to the noisy seismogram and the water level was set to 0.2% o f the maximum value o f the power spectrum o f the seismogram. A s is evident from Fig. 4.85, the W L T in conjunction with the P P D - W E algorithm did an excellent job in recovering the true reflection series. In addition, we can assume that a stationary source wavelet is present due to the nearly identical shape of the estimated reflection coefficients. 1.0 0.8 0.6 £ 0.4 3 0.2 I '• if! 0 f UH I ti -0.2 'A iii) nil -H *i f 1 -0.4 0 _1 I I I • ' 50 _l 100 I I l_ 150 -1—I—I—1—I 200 250 I I I I 300 I I I I 350 1 I I I 400 Time [ms] Figure 4.85. Superposition of the true reflection series (black line) onto the estimated reflection series (dotted line). 110 i) Monte Carlo Simulations There two important Monte Carlo simulation ( M C S ) considerations in the implementation of the P P D - W E algorithm. 1) Step 4 o f Table 4.4, where, i f y\ =y\ , then states xl k and x3 k are jittered by means o f a random number generator. 2) The response o f the algorithm to varying noise realizations. Figure 4.86 shows ten P P D - W E estimates for the noisy synthetic seismogram shown in Fig. 4.79. The response identified b y the dotted line differs from the other nine estimates considerably and is, therefore, ignored. The remaining nine responses are averaged to give the source wavelet estimate. 40 50 Time [ms] Figure 4.86. Ten P P D - W E estimates for noisy synthetic seismogram shown in Fig. 4.79. The estimate illustrated by the dotted line is explained in the text. The bold light grey line identifies the true source wavelet. Ill The ability o f P P D - W E algorithm to respond to varying noise realizations is also o f concern. Figure 4.87 illustrates an example o f a synthetic seismogram where the noise free seismogram o f Fig. 4.79 has additive Gauss-Markov noise with a time constant o f 0.001 ms and a variance o f 0.005 units . 2 1:1 i 0 10 20 30 40 50 , . . . 60 , i 8 70 80 90 Time [ms] Figure 4.87. Synthetic seismogram where the noise free seismogram of Fig. 4.79 has additive Gauss-Markov noise with a time constant of 0.001 ms and a variance of 0.005 units . 2 Figure 4.88 illustrates the output o f the P P D - W E algorithm for ten additive Gauss-Markov noise (time constant o f 0.001 ms and a variance o f 0.005 units ) realizations. A s is evident 2 from Fig. 4.88, the P P D - W E had high repeatability when processing the ten synthetic seismograms. 112 i 100 0.8 i : Time [ms] Figure 4.88. P P D - W E source wavelet estimates for ten additive G a u s s - M a r k o v noise (time constant of 0.001 ms and a variance of 0.005 units ) realizations. T h e bold light 2 grey line identifies the true source wavelet. 113 j) Sensitivity Analysis Due to the fact that the H M M filters are utilized to refine inaccurately specified initial phase estimates, the sensitivity limitation o f the P P D and P P D - W E algorithms resides in their ability to separate source wavelets which have similar phases or have phase differences near 180°. Figure 4.89 illustrates two overlapping Berlage source wavelets where the first wavelet has an arrival time o f 99 ms and the second Berlage wavelet has an arrival time o f 108 ms. These two source wavelets have associated phases o f 190° and 20°, respectively. This results in a phase difference o f 170° that is 10° from the indiscernible value o f 180°. Figure 4.90 shows the output o f the P P D algorithm where the algorithm was able to separate the two source wavelets, but has some difficulty in estimating the modulating amplitude term. This is highlighted i n Figs. 4.91 and 4.92. In Fig. 4.91 the error residual between the estimated P P D Beralge wavelet and the true Berlage wavelet for the source wavelet with arrival time o f 99 ms is shown. Figure 4.92 illustrates the error residual between the estimated P P D Beralge wavelet and the true Berlage wavelet for the source wavelet with an arrival time o f 108 ms. Due to the nature o f these and similar results, it is believed that the current formulation o f the P P D and P P D - W E algorithms cannot sufficiently separate overlapping source wavelets which have phase difference less than 10° to 20° or between 160° to 180°. 0 20 40 60 80 100 120 140 160 180 200 Time (ms) Figure 4.89. The output after overlapping two Berlage source wavelets with additive measurement noise. These two source wavelets have a phase difference of 170°. 114 Bayesian Recursive Estimation 100 0 20 40 60 80 100 120 140 160 180 200 Time (ms) Figure 4.90. The output of the PPD algorithm, where the algorithm was able to separate the two source wavelets, but has some difficulty in estimating the modulating amplitude term. 1 0 0 3 80; 60; ~ 40; ~ 20; o; J -20= 40; -60; 0 20 40 60 80 100 120 140 160 180 200 220 Time (ms) Figure 4.91. The error residual between the estimated PPD Beralge wavelet and the true Berlage wavelet with an arrival time of 99 ms. 115 100 Time (ms) Figure 4.92. The error residual between the estimated P P D Beralge wavelet and the true Berlage wavelet with an arrival time of 108 ms. The P P D and P P D - W E phase resolution limitation window is referred to as the null space. A ±20° resolution is equivalent to a theoretical reflection coefficient resolution o f ±1.1 ms for a source wavelet with a 50 H z dominant frequency (i.e., it is possible to separate overlapping source wavelets which are separated by a minimum of 1.1 ms). Figure 4.93 illustrates the null spaces for a 50 H z source wavelet (T = 20 ms) and null space resolution o f ±20°. In general terms, overlapping wavelets which arrive within the null spaces cannot be separated from the source wavelet to be extracted. A n extracted source wavelet which has overlapping wavelets occupying the null spaces might be evident by a peculiar looking A M T and corresponding extracted source wavelet. 116 Figure 4.93. Illustration of the null spaces with a 50 H z (T = 20 ms) source wavelet. Overlapping wavelets which arrive within the null spaces cannot be separated from the source wavelet to be extracted. 117 3) Implementation of the P P D - W E Algorithm on Real Data The P P D - W E algorithm was implemented on real data acquired during a seismic cone penetration test (SCPT). The P P D - W E is a general purpose algorithm for blind deconvolution that is particularly well suited to wavelet estimation when no assumptions concerning the reflectivity can be made. For example, a minimal set o f reflection coefficients is not conducive to assigning an associated pdf. Hence real data from a S C P T is analyzed due to the sparseness o f the reflection coefficients (i.e., two to three reflection coefficients within the seismogram). The S C P T has proven to be a very valuable geotechnical tool in facilitating the determination o f low strain (<10~ ) in-situ compression (P) and shear (S) wave velocities. 6 The P-wave velocities (Vp) and S-waves velocities (Vs) are directly related to the soil elastic constants o f Poisson's ratio, shear modulus, bulk modulus, and Y o u n g ' s modulus. The estimation o f elastic constants is essential in geotechnical foundation designs [30]. Another important use o f estimated shear wave velocities in geotechnical design is in the liquefaction assessment o f soils [5]. Details o f the seismic cone, the downhole test procedures, and comparisons with the crosshole results at several sites have been described by Campanella et al. [20]. In general terms, the seismic cone is advanced to the depth o f interest using a hydraulic reactionary pushing. The advance is halted at one meter (or other such increment) intervals. When the cone is at rest, a seismic event is caused at the surface using a hammer blow or explosive charge, causing seismic waves to propagate from the surface through the soil to be detected by seismic sensors installed i n the cone penetrometer. This event is recorded and the penetrometer is advanced another increment and the process is repeated. B y determining the relative seismic arrival times average velocities are calculated over the depth increment under study. The relative arrival times are obtained by crosscorrelating the recorded source wavelets recorded at each depth increment [14]-[16]. The interval velocities are then calculated as the ratio between the relative arrival time differences and the corresponding relative travel path differences between successive depths (i.e., v = Ad/AT, where AT is the relative arrival time 118 difference, Ad is the relative travel path interval velocity). difference and v is the corresponding SCPT It is o f paramount importance to extract first arriving source wavelets from each trace so that crosscorrelation is meaningful when obtaining relative arrival times. The data analyzed in this section was acquired with a triaxial system configuration. A triaxial sensor configuration is utilized so that full waveform analysis can be carried out and the possibility o f rod rotation can easily be taken into account [17]. Details o f the S C data acquisition system and seismic sensors utilized are given i n [18]. For the S H wave analysis, the preprocessing o f the seismic time series data captured on the X and Y axes is three-fold. 1) A p p l y bandpass frequency filter. 2) Rotate the X and Y axes responses onto the full waveform axis utilizing hodograms and polarization analysis [10], [17]. 3) A p p l y exponentially decaying windows to the front and end o f the seismic responses on the recorded time series so that the S / N is increased. Fig. 4.94 illustrates a typical S C P T configuration where a S H source is located at the outriggers o f the in-situ testing vehicle. A triaxial S C sensor configuration is utilized and seismic source wavelets are acquired at constant depth intervals so that the in-situ S H wave interval velocities can be calculated. In F i g . 4.94 the interval velocity travel time ATi is calculated by crosscorrelating the source wavelets captured at each depth increment. The seismic cone data analysed in this section is referred to as SC 64 and it was acquired in N e w Zealand on February 13 , 2006 near the Tauranga Eastern Motorway by Perry Drilling th Ltd. o f Tauranga N e w Zealand. The seismic source utilized was a horizontal shear (SH) hammer source. The (SH) source wavelets are generated at the outriggers which were positioned 2.3 metres from the centre o f the rod strings. The source consists o f a sledge hammer horizontally impacting point source steel beams located underneath the outriggers and an electrical contact trigger is utilized. T w o stacked S H sources were generated for each depth increment. Fig. 4.95 illustrates raw X axis time series data captured from S C P T test site SC 64. F i g . 4.96 illustrates raw Y axis time series data captured from S C P T test site SC 64. A s is evident from Figs. 4.95 and 4.96, there are source multiples present within the captured time series data. The overlapping source wavelets present significant difficulties in implementing the crosscorrelation technique. To address this problem the P P D - W E algorithm is utilized to extract the first arriving source wavelet and subsequent to this the crosscorrelation 119 technique is implemented. F i g . 4.97 illustrates the vertical seismic profile (VSP) after processing the time series data shown in Figs. 4.95 and 4.96. Figure 4.94. T y p i c a l S C P T configuration. 120 Figure 4.95. R a w X axis V S P from S C P T hole SC 64. 121 1 , , , . 20 40 60 , , . 80 100 , , ; 120 140 160 180 Time [ms] Figure 4.96. Raw Y axis V S P from SCPT hole SC 64. 122 200 220 240 ! 260 280 300 18 20 40 BO BO 100 120 140 160 Time [ms] 180 200 220 240 260 280 Figure 4.97. V S P for the processed seismic data captured at test site SC 64 for depths 2 m to 17 m. T h e r e is clear evidence of overlapping source wavelets. 123 300 In the implementation o f the P P D - W E technique there are slightly differing possible source wavelet estimation realizations depending upon the specification o f the initial P P D - W E filter parameters o f / (zero crossing), /* (time to source wavelet overlap), and R , (variance o f k measurement noise) and the utilization o f Monte Carlo techniques to obtain realizations o f y k ~ I y{-i) ^ an Ji t t e r m e particles (i.e., Step 4 o f Table 4.4). For this reason, a Monte Carlo technique is utilized which allows the investigator to vary the input filter parameters and subsequently obtain many estimates o f the source wavelet. The source wavelet is then derived by averaging the subsequent P P D - W E source wavelet estimates [13]. In this P P D - W E technique ( P P D - W E M C ) an additional time parameter, t , is introduced. 2 Unfortunately, there are some situations where the A M T o f the source wavelet w i l l slowly start to diverge from the true value due to the degeneracy check. T o mitigate this effect, the degeneracy check is turned off after time t which results in a noisier estimated source 2 wavelet but the diversity o f the particles is maintained. In the P P D - W E M C algorithm the investigator initially specifies minimum t^ min seismogram zero crossing) and t 2min (first parameters. These parameters are modified within each iteration o f the P P D - M E M C algorithm (source wavelet estimate) according to the following two equations: where t ~ N(0,0.4) and t ~ N(0,225). In addition, the measurement noise variance is 2 increased from an initial user specified minimum, R , according to a specified increment min value, R inc (i.e., R = R min + R ). inc The frequency window specified for the H M M - F E filter portion o f P P D - W E algorithm was estimated based upon reviewing the amplitude spectrum o f the acquired seismograms. For example, Fig. 4.98 shows the amplitude spectrum o f the seismogram acquired at 2 m. A s is evident from Fig. 4.98, the dominant seismic bandwidth resides between 45 H z to 60 H z . 124 0.04 T3 =5.0.02E < 40 50 60 Frequency [Hz] Figure 4.98. T h e amplitude spectrum of the seismogram captured at 2 m. In the implementation o f the P P D - W E algorithm on the S C data acquired at test site SC 64, the minimum zero crossing was readily available from the processed seismogram as is illustrated in Fig. 4.99. The H M M - F E frequency window was set to range from 45 H z to 60 Hz. There were thirty five primary source wavelet estimates for each run o f the P P D W E M C algorithm and the minimum noise variance, R min , was set to 0.1 with an increment increase o f 0.1. This results in a measurement noise variance o f 3.6 for the final run o f the P P D - W E M C algorithm. 25 Figure 4.99. Estimating t min 30 35 Time [ms] 40 45 for the seismogram acquired at 12 m. 125 50 55 60 Fig. 4.100 illustrates the extracted (utilizing P P D - W E technique) primary source wavelets V S P for depths 2 m to 17 m o f S C P T hole SC 64. The time series recorded at 2.0 m was assumed to contain no source multiples. Appendix H provides illustrations o f the P P D - W E estimated source wavelets, averaged source wavelets and corresponding residual wavelets for test site SC 64 and depths 3 m to 17 m. A s is evident from F i g . 4.100, the P P D - W E technique did an impressive job in extracting the primary source wavelets. The crossocorrelation technique can now be implemented so that accurate interval S H wave velocities can be obtained. The estimated interval velocities from the time series data given in Fig. 4.100 are outlined i n Table 4.5. The traces shown in F i g . 4.100 have very high calculated interval correlation coefficients as shown i n Table 4.5. Table 4.6 outlines the estimated interval velocities i f the processed seismograms o f F i g 4.97 are crosscorrelated when obtaining the interval arrival times. A s is evident from Table 4.6, there are significant deviations from the estimated interval velocities o f Table 4.5 and the calculated correlation coefficients are corresponding lower. It should be noted that the shear and compression wave velocities are squared in deriving the elastics constants; therefore, variations in the estimated velocities can cause appreciable errors in the calculation o f the elastic constant. The crosscorrelation interval arrival times derived from the processed seismograms o f Fig. 4.97 are also sensitive to the time window specified for the crosscorrelation function. Fig. 4.101 illustrates the residual wavelets derived b y subtracting the estimated primary source wavelets outlined i n Fig. 4.100 from the processed seismograms shown in Fig. 4.97. A s is evident from F i g . 4.101 and the residual wavelet plots shown in Appendix H , the majority o f the seismograms have only one source multiple present. For the residuals at 16 m and 1 7 m there is evidence that there is more than two source wavelet multiples present in the acquired seismogram. For example, F i g . 4.102 shows the estimated reflection coefficients utilizing the W L T i f the estimated source wavelet at 16 m is applied to the processed seismogram recorded at 16 m. F i g . 4.103 shows the estimated reflection coefficients utilizing the W L T i f the estimated source wavelet at 17 m is applied to the processed seismogram recorded at 17 m. 126 Time [ms] Figure 4.100. Extracted (utilizing P P D - W E M C technique) primary source wavelets VSP for depths 2 m to 17 m of SCPT hole SC 64. 127 T A B L E 4.5 INTERVAL V E L O C I T I E S ( P P D - W E E X T R A C T E D S O U R C E W A V E L E T S ) F R O M T H E C R O S S C O R R E L A T I O N T E C H N I Q U E F O R S C P T sc 6 4 Interval Depth (m) Crosscorrelation Interval Velocity Estimate (m/s) 2.0-3.0 106 0.90 3.0-4.0 120 0.96 4.0-5.0 96 0.97 5.0-6.0 118 0.98 6.0-7.0 141 0.96 7.0-8.0 156 0.99 8.0-9.0 150 0.95 Correlation Coefficient 9.0-10.0 139 0.95 10.0-11.0 143 0.96 11.0-12.0 122 0.98 12.0-13.0 110 0.99 13.0-14.0 168 0.96 14.0-15.0 148 0.996 15.0-16.0 166 0.99 16.0-17.0 216 0.995 T A B L E 4.6 INTERVAL V E L O C I T I E S ( P R O C E S S E D SEISMOGRAMS) FROM THE CROSSCORRELATION T E C H N I Q U E F O R SCPT SC 64 Interval Depth (m) Crosscorrelation Interval Velocity Estimate (m/s) 2.0-3.0 93 0.90 Correlation Coefficient 3.0-4.0 104 0.86 4.0-5.0 92 0.93 5.0-6.0 114 0.97 6.0-7.0 121 0.81 7.0-8.0 164 0.97 8.0-9.0 152 0.84 9.0-10.0 167 0.94 10.0-11.0 144 0.97 11.0-12.0 118 0.77 12.0-13.0 108 0.91 13.0-14.0 154 0.80 14.0-15.0 167 0.96 15.0-16.0 148 0.85 16.0-17.0 287 0.79 128 "I 25 • •• ! ! 50 1 i ! i 1 75 ! •• I ! 1 100 1 1 1 : 1 : 125 Time [ms] i ! ! ; 150 1 , , , 1 175 , , , , 1 200 Figure 4.101. Residual wavelets V S P for depths 3 m to 17 m of SCPT hole SC 64. 129 , , : , 1 225 150 200 250 300 350 Time [ms] Figure 4.102. Estimated reflection coefficients for depth 16 m utilizing the W L T . 150 200 250 300 Time [ms] Figure 4.103. Estimated reflection coefficients for depth 17 m utilizing the W L T . 130 350 Fig. 4.104 shows the superposition o f the P P D - W E estimated source wavelets outlined i n Fig. 4.100 (depths 4 m to 17 m) and corrected zero time offset. F i g . 4.105 illustrates the source wavelet obtained by averaging the estimated source wavelets illustrated in Fig. 4.104. The estimated source wavelets outlined in Figs. 4.104 and 4.105 are scale to that o f the central trough. 0 25 50 Time [ms] Figure 4.104. Superposition of estimated p r i m a r y source wavelets for depths 4 m to 17 m. 0 25 Time [ms] Figure 4.105. Estimated source wavelet (black line) obtained by averaging estimated PPD-WE source wavelets shown in F i g . 4.104. T h e light grey lines signify m a x i m u m and m i n i m u m bounds based upon the response outlined in F i g . 4.104. 131 the Test Site SC 64 was located approximately 7 m to 9 m to the east o f a pond that is 1 m to 2 m deep. The test was conducted 1 m below the bund surrounding the pond which is used as an access road. The bund is constructed o f dense material and the S H source for S C P T SC 64 was positioned approximately 4m from the bank o f the bund. Table 4.7 outlines the estimated time offsets for the secondary arriving source wavelet for test site SC 64. These values were estimated visually by reviewing the figures in Appendix H which illustrate the estimated P P D - W E source wavelet superimposed upon the residual wavelet. T A B L E 4.7 E S T I M A T E D S E C O N D A R Y A R R I V I N G S O U R C E W A V E L E T T I M E O F F S E T F O R S C P T S C 64 Depth (m) 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 Estimated Offset of Secondary Arriving Source Wavelet (ms) 22 25 25 26 28 28 28 28 27 29 27 31 30 30 29 The estimated time offsets outlined in Table 4.7 are consistent with the physics o f a primary source wavelet traveling directly to the S C and a secondary source wavelet reflecting off the bund and then traveling directly to the S C . For example, the estimated arrival time o f the primary arriving source wavelet at depth 17 m is 132 ms. This results i n an estimated average velocity o f 130 m/s (for a radial source offset o f 2.3 m). The extra travel distance for the source wavelet to travel to the bund, reflect and then travel to the S C is approximately 3.8 m. This calculation is based upon a S H source bund offset o f 4 m, average in-situ velocity o f 130 m/s and a 2 m estimated depth o f reflection off the bund. The associated 132 travel time offset for a 3.8 m travel distance is 29 ms (i.e., 3.8 m / 130 m/s) which is identical to the value estimated in Table 4.7. A s previously stated, the P P D - W E is a general purpose algorithm for blind deconvolution that is particularly well suited to wavelet estimation when no assumptions concerning the reflectivity can be made. In addition to the previously outlined ground breaking work o f Mendel and his colleagues [42], there have been other researchers who have implemented BRE techniques such as S M C filtering in attempting to carry out blind seismic deconvolution. These techniques rely upon a priori assumptions about the reflection coefficients. For example, in the work outlined by Cheng et. al. [23] the reflectivity sequence is modeled as a Bernoulli-Gaussian white sequence process. This statistical model was first introduced by Mendel [42]. Rosec and his colleagues [46] describe a blind marine seismic deconvolution technique where the reflectivity sequence is modeled as a Gaussian mixture depending upon three parameters (high and low reflector variances and reflector density), on the wavelet impulse response and on the observation noise variance. It is required to estimate these parameters from the recorded time series before the deconvolution can be carried out. Kaaresen & Taxt [39] give a thorough outline o f the many different approaches for carrying out blind seismic deconvolution and their associated advantages and disadvantages. A major component o f the outlined techniques is a priori assumptions about the reflection sequence. In the blind seismic deconvolution technique outlined by Kaaresen and Taxt the reflection sequence is again assumed to be a Bernoulli-Gaussian distribution. They also make additional assumptions about the reflectivity sequence such as it tends to be continuous in the horizontal direction and the wavelet is common to several traces. 133 V. Conclusions This thesis has focused upon the application o f the relatively new signal processing technique referred to as sequential Monte Carlo filtering or, more commonly, particle filtering to seismic problems. Particle filtering is a form o f Bayesian recursive estimation, where it is required to solve optimal estimation problems o f physical systems that are nonlinear and/or non-stationary. In B R E , the posterior density function is estimated, so that the conditional mean estimates o f desired parameters or states can be obtained. B R E has been referred to as a complete solution [33] to the estimation problem, since the posterior density function embodies all available statistical information (i.e., prior, likelihood and evidence). In particle filtering, the required posterior density function is represented by a set o f random samples with associated weights where state estimates are computed based on these samples and weights. Section II o f this thesis outlined the mathematical background o f B R E , where the requirement o f modeling the physical problem in a state-space formulation was addressed and the fundamental Chapman-Komolgorov equation was outlined. Additional topics covered included the Kalman filter formulation, jump Markov-linear Gaussian systems, hidden Markov models, particle filtering, and Rao-Blackwellised particle filtering. Section III o f this thesis addressed the two important seismic signal processing problems o f real-time, non-linear, non-stationary passive seismic event detection and blind seismic deconvolution. A passive seismic monitoring (PSM) system is designed to acquire and analyze, i n real time, the acoustic signals collected by an array o f appropriate seismic transducers. Seismic activity is often observed in the vicinity o f underground excavations, deep open pits and quarries, around and below large reservoirs where fluids are being injected into, or removed from, permeable subsurface formations and adjacent to the sites o f large underground explosions. Section I V ( A ) outlined a B R E algorithm (PS-SEED), which utilizes a hybrid R B P F and a H M M filter for the purpose o f identifying events during passive seismic monitoring. The event- detection algorithm builds upon previous designs, where the major improvements consist o f modeling the problem as a J M L G S and the implementation o f a 134 HMM filter for quantifying the phase o f the seismic events. Simulation results were presented where source wavelets with additive Gauss-Markov background noise were analyzed and the results evaluated in terms o f the improvement in signal-to-noise ratio. There is up to an 80-fold improvement in signal-to-noise ratio after implementation o f the P S - S E E D algorithm. Although the described filter has been utilized for automating the identification o f seismic events during passive seismic investigations, the filter could be applied to many other acoustic-emission problems requiring real-time event detection. Seismic deconvolution is the most widely utilized signal enhancement technique in seismic signal processing. Ideally, by deconvolving the source wavelet from the recorded time series data, only the reflection coefficients remain. A very challenging problem in seismic deconvolution is blind deconvolution. In this case, both the source wavelet and reflection coefficients are assumed unknown. Section I V ( C ) o f this thesis outlined two novel P F algorithms for solving the B S D problem. In these formulations, the source wavelet is again modeled as an amplitude modulated sinusoid and the deconvolution is carried out by determining the seismogram's principle phase components or by sequentially and chronologically extracting the overlapping source wavelets from the seismogram under analysis. These B S D algorithms are referred to as principle phase decomposition (PPD) and principle phase decomposition - wavelet extraction ( P P D - W E ) , respectively. The P P D and P P D - W E algorithms are similar to that o f the P S - S E E D algorithm with the added complexity o f modeling and reconstructing several overlapping source wavelets. The P P D and P P D - W E algorithms utilize R B P F , J M L G S , and H M M filter formulations i n order to separate the overlapping wavelets according to their distinct phase components. Some o f the demonstrated advantages o f the P P D and P P D - W E algorithms compared to other techniques are outlined as follows: • Simple filter formulation with minimal parameter specification. • Assumption o f minimum phase source wavelet is not required. • Avoids problems associated with band-limited source wavelets as in the case o f frequency domain deconvolution. This is due to the fact the P P D algorithm is separating overlapping source wavelets and not estimating high bandwidth reflection coefficients. 135 • It is not required to incorporate applied frequency filters within the state-space formulation as would be necessary when implementing the K F S D algorithm. • Easily handles non-stationary source wavelets. • Provides for time-variant estimations o f the source wavelet. This information can then be utilized within non-blind deconvolution techniques. • The P P D - W E is a general purpose algorithm for blind deconvolution that is particularly well suited to wavelet estimation when no assumptions concerning the reflectivity can be made. • A whiteness assumption governing the reflection coefficient series is not required. • In general the P P D - W E algorithm does not like so many B S D techniques require assumptions concerning the character o f the reflectivity. There is much room for future enhancements for both the P S - S E E D and P P D algorithms. For example, the P S - S E E D algorithm could be updated, so that it not only gives real-time phase estimates but also quantifies the frequency content o f the seismic events. A future possible enhancement to the P P D - W E algorithm is to address the limitation that any errors generated during the wavelet extraction process will propagate as the seismogram is sequentially and chronologically processed. 136 References 1) Aldrige, D . F . (1990) The Berlage wavelet. Geophysics. 55, N o . 11,1508-1511. 2) Allen, R . V . (1978) Automatic Earthquake Recognition and Timing from Single Traces. Bulletin Seismic Society of America, 1521 -1532. 3) Arulampalam, M . S . , Maskell, S., Gordon, N . J . , & Clapp, T. (2002) A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking. IEEE Transactions on Signal Processing, 50, N o . 2, 174-188. 4) Akashi, A . & Kumamoto, H . (1977). Random sampling approach to state estimation in switching environments. Automatica, 13, 429-434. 5) Andrus, R . D . , Stokoe, K . H . , & Chung, R . M . (1999). Draft guidelines for evaluating liquefaction resistance using shear wave procedures. NISTIR 6277. National velocity measurements and Institute of Standards and simplified Technology, Gaithersburg, M d . 6) Bayless, J.W. & Brigham, E . D . (1970) Application of the Kalman filter to continuous signal restoration. Geophysics, 35, 2-24. 7) Baziw, E . & Ulrych, T.J. (2006). Principle Phase Decomposition - A N e w Concept in Blind Seismic Deconvolution. IEEE Transactions on Geoscience and Remote Sensing, 2, N o . 4,418-422. 8) Baziw, E. (2005) Real-time seismic signal enhancement utilizing a hybrid Rao-Blackwellised particle filter and hidden Markov model filter. IEEE Geoscience and Remote Sensing Letter, 2, N o . 4, 418-422. 137 9) Baziw, E . (2004) State-space seismic cone minimum variance deconvolution. In Proceedings of the 2 nd International Conference on Geotechnical Site Characterization (ISC-2): V o l . 1. (pp. 835-84) Porto, Portugal: Millpress Science Publishers. 10) Baziw, E . , Nedilko, B . , & Weir-Jones, I. (2004) Microseismic event detection Kalman filter: derivation o f the noise covariance matrix and automated first break determination for accurate source location estimation. Pure Applied Geophysics, 161, 303-329. 11) Baziw, E . & Weir-Jones, I. (2002) Application o f Kalman filtering techniques for microseismic event detection. Pure Applied Geophysics, 159, 449-474. 12) Baziw, E . & Ulrych, T.J. (2004) A Rao-Blackwellised type algorithm for passive seismic event detection. In Proceedings of the 8 th Annual Consortium for the Development of Specialized Seismic Techniques Technical Meeting: V o l . 1. (pp. 135164) Vancouver, Canada: C D S S T Publishing. 13) Baziw, E . (2007). Implementation o f the Principle Phase Decomposition Algorithm. To appear in the 2007 June issue of IEEE Transactions on Geoscience and Remote Sensing. 14) Baziw, E . (1993). Digital filtering techniques for interpreting seismic cone data. Journal of Geotechnical Engineering, ASCE, 119(6), 98-1018. 15) Baziw, E . , Campanella, R . G . , & Sully, J.P. (1989). Interpretation o f Seismic Cone Data Using Digital Filtering Techniques. In Proceedings of the 12 International Conference th on Soil Mechanics and Foundation Engineering, R i o de Janeiro, 13-18 A u g . A . A . Balkema, Rotterdam. 16) Baziw, E . (2002). Derivation o f Seismic Cone Interval Velocities Utilizing Forward Modeling and the Downhill Simplex Method. Can. Geotech. J. 39, 1-12. 138 17) Baziw, E . (2004). T w o and three dimensional imaging utilizing the seismic cone penetrometer. In Proceedings of the 2nd International Conference on Geotechnical Site Characterization (ISC-2), Porto, Portugal, 19-22 Sept. Millpress Science Publishers, pp. 1611-1618. 18) Baziw, E . , Tichy, J, & de Caprona, G (2000). Data Acquisition i n Seismic Cone Penetration Testing. In Proceedings of the 3rd International Symposium on Integrated Technical Approaches to Site Characterization (ITASCE), Argonne, I L , 11-14 Sept. 2000. Argonne National Laboratory, pp. 69-72. 19) Box, G . E . , Jenkins, G . M . , & Reinsel, G . C . (1994) Time series analysis: forecasting and control (3 ed.). N e w Jersey: Prentice-Hall. rd 20) Campanella, R . G . , Robertson, F . T . C . , & Gillespie, D . (1986). Seismic cone penetration test. In Proceedings of INSITU86. American Society of Civil Engineers (ASCE) Geotechnical Special Publication. N o . 6, pp. 116-130. 21) Casella, G , & Robert, C P . (1996). Rao-Blackwellisation o f sampling schemes. Biometrika, 83, N o . 1, 81-94. 22) Crump, M . D . (1974) A Kalman filter approach to the deconvolution o f seismic signals. Geophysics, 39, 1-14. 23) Cheng, Q . , Rong C . & L i T.(1996). Simultaneous Wavelet Estimation and Deconvolution o f Reflection Seismic Signals. IEEE Transactions on Geoscience and Remote Sensing Letter, 34, N o . 4, 377-384. 139 24) de Freitas, N . (2002) Rao-Blackwellised particle filtering for fault diagnosis. In the 2002 Proceedings of the. IEEE Aerospace Conference: V o l . 4. (pp. 1767-1772) B i g S k y Montana, U S A : I E E E Publishers. 25) D e l Moral, P. (1998) Measure valued processes and interacting particle systems. Application to nonlinear filtering problems. Ann. Appl. Probab., vol. 8, no. 2, 438-495. 26) Doucet, A . , Godsill, S., & C Andrieu, C . (2000). O n sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and Computing, 10, N o . 3, 197-208. 27) Doucet, A . , de Freitas, N . , & Gordon, N . (2001) editors. Sequential Monte Carlo Methods in Practice. Springer Verlag, 2001. 28) Doucet, A . , Gordon, N . J . , & Krishnamurthy, V . (2001) Particle filters for state estimation o f jump Markov linear systems. IEEE Transactions on Signal Processing, 49, N o . 3, 613-624. 29) Doucet, A . , & Andrieu, C . (2001) Iterative algorithms for state estimation o f jump Markov linear systems. IEEE Transactions on Signal Processing, 49, N o . 6, 1216-1227. 30) Finn, W . D . L . (1984). Dynamic response analysis of soils in engineering practice. In Mechanics o f engineering materials. John W i l e y & Sons Ltd., N e w Y o r k . Chapter 13. 31) Frieden, B . R . (1983) Probability, Statistical Optics, and Data Testing (1 ed). Berlin, st Germany: Springer-Verlag Press. 32) Ge, M . , & Kaiser, P . K . (1992) Interpretation o f physical status o f arrival picks for microseismic source location, Bulletin Seismic Society ofAmerica, 80, 1643-1660. 33) Gelb, A . (\91 A) Applied Optimal Estimation (4 ed.). Cambridge, Mass: M I T Press. th 140 34) Gibowicz, S.J. & Kijko, A . (1994) An Introduction to Mining Seismology, San Diego, C A : Academic Press. 35) Gustafsson, F., Gunnarsson, F., Bergman, N . , Forssell, U . , Jansson, J., Karlsson, R., & Nordlund, P. (2002) Particle filters for positioning, navigation, and tracking, IEEE Transactions on Signal Processing, 50, 425- 435. 36) Howie, J . A . & A m i n i , A . (2005). Numerical Simulation o f Seismic Cone Signals. Canadian Geotechnical Journal. V . 42, N o . 2, 574-586. 37) Hsueh, A . C . & Mendel, J . M . (1985) Minimum-variance and maximum-likelihood deconvolution for noncausal channel models. IEEE Transactions On Geoscience and Remote Sensing, 23, N o . 6. 797-808. 38) Kanasewich, E . R . (1981) Time Sequence Analysis ( 3 ed.) Edmonton, Alberta (Canada): rd University o f Alberta Press. 39) Kaaresen, K . F . & Taxt, T. (1998). Multichannel Blind Deconvolution o f Seismic Signals. Geophysics, 63, N o . 6, 2093-2107. 40) L i u , J.S. & Chen, R. Blind (1995) Deconvolution via sequential imputations. Journal of the American Statistical Association, 90, N o . 430, 567-576. 41) Mendel, J . M . (1995) Lessons in Estimation Theory for Signal Processing, Communications, and Control. Englewood Cliffs, N . J : Prentice Hall. 42) Mendel, J . M . (1983) Optimal Seismic Deconvolution an Estimation-Based Approach. New York, N Y : Academic Press. 141 43) Muhlich, M . (2003, March). Particle filters: a tutorial. Bucuresti Filter-Workshop, Institut fur Angewandte Physik, J.W. Goethe-Universitdt Frankfurt [online]. From www.uni-frankfurt.de/~muehlich/ sci/TalkBucurestiMar2004.pdf 44) Ogata, J. (1987) Discrete Process Control (1 ed.) N e w Jersey: Prentice-Hall. st 45) Papoulis, A . (1965) Probability, Random Variables, and Stochastic Process ( 1 ed.) st N e w York, N Y : M c G r a w - H i l l . 46) Rosec O., Boucher J., Nsiri B . & Chonavel T . (2003). Blind Seismic Deconvolution Using Statistical M C M C Methods. IEEE Journal of Oceanic Engineering, 2, N o . 4, 418422. 47) Sheriff, R . E . & Geldart, L . P . (1982) Exploration Seismology, V o l . 1, ( 2 nd ed.) nd ed.) Cambridge, U K : Cambridge University Press. 48) Sheriff, R . E . & Geldart, L . P . (1983) Exploration Seismology. V o l . 2. ( 2 Cambridge, U K : Cambridge University Press. 49) Talebi, S. & Boone, T.J. (1998) Source parameters o f injection-induced microseismicity. Pure Applied Geophyics, 153, 113-130. 50) Talebi, S., Ge, M . , Rochon, P. & Mottahed, P. (1994) Analysis o f induced seismicity i n a hard-rock mine in the Sudbury basin. In Proceedings of the first North American Rock Mechanics Symposium (NARMS): V o l . 1 (pp. 937-944) University o f Texas at Austin, Rotterdam: Balkema Publishers. 51) Ulrych, T.J., and Sacchi, M . D . (2005) Information-Based Inversion and Processing with Applications (1 ed.) , Amsterdam, The Netherlands: Elsevier B . V . st 142 Appendix A. Derivation of Chapman-Kolomogorov Equation The Chapman-Kolmogorov equation is derived based upon the transitional densities o f a Markov sequence and the proof o f (8) is outlined as follows [45]: Proof For any three random variables (r.v.) x i , X2, x , we have the following: 3 p(x \x ) l = j™ p(x 2 x \x ,x x 2 I x )dx )p(x 3 2 2 2 (Al) for a Markov sequence o f order one we have the following: p(x | x ,x ) x 2 Substituting xj = x , x = xk-i, and x k 2 3 = p(x 3 x = zi:k-i | x) ( ) 2 A2 into (A2) and then substituting this result in ( A l ) gives the (8). B. Derivation of the BRE Update Equation The derivation o f (9) is outlined as follows: Derivation P( \:k \ k)P( k) z P( k x \ \:k) = z x x P(Zl; ) ( > A3 k Equation (A3) is from Bayes' rule. Separating the pdf p(zi:k) 143 into p(zk, zi.-k-i) and utilizing factorized joint probability results in the following: 2 P("k \zi:k)= ' » P(Zk>Z\:k-\) » ^ = P(Zk I Zl.-t-l. k )P( \:k-\ x I k )P( k z p(t <A4) x ) x \Z\ -\)p(Z\ -\) k :k :k From Bayes' rule (A4) becomes as follows: P( k x \ z\k)= ^ Z k * Z l ~ ' ^ P(Z | Z\ -l k l X k k :k X ^ ~ ^ ~ ^ )p(Zl.yt-l )P( k ) k Z l : k l Z l ± 1 X k ^ (A5) X Canceling out terms i n (A5) and using the assumption that we have independent measurements (white measurement noise) results in the update equation (9). 3 C . Derivation of the Recursive Weights Update Equation The derivation o f (30) is outlined as follows: Derivation From (28), we have the following: k W K ^P*( Q: \zi ) >r x k Yw . N, i=l :k , , where k i=l 2 3 q(x \z\ ) l 0:k :k p(a,b|c) = p(a|b,c)-p(b|c) and p(a,b) = p(a|b)-p(b) P( Zk\Z\ k-\ > k) = x : P(Zk\ k) x 144 ~i w\= k P ( 0k \ \-k — — — j 4( 0:k\ j lk) X z x z (Af,\ (A6) Note: ~i k= P*( l-kK-k)„P( O:k\ l:k) , t . , t . t , x , W t x z (A7) x If the importance density function is chosen to factorized such that: ^ 0:k I \:k ) = 4( k\ X Z then one can augment old particles particles 0:k-\ > \:k M 0:k-l x X Z x Q. _i by l I \:k > x ~(k x a k k ( ) Z A8 \ ok-i' \k^ x x z t § 0 et n e w x .. l Q k From (A5) we have the following: p(x . \z . )= 0 k P ( Z ~' k 1 Zl:k 1 k l °' p(Z X ° * ~ ~ | Z -\ )p(Zi; -l )P( 0:k ) k X : h Z l k l ) P ( Z l : k 1 ) P ( X ° : ^ k ( > A9 x k l:k k Canceling out terms i n and using the assumption that we have independent measurements results i n the following: p(x Q . \z . )= k x P ( Z k U Q : k k ' ° P(Zk\Z\ -\) ) p ( X k X : k x U x : k - l (A10) ) :k Utilizing factorized joint probability on (A 10) gives the following: P( 0-k \z\ )= x k P ( Z k * X k ' X 0 : k ~ l ^ * * ° ~' ~ ) P(Zk I Zi.*-i ) 145 k X :k 1 Z l : k l p ( °~ X :k l * ~ ) ( A l 1) Zl:k 1 Using the fact that we have a first order Markov process results in the following: P( 0:k x I Zl.* ) = P(z I \ k )p( k x k x P( z i k k-l) x (A12) P( Q:k-\\Z\:k-\) x ; or P( 0:k x I Zl:k ) (Zk\ k x )P( k x I k-\ x x )P( 0:k-\ x (A13) I Zl;jfc-1) Substituting (A 13) and (A8) into (A7) gives the following: P( k z ~i ~i k= k- w Equation q(x\ (A14) I -\yZ ) x l k k becomes \ 1)P( (30) X X w (AU) <*( k\ O:k-V l:k-0 x if it x is z assumed that q(x l \ x Q. _i>z\. ) l k k k = (only dependent on last state and measurement). D. Derivation of the "Likelihood" Formula Derivation Consider the stochastic processes s(t) and n(t) with the arbitrary probability laws ps(s) and PN(n) at t. Let s(t) define the signal and n(t) define a white noise process. Let d(t) also be defined as a stochastic process which is the addition o f s(t) and n(t) d(t) = s(t) + n(t) 146 (A15) The joint probability density o f d(t) and s(t) is defined as follows: p (d,s) ds = p(d\s)p (s) (A16) s From the form o f (A16) it is evident that i f s is fixed, the fluctuations in d w i l l follow those in n. This is stated probabilistically as follows: p(d\s) = p(n\s) (A17) Since the noise is white and independent o f s we have the following [31]: p(d\s) = p(n\s) = p (n) = p (d-s) n 147 n (A18) E . PS-SEED F S M C D and H M M Filter Implementation A s outlined in Table 4.1, the first step o f the P S - S E E D algorithm is to initialize a bank o f K F s equal in number (N ) to the particles pre-specified. The P S - S E E D algorithm utilizes a s F S M C D (within the J M G L S formulation) to sequentially update the set o f K F s for each time increment as was outlined in Table 3.1. The two states o f the F S M C D y l l are event (p(y )) ? and no-event (p(y )). k ~ p{y | y(_i) l k k A Markov chain is uniquely defined b y the k initial distribution and the transitional probabilities at time k = 0. For the simulations presented i n Section I V . A . l . c , the initial distribution for the no-event and event states was 1 2 set to p(y ) = 0.9 and p(y ) = 0.1, respectively. The initial transitional probabilities were 0 Q specified as p(y\ \y\_ ) = 0.8, p(y\ \y\_ ) = 0.8, p(y\ \y\_ ) = b.2, and 2 1 p(y I y _ ) = 0.2. In general terms, the probability o f a source wavelet being present l k k x x x ("event") at time index k is 20%, the transitional probability from moving from the event to no-event state is 80%, and the transitional probability from moving from the no-event to event state is 20%. The P S - S E E D algorithm behaves robustly to the specification o f the initial F S M C D probabilities as long as a sufficient number o f particles are implemented. A t each time increment k the F S M C D is updated according to the following equation: p(y\) y\\y\-\ y\\y -x p(y\). yl\y\- y\\y\-i 2 k x P(A-0 (A19) p(yli) In addition to updating the pdf o f the F S M C D at each time increment k, a Monte Carlo technique is utilized to obtain realizations o f the F S M C D equal to the number o f event/noevent particles (N ). In this step, a random number generator is utilized to obtain N s s samples o f the uniform distribution If [0,1] (i = 1 to k 148 I f If [0,1] < p(y\ ) then an k event has occurred otherwise no-event is present. In terms o f the J M G L S , the system equation does not change but the measurement equation is updated based upon the event /no-event condition. For the no-event condition the measurement equation is equal to the background noise (e.g., z\ = xl^ +v ) as outlined i n Step 1, Case 1 o f Table 4.1. For the k event condition, the measurement equation given as the addition o f the background noise and source wavelet (e.g., = x l + xl k k sin(cok&. + (p ) + v ) k k as outlined i n Step 1, Case 2 of Table 4.1. is defined by the overlapping source wavelet combination (e.g., (22) for the case o f three overlapping source wavelets). Based upon the output o f the bank o f K F s , the P S - S E E D algorithm utilizes a R B P F to obtain asymptotically optimal estimate o f the modulating amplitude term o f the A M S source wavelet defined by the state vector (Steps 6 and 7 o f Table 4.1). From the output o f the R B P F filter, a H M M filter is utilized to refine the initially specified phase component o f the sinusoid within the user specified phase resolution window (Step 9 of Table 4.1). Once the P S - S E E D algorithm obtains an optimal phase estimate at time index k, (18)-(20) o f Table 3.1 are utilized to update the bank o f K F s . The time index k is then incremented and the P S - S E E D algorithm iterates to Step 4 o f Table 4.1. 149 F. P P D F S M C D and H M M Filter Implementation A s outlined in Table 4.3, the first step o f the P P D algorithm is to initialize a bank o f K F s equal in number (N ) to that defined by (65). The P P D algorithm utilizes a F S M C D (within s the J M G L S formulation) to sequentially update the set o f K F s for each time increment as was outlined i n Table 3.1. The two states o f the F S M C D 1 (p(y )) 2 and no-event (p(y ))- k k yj^ ~ p[y k \ y{_\) a r e event A Markov chain is uniquely defined by the initial distribution and the transitional probabilities at time k = 0. For the simulations presented i n Section IV.B.2.e, the initial distribution for the event and no-event states was set to 1 2 p(y^) = 0.4 and p(y^)-0.6, specified 2 as respectively. The initial transitional probabilities were p(y\ \y{_ ) = 0A, l p(y\ \ y _ ) = OA, k x p(y\ \ y\_ ) = 0.6, and x 1 P(yk I yk-l) = 0-6- I n general terms, the probability o f an overlapping source wavelet being present ("event") at time index k is 40%, the transitional probability from moving from the event to no-event state is 60%, and the transitional probability from moving from the no-event to event state is 40%. The P P D algorithm behaves robustly to the specification of the initial F S M C D probabilities as long as a sufficient number o f particles are implemented. A t each time increment k the F S M C D is updated according to (A19). In addition to updating the pdf o f the F S M C D at each time increment k, a Monte Carlo technique is utilized to obtain realizations o f the F S M C D equal to the number o f event/noevent particles ( A ^ ) . In this step, a random number generator is utilized to obtain N$ samples o f the uniform distribution If[0,1]k (i = 1 to N ). If If[0,1]\ < p(y\^then f s an event has occurred otherwise no-event is present. For the event condition the number o f K F s sequentially updated is equal to the calculated number o f combinations defined by (64). In terms o f the J M G L S , the system equation (26) does not change but the measurement equation is updated based upon the event /no-event condition. For the no-event condition the 150 measurement equation is equal to the background noise (z -xA ). k k For the event condition, the measurement equation is defined by the overlapping source wavelet combination (e.g., (63) for the case o f three overlapping source wavelets). Based upon the output o f the bank o f K F s , the P P D algorithm utilizes a R B P F to obtain asymptotically optimal estimates o f the modulating amplitude terms o f the A M S source wavelets defined by the state vector (Steps 6 and 7 o f Table 4.3). From the output o f the R B P F filter, H M M filters are utilized to refine the initially specified phase components o f the sinusoids within the user specified phase resolution window (Step 9 o f Table 4.3). A s previously stated, a H M M filter is implemented for each A M S source wavelet. For the simulations presented in Section IV.B.2.e, the initialization o f the pdf o f the A M S source wavelet phases was set to the uniform distribution, while the fixed-grid transitional p d f p(<p \<p k-\) was set quite high (e.g., 0.996) and the remaining values l l k were set to have a uniform distribution. This is due to the fact that it is expected the A M S source wavelets have a constant phase value with the H M M filters utilized to obtain optimal phase estimates within the respective phase resolution windows. For the simulations presented in Section IV.B.2.e, each H M M filter formulation had 6 (N ) phases specified for p each phase resolution window to reflect the possible phase shifts o f ±3° in 1° increments. Once the P P D obtains optimal phase estimates at time index k, (18)-(20) o f Table 3.1 are utilized to update the bank o f K F s . The time index k is then incremented and the P P D algorithm iterates to Step 4 o f Table 4.1. 151 G . Implementation of the Three State F S M C D within the PPD-WE Algorithm The three states o f the P P D - W E F S M C D present at time index k (p(y )), measurement noise k the source wavelet understudy plus measurement noise present at time index k(p(y )) and the source wavelet understudy is overlapped with other k time series data plus measurement noise at time index k(p(y )). k A s previously stated, a Markov chain is uniquely defined by the initial distribution and the transitional probabilities at time k = 0. For the simulations presented i n Section IVB.2.f, the initial distribution for the three state F S M C D was set to p(y ) = 0.5, p(y l 2 Q ) = 0.4, and p(y^) = 0.1, respectively. This initial distribution is based upon the fact that at the start o f the time series under analysis only measurement noise or the source wavelet understudy w i l l be present (ability to window in on initial seismic wavelet recordings). The initial transitional probabilities were specified as y\\y l 0 yf\y o l [0.1429 y\\yi y\\y\ y\\y\ y\\y\ y\\y\ y\\yl y\\y\ = 0.5714 L- 0 2857 0.0357 0.0909' 0.3214 °- 6429 0.0909 °- 8182 (A20) These initial transitional probabilities were defined based upon the fact that the majority o f the time series data under analysis would reflect the case where the source wavelet understudy is overlapped with other time series data plus measurement noise at time index k. It should also be noted that it is required that the columns o f (A20) add up to 1.0 according to the law of total probability. A t each time increment k the pdf o f the three state F S M C D is updated according to the following equation: 152 ~y\ P(y\) p(y ) 2 k = y\ \yU yl \yU y\ p(y\) yU The calculated probabilities p(y\) y\ \yU p(y\-0 \yU yl \yU yl \yU yl \yU yl p(y\-\) y\-\ (A21) p(y\-x) o f (A21) are sorted from lowest to highest. The Inverse Transform Method is then utilized to obtain realizations o f the three state F S M C D equal to the number o f particles, N , specified b y the investigator. In this step, a random number s generator is utilized to obtain N samples of the uniform distribution If [0,1] (i - 1 to N ). s k s The realizations o f the three state F S M C D are then calculated as follows: y\ if realizations 1 in < (A22) p{y\)+p{yl) (A22) are based upon the assumption 3 that p(y )<p(y )<p(y ). k k otherwise outlined 2 k ifP yk yl The i (yl)< uUo.iJk 2 yk U [QX] <p{y ) k In terms o f the J M G L S , the system equation (67 does not k change but the measurement equation is updated based upon the estimated value o f y' . For k the case (z yjj. =y\ the measurement = x5 ). For the case y l k k k = v | the measurement equation is set as z k For the case y' = k equation is set equal to the background noise =xl k sin((okA + (p\ ^ ) + x 5 (A23) k the measurement equation is set to (68 and states xl and x3 are k k jittered (i.e., add on U[-1.5,+1.5]) to facilitate a diversity o f particles and minimize sample impoverishment [35]. 153 H . PPD-WE Estimated Source and Residual Wavelets for S C P T Test Hole SC 64. 1) Estimated Source Wavelet at 3m 1 2 . _ Time [ms] Figure H . l . Thirty five PPD-WE primary source wavelet estimations at depth of 3.0 m. 1:. o.a;v 6 10 20 30 40 SO T i m e 60 70 SO 90 "10O [ m s ] Figure H.2. The Averaged Estimated First Arriving Source Wavelet at 3.0m. 1 ; o.a; -1 "l O 10 • ••, 20 30 40 50 Time [ms] 60 r-> TO . . . . . . SO 90 , I 10O Figure H.3. The estimated first arriving source wavelet superimposed upon residual wavelet at 3.0 m. 154 2) Estimated Source Wavelet at 4 m 1:1 Time [ms] Figure H . 4 . T h i r t y five P P D - W E p r i m a r y source wavelet estimations at depth of 4.0 m 1:: Time [ms] Figure H . 5 . T h e Averaged Estimated First A r r i v i n g Source Wavelet at 4.0m. o.s; Time [ms] Figure H . 6 . T h e estimated first arriving source wavelet superimposed upon residual wavelet. O n e can have high confidence in the results due to the fact that the P r i m a r y and Secondary source wavelets are very similar in form. 155 3) Estimated Source Wavelet at 5m Time [ms] 100 Figure H.8. The Averaged Estimated First Arriving Source Wavelet at 5.0m. 1 Time [ms] Figure H.9. The estimated first arriving source wavelet superimposed upon residual wavelet at 5.0 m. 156 4) Estimated Source Wavelet at 6m 1:1 -1 Time [ms] Figure H.10. Thirty five P P D - W E primary source wavelet estimations at depth of 6.0 m. 1:: T i m e [ms] Figure H . l l . The Averaged Estimated First Arriving Source Wavelet at 6.0m. 1-. — — °-8: t O • * • < • ;< 25 < < : > i ^ T i m e50[ms] ; j , , ! ! , , • , 75 Figure H.12. The estimated first arriving source wavelet superimposed upon residual wavelet at 6.0 m. 157 f 100 5) Estimated Source Wavelet at 7m 50 Time [ms] Figure H.13. Thirty five P P D - W E primary source wavelet estimations at depth of 7.0 m. Time [ms] Figure H.14. The Averaged Estimated First Arriving Source Wavelet at 7.0m. Time [ms] Figure H.15. The estimated first arriving source wavelet superimposed upon residual wavelet at 7.0 m. 158 6) Estimated Source Wavelet at 8m Time 100 [ms] Figure H.16. Thirty five P P D - W E primary source wavelet estimations at depth of 8.0 m. -in o-- 1 O 1 1 ' ; 1 i 25 •• •• ' •• 1 1 • 50 T i m e [ms] •• 1 ' 1 1 1 • • • : 75 Figure H.17. The Averaged Estimated First Arriving Source Wavelet at 8.0m. 1:: Time [ms] Figure H.18. The estimated first arriving source wavelet superimposed upon residual wavelet at 8.0 m. 159 1 IOO 7) Estimated Source Wavelet at 9m '1 O : : . . . 1 25 . . . . , i , 50 Time [ms] , , , , , , , , , , 75 Figure H.19. Thirty five P P D - W E primary source wavelet estimations at depth of 9.0 m. Time [ms] Figure H.20. The Averaged Estimated First Arriving Source Wavelet at 9.0m. Time [ms] Figure H.21. The estimated first arriving source wavelet superimposed upon residual wavelet at 9.0 m. 160 ( 100 8) Etimated Source Wavelet at 10m 50 Time [msj Figure FJ.22. Thirty five P P D - W E primary source wavelet estimations at depth of 10.0m. 50 Time [ms] Figure H.23. The Averaged Estimated First Arriving Source Wavelet at 10.0m. 50 Time [ms] Figure H.24 The estimated first arriving source wavelet superimposed upon residual wavelet at 10.0 m. 161 9) Estimated Source Wavelet at 11m: O 25 50 Time [ms] 75 100 Figure H.25. T h i r t y five P P D - W E primary source wavelet estimations at depth of 11.0m. 1 ;; — .8 Time [ms] Figure H.26. T h e Averaged Estimated First A r r i v i n g Source Wavelet at 11.0m. yesian Recursive Estimation 1. _ _ 0.8- -1 \ O , , , . , , ^ 25 , , , , , , 50 Time [ms] , , , . 75 . , Figure H.27. T h e estimated first arriving source wavelet superimposed upon residual wavelet at 11.0 m. 162 100 1 10) Estimated Source Wavelet at 12m Time [ms] Figure H.28. Thirty five P P D - W E primary source wavelet estimations at depth of 12.0m. 50 Time [ms] Figure H.29. The Averaged Estimated First Arriving Source Wavelet at 12.0m. Time [ms] Figure H.30. The estimated first arriving source wavelet superimposed upon residual wavelet at 12.0 m. 163 11) Estimated Source Wavelet at 13m 50 T i m e [ms] Figure H.31. Thirty five P P D - W E primary source wavelet estimations at depth of 13.0m. 50 T i m e [ms] Figure H.32. The Averaged Estimated First Arriving Source Wavelet at 13.0m. 50 T i m e [ms] Figure H.33. The estimated first arriving source wavelet superimposed upon residual wavelet at 13.0 m. 164 12) Estimated Source Wavelet at 14m 1:1 Time [ms] Figure H.34. Thirty five P P D - W E primary source wavelet estimations at depth of 14.0m O 25 50 Time [ms] 75 Figure H.35. The Averaged Estimated First Arriving Source Wavelet at 14.0m. 1:1 Time [ms] Figure H.36. The estimated first arriving source wavelet superimposed upon residual wavelet at 14.0 m. 165 100 Time [ms] Figure H.38. The Averaged Estimated First Arriving Source Wavelet at 15.0m. O 25 50 T i m e [ms] 75 100 Figure H.39. The estimated first arriving source wavelet superimposed upon residual wavelet at 15.0 m. 166 14) Estimated Source Wavelet at 16m O 25 50 T i m e [ms] 75 100 Figure H.40. Thirty five P P D - W E primary source wavelet estimations at depth of 16.0m. o3 O 25 50 T i m e [ms] 75 100 Figure H.41. The Averaged Estimated First Arriving Source Wavelet at 16.0m. 1 o o -• -i O •• • 1 1 1 1 25 > . . . > 1 50 Tim© . [ms] , , : ! 75 1 > 1 ! . Figure H.42. The estimated first arriving source wavelet superimposed upon residual wavelet at 16.0 m. 167 1 100 15) Estimated Source Wavelet at 17m 50 Time [ms] Figure H.43. Thirty five P P D - W E primary source wavelet estimations at depth of 17.0m. Time [ms] Figure H.44. The Averaged Estimated First Arriving Source Wavelet at 17.0m. 50 Time [ms] Figure H.45. The estimated first arriving source wavelet superimposed upon residual wavelet at 17.0 m. 168
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Application of Bayesian recursive estimation for seismic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Application of Bayesian recursive estimation for seismic signal processing Baziw, Erick 2007
pdf
Page Metadata
Item Metadata
Title | Application of Bayesian recursive estimation for seismic signal processing |
Creator |
Baziw, Erick |
Publisher | University of British Columbia |
Date Issued | 2007 |
Description | Bayesian recursive estimation (BRE) requires that the posterior density function be estimated so that conditional mean estimates of desired parameters or states can be obtained. BRE has been referred to as a complete solution to the estimation problem since the posterior density function embodies all available statistical information (i.e., prior, likelihood and evidence). Until recent advances in BRE, most applications required that the system and measurement equations be linear, and that the process and measurement noise be Gaussian and white. A Kalman filter, KF, (closed form solution to the BRE) could be applied to systems that met these conditions. Previous applications of the KF to solve seismic signal processing problems (e.g., deconvolution) have had very limited success and acceptability in the geophysics signal processing community due to the restrictive nature of the KF. The recently new BRE development of sequential Monte Carlo (SMC) techniques for numerically solving non-stationary and non-linear problems has generated considerable interest and active research within the last decade. This thesis focuses upon the implementation of SMC techniques (e.g., particle filtering) for solving seismic signal processing problems. All the associated filters of BRE (hidden Markov model filter, KF, particle filter, Rao-Blackwellised particle filter, and jump Markov systems) and a new and highly robust and unique model of the seismic source wavelet are implemented in two innovative algorithms for solving the important problems of passive seismic event detection and blind seismic deconvolution. A ground-breaking concept in blind seismic deconvolution referred to as principle phase decomposition (PPD) is outlined and evaluated in this thesis. The PPD technique estimates and separates overlapping source wavelets instead of estimating high bandwidth reflection coefficients. It is shown that one can then easily generate reflection coefficients from the separated source wavelets. In this thesis many advantages of the PPD are outlined. Simulated seismogram data with low signal-to-noise ratios is blindly deconvolved where non-stationary, mixed-phase, and zero-phase source wavelets are present. I believe that there are currently no existing blind seismic deconvolution techniques which could obtain comparable performance results of the PPD technique. The work in this thesis has resulted in three IEEE publications and one peer reviewed conference publication. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-01-20 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0052619 |
URI | http://hdl.handle.net/2429/30715 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2007-266762.pdf [ 7.03MB ]
- Metadata
- JSON: 831-1.0052619.json
- JSON-LD: 831-1.0052619-ld.json
- RDF/XML (Pretty): 831-1.0052619-rdf.xml
- RDF/JSON: 831-1.0052619-rdf.json
- Turtle: 831-1.0052619-turtle.txt
- N-Triples: 831-1.0052619-rdf-ntriples.txt
- Original Record: 831-1.0052619-source.json
- Full Text
- 831-1.0052619-fulltext.txt
- Citation
- 831-1.0052619.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052619/manifest