A P P L I C A T I O N O F B A Y E S I A N R E C U R S I V E E S T I M A T I O N F O R SEISMIC S I G N A L P R O C E S S I N G by E R I C K B A Z I W B . A . S c , The University of 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 O F T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF D O C T O R OF P H I L O S O P H Y 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 B R I T I S H C O L U M B I A Apr i l 2007 © Erick Baziw, 2007 Abstract 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. 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 in 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 of the K F . The recently new B R E development o f 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 S M C techniques (e.g., particle filtering) for solving seismic signal processing problems. A l l the associated filters of 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 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 P P D 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 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 of the P P D technique. The work in this thesis has resulted in three IEEE publications and one peer reviewed conference publication. n Table of Contents Abstract i i Table of Contents i i i List of Tables v List o f Figures v i List o f Abbreviations x i i 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 Model 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) PS-SEED Filter Outline 25 a) State-Space Formulation 25 b) PS-SEED 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 of 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 Bl ind Seismic Deconvolution 45 a) Amplitude Modulated Sinusoid 47 b) P P D Algorithm Outlined 49 c) PPD Kalman Filter and Jump Markov Linear Gaussian System Formulation....52 d) PPD 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 of the Wavelet to be extracted within the P P D - W E Algorithm 89 h) Utilization of the P P D - W E Algorithm within Standard Frequency Domain Deconvolution Techniques 104 i) Monte Carlo Simulations I l l j) Sensitivity Analysis 114 3) Implementation of the P P D - W E Algorithm on Real Data 118 m V . Conclusions 134 References 137 Appendix 143 A . Derivation of Chapman-Kolomogorov Equation 143 B . Derivation of the B R E Update Equation 143 C. Derivation of the Recursive Weights Update Equation 144 D . Derivation of the "Likelihood" Formula 146 E. P S - S E E D F S M C D and H M M Filter Implementation 148 F. PPD F S M C D and H M M Filter Implementation 150 G. Implementation of 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 SCPT Test Hole SC 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 3.1 K F Governing Equations for J M L G S 10 Table 3.2 H M M Filtering Algorithm 12 Table 3.3. Typical Particle Filtering Algorithm 16 Table 4.1 PS-Seed Filter Formulation 27 Table 4.2 K F Fixed-Interval Smoother Governing Equations 35 Table 4.3 P P D Filter Formulation 55 Table 4.4 P P D - W E Filter Formulation 72 Table 4.5 Interval Velocities ( P P D - W E Extracted Source Wavelets) from the Crosscorrelation Technique for SCPT SC 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 SC 64 .... 132 v List of Figures Figure 3.1. Schematic illustrating the bootstrap re-sampling methodology 19 Figure 4.1. Source wavelet embedded with varying types of Gauss-Markov background noise 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 of 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 t h order A R M A model 42 Figure 4.6. The reflection coefficients utilized to test the performance of the K F S D algorithm 422 Figure 4.7. the output after convolving A R M A wavelet with reflection coefficients shown in Fig. 4.6 43 Figure 4.8. The K F S D estimated reflection coefficients 433 Figure 4.9. Seismic data of Fig. 4.7 with Gauss-Markov noise of T c = 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 of 96° 48 Figure 4.12. Amplitude modulating term for estimating Ricker wavelet of 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 of 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 of a 50 H z sinusoid with phase of 100° at time of 104 ms 60 Figure 4.20. The output of PPD 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 of Fig. 4.17 and the estimated Berlage wavelets of 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 of PPD algorithm where the algorithm did an impressive job in extracting the dipole Berlage source wavelets 63 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 v i Figure 4.26. The output after superimposing Berlage wavelet with arrival time of 99 ms and zero-phase Ricker wavelet of Fig. 4.11 64 Figure 4.27. The output of 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 of Fig. 4.26 and the estimated source wavelets of Fig. 4.27 65 Figure 4.29. Mixed 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 units2 and time constant = 0.001 ms). .67 Figure 4.31. The addition of a 40 H z sinusoid with phase of 100° at time of 104 ms to the sinusoid shown in Fig. 4.18 67 Figure 4.32. The output of 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 of Fig. 4.30 and the estimated source wavelets of Fig. 4.32 68 Figure 4.34. The output after convolving Berlage wavelet of F ig 4.15 with reflection coefficients illustrated in Fig. 4.16 and adding Gauss-Markov measurement noise with a variance of400 units 2 and a time constant of 0.01 ms 74 Figure 4.35. The output after applying a 100 Hz, 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 of 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 in F ig 4.16 77 Figure 4.38. The estimated and true Berlage source wavelets for the second reflection coefficient illustrated in F ig 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 of Fig 4.15 with reflection coefficients illustrated in Fig . 4.43 and adding Gauss-Markov measurement noise with a variance of400 units 2 and a time constant of 0.01 ms 80 Figure 4.45. The output after applying a 100 Hz, 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 of P P D - W E algorithm on data shown in Fig. 4.44 82 vn Figure 4.47. Estimated and true Berlage source wavelets for the first reflection coefficient illustrated in Fig 4.44 83 Figure 4.48. Estimated and true time series residual for the last two reflection coefficients illustrated in F ig 4.43 convolved with the Berlage source wavelet of 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 / 3k phase estimates for first pass of 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 in Fig 4.44 86 Figure 4.53. The estimated and true Berlage source wavelets for the third reflection coefficient illustrated in Fig 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 of the filtered synthetic seismogram of Fig. 4.44 with the A M S - E sinusoid and corresponding initial zero crossing of 4.2 ms shown 90 Figure 4.59. Illustration of the output of the H M M - F E filter after processing the synthetic seismogram shown in Fig. 4.44 91 Figure 4.60. Illustration o f BSW1 with parameters/= 55 Hz, n - 2, a = 170 and (j) = 60° specified 92 Figure 4.61. Illustration of B S W 2 with parameters/= 50 Hz, n = 2, a = 170 and <j) = 60° specified 93 Figure 4.62. Illustration of B S W 3 with parameters/= 45 Hz, n = 2, a = 170 and tj) = 60° specified 93 Figure 4.63. Illustration of 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 Hz , 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 of the output of 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 of 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 in 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 of Fig. 4.74 with the time parameters t* and t' set to 17 ms and 15 ms, respectively 101 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 102 Figure 4.77. The estimated and true BSW3 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 of seismogram illustrated in Fig. 4.79 105 Figure 4.81. Superposition of 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 of Fig. 4.83 with the time series beyond 105 ms set to zero 109 Figure 4.85. The superposition of 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 I l l 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 112 Figure 4.88. The P P D - W E source wavelet estimates for ten additive Gauss-Markov noise (time constant of 0.001 ms and a variance of 0.005 units2) 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 of 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 ix Figure 4.91. The error residual between the estimated PPD Beralge wavelet and the true Berlage wavelet with arrival time of 99 ms 115 Figure 4.92. The error residual between the estimated PPD Berlage wavelet and the true Berlage wavelet with arrival time of 108 ms 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 Figure 4.94. Typical S C P T configuration. 120 Figure 4.95. Raw X axis V S P from S C P T hole SC 64 121 Figure 4.96. Raw Y axis V S P from S C P T hole SC 64 122 Figure 4.97. V S P for the processed seismic data captured at test site SC 64 for depths 2 m to 17 m. There is clear evidence of overlapping source wavelets 123 Figure 4.98. The amplitude spectrum of 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 of SCPT hole SC 64 127 Figure 4.101. Residual wavelets V S P for depths 3 m to 17 m of 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 by 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 of 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 of 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 of 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 a t5 .0m 156 Figure H.10. Thirty five P P D - W E primary source wavelet estimations at depth of 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 of 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 of 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 of 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 of 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 of 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 of 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 of 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 of 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 of 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 of 17.0 m . 168 Figure H.44. The Averaged Estimated First Arr iving 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 . A A M S A M T A M S - E A R M A B B R E B S W B S D BIS D D F F F S M C D F M D S M H H M M H M M - F E J J M L G S K K F K F S D Meaning amplitude modulated sinusoid amplitude modulating term A M S source wavelet to be extracted autoregressive moving average Bayesian recursive estimation Berlage source wavelet blind seismic deconvolution Bayesian importance sampling dominant frequency finite state Markov chain distribution Forward Modeling / Downhil l Simplex hidden Markov model H M M frequency estimation jump Markov linear Gaussian systems Kalman filter Kalman filter seismic deconvolution xn p PF particle filters P D F probability distribution function P S M passive seismic monitoring P S - S E E D passive seismic signal enhancement and event detection PPD principle phase decomposition P P D - W E PPD filter wavelet extraction P P D - W E M C PPD filter wavelet extraction Monte Carlo R R B P F Rao-Blackwellised particle filter S S/N signal to noise ratio S M C sequential Monte Carlo SIR sampling importance re-sampling S T A / L T A short term averaging / long term averaging S W E source wavelet to be extracted SC seismic cone S C P T seismic cone penetration test V V S P vertical seismic profiling V p P-wave velocity V s S-waves velocities W W L T water level technique x i i i Acknowledgments The work outlined in this thesis would not have been possible without the guidance of 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 of my research interests. Dr. Vikram Krishnamurthy, of the University of British Columbia's Electrical and Computer Engineering Department, is an expert in the design and implementation of Bayesian recursive estimation filters, and I greatly appreciate his agreement to act as one of my committee members. Dr. Krishnamurthy is responsible for introducing me to the relatively new and exciting technique of sequential Monte Carlo filtering or more commonly known as particle filtering. Dr. Michael Bostock's expertise in seismology was fundamental for a thorough review of my work. In addition, Dr. Bostock's assistance proved essential in the completion o f my thesis and its subsequent submission to the university's Department of 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 A p p l i c a t i o n o f B a y e s i a n R e c u r s i v e E s t i m a t i o n for S e i s m i c S i g n a l P rocess ing I. Introduction There is considerable interest in the engineering community in the relatively new implementation of sequential Monte Carlo (SMC) 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 of Bayesian recursive estimation (BRE) 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 of 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 of this thesis outlines the mathematical background of B R E where the requirement of 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 of Bayless, Brigham, Crump, and Mendel the Kalman filter (KF) set of equations were utilized. Kalman filtering is a particular restrictive form of B R E where certain conditions are required. These special conditions consist of the case where the measurement and process noise are zero mean independent Gaussian white noise processes, the system equation is a linear function of the state vector and process noise, the measurement is a linear function of the state vector and measurement noise, and the initial estimate of state vector has a Gaussian distribution. The focus of my research has been to structure geophysical optimal estimation problems 1 into B R E formulations, so that the relatively new and powerful PF 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 (PSM) system is an assembly of hardware and software components designed to acquire and analyze, in real time, the acoustic signals collected by an array of 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) in the state-space system and measurement equations are due to the occurrences and losses of events within the measurement noise. The Rao-Blackwellised particle filter (RBPF) 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 of 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 of the phase of 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 of the source wavelet and a wide variety of additive measurement noise models can readily be incorporated into the model. There are two types of seismic deconvolution problems. The first consist of 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 PPD algorithm is shown to have many advantages such as simple filter formulation with minimal parameter specification, conducive to blind seismic deconvolution, assumption of 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 of 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 PPD formulations, an innovative model of 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 approximation for many analytical representations of seismic source wavelets such the exponentially decaying cyclic waveform, mixed-phase 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 PPD 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 of physical problems. Application of this type of filter requires that the dynamics of the system and measurement model that relates the noisy measurements to the system state equations be describable in 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: xk= fk-\(xk-\>uk-\) p(xk\xk-l) (1) In (1), the vector ft is a non-linear function of the state vector Xk and the process or system noise «£. (1) describes a Markov process of order one. The sampled, potentially non-linear, measurement equation is given as follows: Zk=hk(xk>vk) <-> P(zk\xk) (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 of utilizing a state-space formulation in describing physical problems is three-fold. 1) Time variance of 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 in the mathematical model of 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, in 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 of earth and instrument responses. p(t): is the reflectivity of 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 of a horizontally layered, one dimensional earth. The impedance of the Ath layer of the pancake earth model is defined as where andF^ are, respectively, the density and velocity in the M i layer. The relationship between and p^, the reflection coefficient (assuming only primary reflection) of the kth layer, is sk = Pkvk (4 ) f*k = sk+\~sk sk+\ +sk (5) 5 Rearranging (5) gives £k+\ ~€k k = *i n 7=1 1-M (6) J J Therefore, it is required to estimate the reflection series, 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 of (3) is given as z(t) = S(t)*fi(t)<^Z(z) = S(z)xV(z) (y) where ^ (z) denotes the Z transform of 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 of 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 of an A R and M A process. A s outlined in Section IV(B), Mendel [42] fitted the A R M A representation of the basic seismogram into a linear state-space formulation. Mendel then applied a linear Kalman filter smoother (a particularly restrictive form of 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 of acoustic events during passive seismic monitoring. The solution to this problem also requires a non-linear approach. In order to address the issues of blind seismic deconvolution and passive seismic event detection within the state-space and B R E contexts, the relatively new techniques of sequential Monte Carlo (SMC) filtering are both reviewed and adapted in this thesis. In general terms, S M C filters are a form of Bayesian recursive estimation (BRE) where a set of particles define the possible values of 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 of particles and their corresponding weights allows one to numerically calculate the posterior density function and subsequently derive a conditional mean estimate of the variable of interest. 7 III. Mathematical background A. Bayesian Recursive Estimation In the Bayesian approach to optimal estimation, it is attempted to construct the posterior estimate of the state given all available measurements [3]. In general terms, it is desired to obtain estimates of the discretized system equation states xk based on all available measurements up to time k (denoted as zi:k) by constructing the posterior p(xk \zi.-k) and having the initial (prior) pdf of 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 of prediction and update [3]. The prediction step involves using the system equation defined by (1) to obtain the prior pdf of the state at time k via the Chapman-Kolmogorov equation that is given as follows: P(xk I z i . * - l )=\p(xk\ xk-\ )P(xk-\ I z i . * - i )<bck-\ (8) The Chapman-Kolmogorov equation is derived based upon the transitional densities of a Markov sequence. The derivation of 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 pdf is updated via Bayes' rule as follows: n / v \r )-P(Zk\Xk)P(Xk\Z\±-\) , m P( xk I *i.-jfc) ; — : ; (9) where the normalizing constant is given as follows: p(zk I zi-k-i) = \P(zk I xk )p(xk | z\:k-\ )dxk 8 Equation (4) can also be represented as follows: likelihood • prior posterior = 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 x0 has a Gaussian distribution [33]. 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 a finite state 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 T A B L E 3.1 K E G O V E R N I N G E Q U A T I O N S F O R J M L G S Eq. Description Mathematical Representation 10 Finite state Markov chain. 4 ~ 4 i i ^ J 11 System equation. xik=F(yikK-i+G(yikK-i 12 Measurement equation. 4=^4>4 + 4 13 State estimate extrapolation. xk\k-\ = F (ylk )x'k-\\k-\ 14 Error covariance extrapolation. pi\k-i=F(yik)pUk-iF(yik)T + G(yik)QUk-iG(yik)T 15 Measurement extrapolation. 4 = H(ylk 16 Innovation. 4=4-4 17 Variance of innovation. sik=H(yik)pik\k-xH(yik)T +Ri 18 Kalman gain matrix. 19 State estimate update. xl - xl + Ki A' xk\k xk\k-l +Akak 20 Error covariance update. In (11) and (12) and Uk are mutually uncorrected i.i.d Gaussian zero mean white noise processes with variances of Qk and Rk, respectively (i.e., vk ~ N(0,Rt) and uk~N(0,Qk)). Matrices F and G define the recursive nature o f the system and process noise, respectively. Matrix F J defines the relationship between the system states and the available measurements. 10 In terms of the B R E , the Kalman Filter can be viewed as the following recursive relationship: P(xk-\ I Z\.-k-\) = N(xk-l >' xk-l\k-l > pk-\\k-\) (21a) P(xk I Z\:k-\ ) = N{xk: xk\k-\ • pk\k-\) (21b) P(xk\Z\:k) = N{xk:xk\k>Pk\k) (21c) The K F has the sufficient finite dimensional statistics of 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 of states. In the H M M filter the posterior pdf is represented by the delta function approximation as follows: p(xk-\ I ) = I ^i-M-i^v^it-i - xk-l j (22) i=\ where and tvlk-l\k-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 by substituting (22) into the Chapman-Kolmogorov equation (8) and the posterior pdf 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 F I L T E R I N G A L G O R I T H M Step Description Mathematical Representation 1 Initialization (k=0) - e.g., w[ ~\/Ns ,i = l, ...,NS-initialize particle weights. 2 Prediction - predict Ns t . . \ the weights. ™k\k-\ = Z I xl-l) 3 Update-update the wj^lP(zk \ xj) weights. Wk\k = -jr, T — ( \ iwk\k-Azk\xk) 7=1 4 Obtain optimal Ns . minimum variance xk ~ S ^Arl/t-^/t & estimate of the state i = 1 vector and p _ ' / . / ^ w i \T corresponding error xk ~ L*>k\k<xk xkAxk xk) covanance. 5 Let k = k+1 & iterate to step 2. In the above equations it is required that the likelihood pdf p(ik \ x k ) and the transitional probabilities p(x'k \ x[_x) be known and specified. 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 of filters that rely upon sequential Monte Carlo (SMC) methods have been made popular within the last decade. This family of 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 of the posterior pdf. If there was perfect Monte Carlo sampling ( x'k ~ p(xk | z, :k) ) then: according to the strong law of large numbers. For non-linear, asymptotically optimal estimation problems, it is nearly impossible to sample i\omp(xk \zx.k); therefore, the weights in (22) are obtained using Bayesian importance sampling (BIS) [3]. The basic idea in BIS is to represent the posterior density by a set of 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 of p(x) is defined as follows: Ns -> oo I TV lim E [ x k I Z\:k ] Bayesian Importance Sampling (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 xl can easily be obtained1 is then defined as follows: 1 Ns (24) q(x) = —I,S(x-xl) K ) Ns ,=i The numerator and denominator of (23) are multiplied by q(x) to give the following: Y q(x) Y q(x) Ns i = x Based upon the results of (25) it is clear that the constant 7 must be given as 2 shown: y.. 1 I'P(X') Substituting (26) into (25) gives the following: ^ N f s [ x - x ' q(x!j i=i t— w N^P (xl) 1=1 q(xi ) Comparing (27) with (23) we see that the normalized weights for estimating pdf p(x) are defined as shown: 1 Thus q(x) can be estimated numerically via realizations of x'. 2 The normalization requirement for the probability density function. 14 P w{ = = - *L where iv' = ll^l (28) i=\ q(xl) (=1 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\xl )p(xl) = P(z\xl )p( xl ) _ x i (29) q(xl) p(xl) For the PF set of equations, it can be shown as outlined in Appendix C that the following recursive relationship holds for the importance weights as shown: wk=:wk-l i P(zk\xlk)p(xlk\xlk_x) (30) q(xk \xk_vzk) The posterior density function is then approximated as follows: I yv\8(xk -xlk) P(xk I Zl.-jfc) = — 1 (31) 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)). In addition, a particle degeneracy check is made as indicated in Step 4. 15 T A B L E 3.3 T Y P I C A L P A R T I C L E F I L T E R I N G A L G O R I T H M Step Description Mathematical Representation ~ l Initialization (k=0) - x[ ~ Px*W[ ~\/N s,i = 1,...,NS. generate sample K 0 K particles and set corresponding weights. 2 Update the weights by NS . the likelihood function ™\ = Wk-\P(z J x'k) & wlk = wlk / I w\ and then normalize. Obtain asymptotically NS optimal minimum xk ~ X Wkxk & variance estimate of i = 1 the state vector and „ i , i - \, i - J corresponding error . . K « covanance. Sampling Importance A-r 7 Re-sampling (SIR). N eff ~ ~NS ~ < N? Re-sample i f NEJF <NT £ ( W K ) i=l Prediction: take u\ ~ pUk & propagate the state particles. Let k = k+1 & iterate to step 2. x l k + l =Axl +Bu\,i = 1,...,NS. 16 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 of the importance weights increase over time. A simple statistic to monitor that gives an indication of the degeneracy is the effective sample size NeJf. A n approximation to Nejf is given in Step 4 of Table 3.3, where Ns defines the number of particles utilized and NT is a user specified threshold (e.g., NT = (0.6 to 0.8) Ns). A small value of Nejf indicates severe degeneracy. The standard technique to counter the degeneracy problem is to re-sample the particles utilizing a Bayesian bootstrap technique [3] i f the effective number of particles is less than Nj. The SIR methodology implemented is a Bayesian bootstrap technique. The Bayesian bootstrap technique maps the weighted random measure , wlk j on to the equally weighted random measure jx^, iVJ*j by sampling uniformly with replacement from JJ4 with j t\Ns probabilities jw>£ t = 1 . Fig . 3.1 [43] outlines the general idea behind SIR. Fig. 3.1 illustrates the initial particles with uniform weight. After a pass o f the SIS algorithm a degeneracy check is made. If Nejf < NT, then a re-sampling type algorithm is implemented. In the qualitative example shown in F ig 3.1, Nc particles are generated (with weightNJ^) from the weighted random measure jx^, wk j . Nc is approximately equal to the number of uniform weights contained within wlk (i.e., Nc ~yvlkjN~l and Nc> 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 ig 3.1 is outlined as follows [3]: SIR Algorithm: • Initialize the C D F : cx=w\ • F o r i = 2 : N s - Construct C D F : c, = c,_, + w'k End For • Start at the bottom of the C D F : i=l • Draw a starting point from the uniform distribution: w, ~ U(0,N] ) • Forj = l : N s - Move along the C D F : « y =ul+N; l(j-l) - While Uj > C; - i = i+l - End While - Assign sample: x{ = x'k - Assign weight: w{ = J V ; 1 End For In summary, the SIR technique generates JV; children such that = Ns withE[Nt] = Nsw[ and var[Nj] = Nsw[(\ - w[) . 18 t I f - 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 by selecting an importance density that minimizes the variance of the importance weights so that NeJgr is maximized. Consider the 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 of 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. • To overcome this problem, it is critical to obtain data points from importance regions o f p(x) by explicitly searching for significant regions in the target distribution. The optimal choice of the importance density that minimizes the variance of the importance weights is defined to be as follows [3]: q( xlk | xlk_x, zk )op = P(x\\ x l k l , zk ) (32) Substitution of (32) into (30) gives the following result: ™k =™k-\P(zk I xk-\) = ™k-liP(zk \xk)p(x/k \x1k_x)dx'k (33) Equation (33) is the optimal importance density because for a givenx^. ,wlktakes the same value independent from the sample drawn fromq(xk \ xk_\, zk )op. Therefore, 20 conditional on xlk_x, the variance of different ^resul t ing from different sampled x'^is zero (i.e., var^hvk j = 0). The optimal importance density is advantageous because one can compute the importance weights n^and sample - ^ i n parallel, since w^is independent of x k . Two major drawbacks of the optimal importance density function are as follows: • The need to be able to sample from p( x'k \ x'kl ,zk ) • • The need to compute the integral defined by (33) in closed form. The two cases in which the optimal importance can be used are firstly when xk is a member of a finite set. In this case the integral in (33) becomes a sum and it is subsequently possible to sample from p( xlk \ x ' k l , zk ) . The second case in which the optimal 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 of the importance density is defined as (i.e., the prior) [3],[43]: 9(xlk\xlk_l,zk) = p(xlk\xlk_l) (34) Substitution of (34) into (30) gives: *k = KiP(*k\* ik) w The prior importance density function is easy to implement, and it tends to move samples towards the region of high likelihood. Two 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 (RBPF) allows for the reduction in the number of particles when implementing B R E on non-linear systems. This is highly advantageous because, even though there is a theoretical independence of accuracy on the particle dimension, practice has shown that the number of particles needs to be quite high for high dimensional systems [35]. The R B P F combines a bank of K F s with a PF. In this case, the K F s are utilized for generating a set of particles, where the weights of the particles are calculated with a PF [24]. In essence, the posterior pdf is approximated with a recursive, stochastic mixture of Gaussians [4], [24], [26]. This type of 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 of particles are generated by firstly computing the outlined in Table 3.1) is utilized to compute a set of particles. The posterior pdf of the state vector is then calculated and subsequent asymptotically optimal estimates are obtained [24]. finite state Markov chain distribution, which is denoted Secondly, based upon the samples drawn from a bank of Kalman filters (as 22 IV. Seismic Signal Processing A . Passive Seismic Monitoring There is considerable interest in the engineering community in the real-time identification of "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 of many structures, such as metal and concrete bridges, gas and oil pipe lines, large storage tanks and aerospace vehicles. A particularly important application of acoustics emission analysis is within the field of passive seismology. A passive seismic monitoring (PSM) system is designed to acquire and analyze, in real time, the acoustic signals collected by an array of appropriate seismic transducers. Seismic activity is often observed in the vicinity of 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 of large underground explosions [34]. Extreme examples of 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 of detecting rock failures in the vicinity of underground excavations caused by the sudden release of strain energy resulting from the redistribution of stresses around openings. Various hydrocarbon production sites also benefit from seismic monitoring systems during certain phases of production. Primary or secondary extraction or the injection of material into the reservoir to enhance production can cause significant stress changes. These stress changes can result in failures of the overlying strata and the migration of 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 of 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 of hydroelectric or large irrigation reservoirs, changes in regional loading and pore pressures cause significant stress variations within the surrounding rock mass. These can induce a wide range of micro and macro-seismic events, some of 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 of the cause of seismic events, their reliable detection and placement on a common time base is of critical importance. This is because the arrival times and phase identification of the source wavelets (P- and S-wave) at various detector packages, within a three-dimensional array, provide the basis for the calculation of the location of the seismic event. Imprecision or uncertainty in arrival time and phase determination reduces the precision of the source location operation. For example, in many passive seismic monitoring situations where there is interest in the behaviour of specific geological features, or where events must be related to specific structures in a mining or hydrocarbon extraction environment, the accurate determination of event arrival times is the primary rationale for the installation of the system. In regions where the level of induced seismicity is high, and it is accompanied by significant ambient noise, it is essential that the passive seismic monitoring systems possess the capability of 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 of a large volume of data, and the delivery of 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 of 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 of 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 (PS-SEED) 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 (BRE) 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 of 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 of utilizing a R B P F , which individually weights and subsequently sums a bank of 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 of 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 of variance and time constant. This does not preclude any other types of 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 in an ad hoc manner. To allow for greater flexibility, the PS-S E E D models state x2(t) as a Gauss-Markov process similar to that outlined for the ambient noise model. More sophisticated amplitude models can be implemented, such as the formulation of 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 RBPF 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 of Kalman filters that describe the physics of 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 of 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 PS-SEED filter formulation where Parameters a(l-2)k and b(l-2) k 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 b2k are set by specifying 1/3 of the expected maximum possible variance of the amplitude of the seismic event and the corresponding time constant (i.e., correlation between samples) [11]. Appendix E provides for more detail on the implementation of the P S - S E E D algorithm. 26 T A B L E 4.1 P S - S E E D F I L T E R F O R M U L A T I O N Step Description Mathematical Representation 1 4 5 10 11 Specify and initialize J M L G S system equations and Ns. Note: The measurement noise (vk) variance Ru is set to alk x(sk)2, where (s^2 is the variance of the Gauss-Markov noise. This allows the filter to put more weight on the prior when the Gauss-Markov noise becomes more correlated from sample to sample. Initialize the prior and transitional pdf for the J M L G S . Initialize the prior and transitional pdf for the fixed-grid phase. Draw samples for finite state Markov chain. Propagate sinusoids based current time index (t =k?) and phase estimates. Utilizing (13)-(17) outlined in Table 2.1 and updated sinusoids of Step 5, propagate the system and measurement equations, calculate importance weights for particles and then update and normalize the weights. Obtain asymptotically optimal estimate of the state vector. Sampling Importance Re-sampling (SIR). Re-sample i f Nejr < NT Use xk and H M M filter equations (Table 2.2) to estimate <plk i f ylk = y2k. Utilize (18)-(20) to update the bank of K F s Let k = k+1 & iterate to step 4. System Dynamics ah-\ xk = _x2k_ 0 + 0 a2k-l uh-\ u2k-i X\k_x x2k-l bh-l 0 0 blk_x Measurement Equations Case 1 (y\k - ambient noise): 4 = xh + vk Case 2 (y2 k - ambient noise + event): = x\k + x2k sin(a>kA + <Pk) + vk p(y\),p(yl),&p(yik \yJk_x), p((pi )&p((pi\(PJ ),ij = l,..,,Np. Np = number of fixed-grid values. yik~pk\yU Sinusoid(t) = sin((Ot+ cp'k) w[=Wk_,N(zk\zik,Sik)„ i = l,...,Ns. • Ns • wk =wk/zK xk~z Axk i=l i=i 1) de Freitas [23] demonstrates that the importance weights for y\ are given by the predictive density p(zk\ z,.k_x,y\.k ) = N(zk;z\,S'k ), where Ndenotes a Gaussian distribution. 27 c) Evaluating the PS-SEED with Simulated Data This section outlines the performance of the P S - S E E D for the specific case of 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 of initial amplitude, damping factor and dominant angular frequency [11]. The simulated seismic event had a frequency of 200 Hz , initial amplitude of 160 mm/s 2 and damping factor of 79/s specified. The sampling rate was set at 0.05ms, and a total sampling time of 300 ms was specified. Fig. 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 2 /s 4 , 0.0001 ms), (1000 mm 2 /s 4 , 0.0001 ms), (4000 mm 2/s 4 , 0.0001 ms), (1000 mm 2/s 4 , 0.1 ms), (1000 mm 2 /s 4 , 1.0 ms) and (2000 mm 2 /s 4 , 10 ms), respectively. The source wavelet had an arrival time of 150 ms (f = 0° in (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 of 133 ms (f = 140°), 138.7 ms if = 90°) and 164.4 ms (f = 45°) specified, respectively. In Fig . 4.1, the units of the amplitudes of 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 of the J M L G S system equations (Case 1 and 2 of Step 1 in Table 4.1) was carried out similarly to that outlined in [11]. The initialization of the finite state Markov chain probabilities are based upon the likelihood of an event occurring, and the transitional probability of moving from a non-event to an event. For the simulations presented here, the 1 2 1 1 pdfs and the transitional pdfs were set to p(yo) = 0.9, p(yo) = 0.1, p(yk \ yk_\) = 0.8, P(yl I y[-\) = 0-8, p(yl I yl~\) = 0-2> and p(yl | y\_x) = 0.2, respectively. In general terms, the probability of an event is about 20 percent. The initialization of the pdf of 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 (e-§-> 0-996) and the remaining values were set to have a uniform distribution. There were 90 (Np) phases specified to reflect the possible phase 28 shifts of 0° to 180° in 2° increments. The number of particles (NS) specified for the P S - S E E D filter was 100. Degeneracy threshold parameter NT was set to 0.8 NS. The P S - S E E D filter was implemented on the simulated data illustrated in Fig. 4.1, with the output for state x2k (amplitude component of A M S model) illustrated in Figs. 4.2(b)-(g). Fig. 4.2(a) shows the actual modulating amplitude component of the simulated wavelets. As is illustrated in Figs. 4.2(b)-(g), there is a very impressive real-time S/N improvement after implementation of the P S - S E E D algorithm. In these test cases, the S/N is defined as the ratio of the average maximum amplitude of the signal to the maximum average amplitude of the noise. The P S - S E E D was able to increase the S/N by approximately 80-fold, 80-fold, 30-fold, 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 of the source wavelet, making real-time noise/signal separation more difficult. The H M M filter portion o f the PS-SEED algorithm also did a good job of estimating the phase shifts for the source wavelets illustrated in Fig. 4.1 as indicated by the responses shown in Fig. 4.2. Fig. 4.3 illustrates the estimated probability of an event (P{y2k | z 1 ; £ )) for the test case shown in Fig. 4.1(b). A s expected, the probability of an event increases along with the seismic amplitude estimate illustrated in Fig. 4.2(b). 29 0 50 100 150 200 250 300 Time [ms] Figure 4.1. Source wavelet embedded with varying types of Gauss-Markov background noise. 30 0 50 100 150 200 250 300 Time [ms] Figure 4.2. PS-SEED output results for estimating state x2k after processing Fig. 4.1 data. 31 1 0.8 0.6 O L 0.4 0.2 5 0 1 0 0 1 5 0 Time [ms] 200 2 5 0 3 0 0 Figure 4.3. Estimated probability of an event for the test case shown in Fig . 4.1(b). B . Seismic Deconvolution Seismic deconvolution is one of 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 of great interest to the geophysicist. In seismology, the recorded time series, z(t), is defined to be the linear convolution of the source wavelet, S(t), with the earth's reflection coefficients, fi(t), with additive measurement noise, v(t). The mathematically representation of this relationship is given as shown: z(0 = ltoM(T)S(t-T)dT + v(t) ( 3 ? ) The discrete representation of (38) is given as follows: k z(k)=Y.MO)S(k-(i-l)) + v(k), k = l,2 — N (38) i=l There are many techniques of seismic deconvolution that can be implemented so that an optimal estimate is made of the earth model. A majority of the standard seismic 32 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. Many o f these deconvolution techniques are ad hoc in nature [42], they are effected by the band-limited nature of 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 in 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 of 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 of great interest. PF 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 Kalman Filter Smoother to Seismic deconvolution Mendel [42] is one of the original researchers to carry out work in fitting the seismic convolution model into a state-space representation for the purpose of 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 of 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 of past observations plus a Gaussian random variable. The M A process is generated by a finite linear combination of past and present inputs only. In the Kalman filter seismic deconvolution (KFSD) 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 fixed-interval 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 discrete fixed-interval Kalman 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 the first-pass estimates (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-time filtered estimate. 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 F I X E D - I N T E R V A L S M O O T H E R G O V E R N I N G E Q U A T I O N S Eq. Description Mathematical Representation xk\N =xk\k-l +rk\k-\rk\N rJW = rMUf + Hj [HjPjj + R j \ l [Zj - HjXj{j_x} 39 State smoother tm = + Pi estimate. 40 Backward recursive " J\1S ~ j equation. 41 Smoothing state pp = F ^ [ 7 . K k H k \ transition matrix. k + l 42 Smoothing error Pk\N = pk\k-\ - pk\k-\sk\Npk\k-\ covariance matrix 43 Covariance matrix of rlW. SJ\N = F J + i Sj+\\N pj+i + Hj H )PJ\}-\H'j + Rj 1-1 H In (39) k = N - l , N - 2 , •••, 1 and n x 1 vector r is a n x 1 vector. In (40) j = N , N - l , •••, 1 and r(N+l \N) = 0. In (42) k = N - l , N - 2 , , 1. In (43) j = N , N - l , •••, 1 andS N + 1 {N=0. 35 b) Kalman Filter Seismic Deconvolution (KFSD) Algorithm Outline The first step in the K F S D algorithm is for the user to specify the order of the A R M A process. For example, the Z transform for a fourth order A R M A model is given as follows: V(z) = bxzA + b2z3 + b3z2 + b4zl __X(z) z 4 + a i z 3 + a2z2 + a3z^ + a4 U(z) (44) where the parameters, bi, bj, b$, and b4, define the M A process coefficients, while the parameters, ai, a2, 05, and a4 define the A R process coefficients. V(z) is the Z transform of the source wavelet, X(z) is the Z transform of the seismogram, and U(z) is the Z transform of the in-situ reflection coefficients. In (44), the output Xk+4 is estimated on the basis of xK+3, xK+2, xk+i, xk, p-k+4, Mk+3, jUk+2, and fik+i according to the following equation: xk + 4 + alxk+3 + a2xk+2 + a 3 x k + l + a 4 x k = (45) blMk+4 + b2Mk+3 +b3Mk + 2 + b4^k + \ For the 4 t h order case3, the K F S D state-space formulation is derived by introducing variable Yl(z) into (44) [42] outlined as follows: V(z) = bxz4 +b2z3 +b3z2 +b4zl Y\(z) _X(z) z 4 + a i z 3 +a2z2 + a3zX + a4YXz) U(z) (46) Equating numerator and denominator terms in (46) results in the following two expressions: xk =b\y^k+4 +b2y^k+3 +hylk+2 + b4y^k+\ (47a) 3 Can be generalized to an A R M A model of any order. 36 Mk =ylk+4 +a\y1k+3 +a2ylk+2 + ^ ylk+l + a4y\k B y choosing xlk = ylk, x2k = yh+i, x3k = ylk+2, and x4k = ylk+3, (47a) and (47b) can be written as into a state-space formulation as follows: xlk+l x2k+l x 3 / t+l x4k+l 0 1 0 0 x l k "0" 0 0 1 0 x2k + 0 0 0 0 1 x 3 k 0 ~a3 ~a2 -ai x*k_ 1 (48) where the impedance structure /ik is defined as E\jukjuT\ = Qkd(k-r). juk is a Gaussian white noise process with mean zero and variance Qk. The discrete measurement equation is given as follows: zk = [bA b2 h] xlk x2k x3k x*k (49) The computational sequence for the discrete K F S D algorithm is outlined as follows: 1. Specify the order of the A R M A process and derive coefficients (e.g., ay, a2, a3, a4, bi, b2, b3, and 64) (subsequently addressed). 2. Define the state-space matrix equation (48) and measurement equation (49). 4. Obtain the filtered estimates (e.g., xk\k, Pk\k, Kk) by implementation of equations (9)-(20) outlined in Table 3.1. 4. Obtain smoothed estimates (e.g., xk\N, rk^, Pk\jq, F^p ) by utilizing the values derived in 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 =QkGkrk+\\N (50) A s previously stated, the first step in the K F S D algorithm is for the user to determine the order of the A R M A process and subsequently derive the necessary model parameters. This portion of the K F S D analysis is referred to as system identification. The maximum-likelihood of approximating the true source wavelet with an A R M A model increases monotonically with increasing system order, while the computational cost of increasing the system order is proportional to n 3 where n is the A R M A model order [42]. In general terms, the process of 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 of 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 of the A R M A model parameters by utilizing a least squares method is demonstrated by again considering the 4 t h order (i.e., n = 4) A R M A model given in (44). If the numerator and denominator of (44) are multiplied by z4, we obtain the following: In (51), the output xk is estimated on the basis of x^-; xk.2> xks, xk.4, juk, p,k.i, fik-2, and /uk,3 according to the following equation: X(z) _ bl+b2z 1 +b3z 2+b4z (51) U(z) \ + a\Z 1 + a 2 z 2 + a 3 z 3 + « 4 Z -4 3 8 *k = - a \ X k - l - a2Xk-2 ~ a i X k - 3 ~ a 4 X k - 4 + V * + b2Mk-l + hMk-2 + b4Mk-3 (52) where xk is the estimated value of xk. The error between the estimated output value and actual output value is defined as follows: ek = xk - xk (53) Since xk depends on past data up to n sampling periods earlier, the error ek is defined only for 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 N = C N q N + e N (54) where xN = [x4 x5 ... xN], qN = f-aj, -a2, -a3, -a4, b;, b2, b3, b4], eN = [e4 e5 ...eNJ, and CN is defined as x 3 x 4 x2 x3 x l x2 x0 x l M4 M5 M3 ju4 M2 M3 Ml M2 C N -1 X / V - l xN-2 xN-3 xN-4 MN MN-l MN-2 MN-3 The least squares performance index is defined as follows: JN = X~ X s 2 k = \ s / N s N (56) l k = 4 2 The least squares method involves minimizing (56), such that the A R M A parameter values wi l l best fit the observed data. In Ogata's [36] formulation, it is assumed that the input sequence {p.kj is such that for N>4, c'NCN is nonsingular. Ogata shows that the optimal 39 estimate of qN is given as shown: IN C N C N \ cNyN (57) In (57) it is required that {jukJ is sufficiently time-varying, so that dNCN is non-singular. Equation (57) is a first best estimate (in a least squares sense) of the A R M A parameters. Ogata presents a subsequent recursive formulation for the estimate of A R M A parameters, utilizing (57) as an initial estimate. The recursive least square estimation is defined as follows: 9 N+l ~9N + KN+l LvjV+1 -CN+1xN] (58) where K N+l CNCN\ c 'N+l 1 + c N+l (59) and cN+\ = 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 {juk} is sufficiently time-varying so that CfNCN is non-singular. Initial estimates of the A R M A parameters are obtained by implementation of (57), known /ik, and the convolved output sequence yk. The initial estimates of qN 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 Algorithm with Simulated Finite Difference Data This section presents test test-bed simulation results when implementing the K F S D algorithm previously outlined. The first step in the simulation was to define a seismic source wavelet. Amin i and Howie utilized the finite difference program F L A C to model seismic source wavelets [9], [36]. Fig. 4.4 illustrates the simulated source wavelet generated by Amin 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 of 180 m/s, and a sampling rate of 0.02 ms. I ' • i • • i • • i • • I | i i| ) M , ) | i , , i . , ) n , J | i • j , i | . i | • • , i . , n , M | u 1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 Time (ms) Figure 4.4. The finite difference simulated source wavelet. The K F S D algorithm was then implemented on the data illustrated in Fig. 4.4. In deriving the necessary A R M A parameters, a 5 t h order A R and M A was utilized with a data compression ratio of 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 of the A R M A estimation of the time reversed source wavelet shown in Fig. 4.4 41 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 Fig. 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 300 Time (ms) Figure 4.6. The reflection coefficients utilized to test the performance of the K F S D algorithm. 42 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 of Tc = 0.02 ms and standard deviation of s = 0.002 was specified. The K F S D algorithm derived the output shown in Fig. 4.8. A s is illustrated in Fig. 4.8, the arrival time locations of the reflection coefficients is recovered exactly, but there is some degradation in the estimation of the amplitudes. The K F S D algorithm gives accurate estimates of the relative amplitudes of the reflection coefficients. at 0.8 0.6-0.4-0.2-Q. E < -0.2--0.4--0.6 • i • • i 1 1 i • 1 i • • i ' • i i . i i , , . , | , . , . . | . i , i i 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 Time (ms) Figure 4.8. The KFSD estimated reflection coefficients. 43 The K F S D algorithm was next tested for its ability to derive the reflection coefficients within a high noise environment. Fig. 4.9 shows the seismic data of Fig. 4.7 with Gauss-Markov noise of Tc = 0.02 ms and s =0.08 added. The K F S D algorithm derived the estimates illustrated in Fig. 4.10. A s is shown in Fig. 4.10, the arrival time locations of the reflection coefficients are recovered exactly, but the algorithm fails to determine the amplitudes. The K F S D algorithm gives accurate estimates of the relative amplitudes of 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 . . . . . . . . . . . i . . . . . . . . . . . i . . . i . . . i . . . i . . . i . . . . . . . i > , . , . . . i . . . i 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 Time (ms) Figure 4.9. Seismic data of Fig. 4.7 with Gauss-Markov noise of T C = 0.02 ms and s = 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 KFSD 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 Tc = 0.02 ms and s =0.08 is added. The K F S D algorithm provided very similar results to those illustrated in Fig. 4.10. The main disadvantage of 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 of 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 sl'"''sm (similar to a finite state Markov chain). Fortunately for digital communication problems, the blurring function (i.e., source wavelet) has very few coefficients and there are typically a very limited number of state levels describing the input signal. For example, in the simulated results given by L i u and Chen [40], the number of blurring coefficients was three, with 16 possible state levels for the input signal. This system configuration results in a total of 16 3 = 4,096 possible combinations. Unfortunately, a typical seismic source wavelet has a minimum of 30 coefficients, and the reflections coefficients would require a minimum of 40 state levels 4. This system configuration results in a total of 4 0 3 0 = 10 4 8 combinations. In addition, the input data communication signals are characteristically highly variable, unlike the sparse reflection coefficients typically encountered in seismology. The number of combinations is indicative of the required particles for solving the blind deconvolution problem, and 10 4 8 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 of 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 (PPD-WE) 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 of several source wavelets of differing arrival times. 4 Range of values of+0.2 to -0.2 in steps of 0.01. 46 a) Amplitude Modulated Sinusoid Similar to the P S - S E E D algorithm, a fundamental and unique component of the PPD 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 of the source wavelet has proven to be highly robust and accurate in the case of passive event detection, where it is desired to identify non-stationary 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 of 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 of the A M S model. The Ricker wavelet is mathematical represented in the time-domain as [47]: / ^ ^ ^ ( l - S ^ v i O - ^ o ) 2 ) ^ - ^ ^ 0 ^ , <>t0 (61) where Ao = wavelet maximum amplitude (centered between two flanking lobes), VM = dominant or peak frequency of the Ricker wavelet, and t0 = wavelet arrival time of maximum amplitude. Fig. 4.11 illustrates a zero-phase Ricker wavelet, with the following parameters specified, Ao = 160 units, VM = 50 Hz , and t0 = 120. Superimposed upon the Ricker wavelet is a scaled 50 H z sinusoid with a phase of 96°. If the sinusoid shown in Fig . 4.11 (unsealed) is multiplied by the P S - S E E D numerically estimated amplitude modulating term shown in Fig. 4.12, then the source wavelet illustrated in Fig. 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 by 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 140 150 160 170 110 120 130 Time (ms) Figure 4.12. Amplitude modulating term for estimating the Ricker wavelet of Fig. 4.11. 180 48 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 of the P P D algorithm has similarities to the passive seismic signal enhancement and event detection (PS-SEED) fdter previously outlined, with the added complexity of modeling and reconstructing several overlapping source wavelets, as opposed to identifying a single source wavelet ("event") in real-time. In the PPD 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 PPD 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 PPD algorithm also provides for time-variant estimations of the source wavelet. This information could be utilized within non-blind (known source wavelet) deconvolution techniques. The concept of decomposing a seismogram into overlapping A M S source wavelets, based upon inherent phase differences, is illustrated in Fig. 4.14. Fig. 4.14 shows three sinusoidal wavelets that have an identical dominant frequency and varying phase components of 49 0°, 110°, and 250°. Three source wavelets are generated by modulating the amplitudes o f the sinusoidal wavelets (outlined in bold), as is shown in F ig 4.14. The seismogram can be generated and decomposed by superimposing the A M S source wavelets and separating the A M S wavelets, respectively. A Y A A A A v/ A ^ / u f V A / V Y l\ n y /7 SAA<\ A A A: A A A A A A-/ A \/ A 7 \ / V'/- \ Figure 4.14. An illustration of the convolution operation as the summation of A M S source wavelets. Similar to the PS-SEED filter, the PPD algorithm implements a R B K F that individually weights and subsequently sums a bank of 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 of 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 of 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 of Kalman filters required and the associated definition of the system and measurement equations must reflect that possibility of 50 overlapping source wavelets. In shallow and deep reflection seismology, it is reasonable to assume there wi l l be a maximum of 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 of the source wavelet. For shallow investigations (20 m to 100 m), it is feasible to assume a dominant compression source wave frequency (/) of approximately 100 H z (T = l/f= 10 ms) and a corresponding velocity (Vp) o f 1500 m/s. Based upon these parameters, 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 of 4.75 m is 5 ms (i.e., At = 2*4.75/V). If the period (10 ms) of the source wavelet is divided by the two way travel time, a ratio of 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 of the source wavelet is equal to its period. A similar calculation can be made for a deep reflection (500 m to 5 Km) survey, where Vp « 3000 m/s to 5000 m/s and / ~ 50 Hz . If it is assumed that there is a maximum of three overlapping source wavelets for an event condition, then there are a total of 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 of 200 event/no-event K F s (particles) for each time increment, and it is assumed that there is a 40 percent chance of an event, then the total number of required Kfs (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 of 10 4 8 without the requirement of modeling the reflection coefficients as discrete state levels. 51 In addition to utilizing a R B P F , and similar to the P S - S E E D algorithm, the PPD algorithm implements a set of H M M filters in order to refine the initially specified phase components of the A M S source wavelets. For example, i f an approximate phase of 100° is specified with a fixed-grid phase interval resolution of ±15°, then the H M M filter wi l l adjust the 100° phase within the -85° to 115° phase interval window to give an optimal estimate. c) P P D Kalman Filter and Jump Markov Linear 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 of possible overlapping source wavelets. For example, (62) illustrates the required K F system equation for the case of a maximum of three overlapping source wavelets. The state vector represents the amplitude component of the three source wavelets and the additive background noise (state x4k) that is modeled as a Gauss-Markov process with defining parameters of variance and time constant. A s previously stated for the P S - S E E D algorithm, this does not preclude any other types of 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 of 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 PPD measurement equation for the B S D problem when it is assumed that a maximum of three source wavelets are overlapped at time index, k, for the event condition. 52 x2k+\ x^k+i x4k+\_ a\ 0 0 0 xlk 0 0 0 " u\k 0 a,2 0 0 x2k 0 h 0 0 u2k + 0 0 a 3 0 x3k 0 0 h 0 u3k 0 0 0 an_ x4k_ 0 0 0 bn_ _u4k (62) zk = x^k sin{cokA + (p\k ) + x2k sin(a>kA + (p2^ ) + x3k sinicokA + (p3k ) + x4k (63) In (63), co represents the dominant frequency of 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. As 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 of three source wavelets. The PPD algorithm also synthesizes measurements for the no-event condition. In these cases, synthetic measurements of very low value are assigned to the respective source wavelets. This removes the residual influence of source wavelets that do not contribute to the estimated measurement equation described by (68). d) PPD Parameter Specification and Algorithm Implementation The first step in 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 of 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 PPD filter without the requirement of specifying all the initial phase estimates. This PPD 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 of the current source wavelet to be extracted. The final important parameter to specify within the P P D algorithm is the number of Kalman filters or particles required. A s previously outlined, in the R B P F implemented for the P P D algorithm, a set of particles (N;s) that define the event/no-event5 condition are initially generated by computing a F S M C D . If an event has occurred, then a set of additional K F s or particles are generated, equal to the total number of overlapping source wavelet combinations, which is calculated as follows: N -N • 1 i * NSW'- (64) ^OSC - A \ S T F + 1 + 2, — — i=2 i'-(Nsw-i)! In (64), Nose = the number of overlapping source wavelet combinations and Nsw = the number of source wavelets. The total number of particles (Ns) and corresponding Kalman filters for the PPD algorithm is calculated as follows: / , / \ T J Z 7 W . a r / (65) Ns =N0SCxPExN/s +(\-PE)xN/s In (65), P E = the probability of 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 of 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 P P D F I L T E R F O R M U L A T I O N Mathematical Representation Step Description 1 Specify the number of 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. 2 Calculate Ns and initialize the prior and transitional pdf for the event/no-event J M L G S . 3 For each phase interval, initialize the prior and transitional pdf for the fixed-grid phases. System Dynamics See (57) for the case of 3 source wavelets Measurement Equations See (58) for the case of 3 overlapping source wavelets. p(y l0).p(yZ).&p(yk\yJk-l), ij=l,2 p( (p])&p( <p\ | q>{ ), i = 1, ...,NSw, and j = l,...,Np], whereNpj = number of fixed-grid value (i.e., phase interval resolution). 4 At each time index, draw samples for finite i _ p l i J state Markov chain for the event/no-event k ^ k conditions. If event is present (i.e., y'k = ylk ) generate additional particles (i.e., Kalman filters) based upon the calculated combinations. 5 Propagate sinusoids based current time Sinusoid(1..3) = sin(cot+f (1..3)k) index (t =k?) and phase estimates. 55 T A B L E 4.3 ( C O N T I N U E D ) PPD F I L T E R F O R M U L A T I O N Step Description Mathematical Representation wik=wik_,N(zk\Zik,Sik),, i 1,...,NS. 8 7 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 NEJR < NT N. i=\ 1 9 Use xk and H M M filter equations (Table 3.2)toestimate/f/..j>) /t. 10 Utilize (18)-(20) to update the bank of 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(zk \ z, (t_1,y\.k ) = N(zk;z'k,Sk ), whereNdenotes a Gaussian 56 e) Evaluating the P P D Algorithm with Simulated Data The deconvolution of 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 of the P P D was assessed by deconvolving seismograms that contained mixed phased Berlage wavelets [1]. The Berlage wavelet is defined as follows: 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 of approximately 80 units, and the sampling rate was 0.05 ms. The Berlage wavelet of Fig. 4.15 was then convolved with the closely spaced reflection coefficients shown in Fig. 4.16 to give the output illustrated in Fig . 4.17. The synthetic seismogram shown in Fig. 4.17 and all subsequently illustrated synthetic seismograms in Section IV.B.2.e have additive Gauss-Markov noise with a variance of 50 units 2 and a time constant of 0.001 ms. w(t) = AH(t)tne n -at cos( 2/ft + (/)), (66) 100; 60; S T 0 10 20 30 40 50 60 70 Time (ms) Figure 4.15. A mixed phased Berlage wavelet. 80 90 100 110 120 130 140 57 140 160 180 200 80 100 120 Time (ms) Figure 4.16. Reflection coefficients convolved with Berlage wavelet shown in Fig. 4.15. 220 0 20 40 60 80 100 120 140 160 180 200 220 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 Hz , eighth order, zero-phase shift, low-pass Butterworth frequency filter [31] is applied so that the signal-to-noise ratio of 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 of 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 wi l l remove information associated with the high bandwidth 58 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 in a relatively complicated filter formulation becoming more complicated and unstable. In order to demonstrate the general concept of 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 of 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 of 50 H z and corresponding phase of 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 of Fig. 4.17, a dominant frequency of 50 H z is specified along with an initial phase component of 190°. Time (ms) Figure 4.18. Estimating the dominant frequency and initial phase component interactively. 59 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 of 50 H z and a corresponding phase of 190°. For these reasons, an additional scaled 50 H z sinusoid with a phase of 100° is interactively introduced and summed with the sinusoid of 190° at 104 ms. A s is shown in Fig. 4.19, the two sinusoids add such that they nearly have identical oscillations as the synthetic seismogram; therefore, the PPD algorithm is initialized with a dominant frequency of 50 H z and two phases of 190° and 100° and corresponding phase interval resolution of ±3° (i.e., 187° to 193° and 97° to 103°, respectively). Figure 4.20 illustrates the output of the P P D algorithm, where the two Berlage wavelets have been separated exactly. This is demonstrated by the error residual illustrated in Fig. 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 of 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. 140 160 180 200 80 100 120 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. 220 61 The PPD 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 in Fig. 4.22 to give the output (with additive Gauss-Markov noise) as illustrated in Fig . 4.24. 1 | : 1 o> 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 and adding the Gauss-Markov measurement noise. 62 Figure 4.24 illustrates the output of 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 of Fig. 4.23 and the estimated Berlage wavelets of Fig. 4.24. u u < . . . . i , . , , , 0 20 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 220 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 The next two sets of 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 of Fig. 4.11 is superimposed upon the Berlage wavelet illustrated in Fig. 4.15 (with an arrival time of 99 ms) to give the output (with additive measurement noise) shown in Fig. 4.26. Figure 4.27 illustrates the output of 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 of Fig. 4.26 and the estimated wavelets of Fig. 4.27. 100 120 Time (ms) 140 160 180 200 220 Figure 4.26. The output after superimposing a Berlage wavelet with an arrival time of 99 ms and the zero-phase Ricker wavelet of Fig. 4.11. 64 180 | . . I I . | . I I : i | I I I . I | I I I I I | I I I I I | | I ! I I , i , . I , , , I I | i 1 | , | , I I | | | | , , | 0 20 40 60 80 100 120 140 160 180 200 220 Time (ms) 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 of the PPD algorithms to deconvolve two Berlage wavelets, where the first wavelet has a frequency of 50 H z and the second wavelet has a frequency of 40 Hz . This test illustrates the effects of 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 of 104 ms. This source wavelet was then superimposed upon the Berlage wavelet of Fig. 4.15 with arrival time of 99 ms. Figure 4.30 illustrates the result of 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. As outlined in Fig. 4.18, the first Berlage source wavelet was determined to have a dominant frequency of 50 H z and corresponding phase of 190°. Figure 4.31 illustrates an additional scaled 40 H z sinusoid with a phase of 100°, which is interactively introduced and summed with the 50 H z sinusoid of 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 of 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 units2 and time constant = 0.001 ms). 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 140 160 180 100 120 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. 220 68 f) Proposed Methodology in Processing Seismograms with the P P D Algorithm A s previously outlined, the main challenge in implementing the PPD algorithm is the specification of the principle phase components of 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 of pre-specifying the principle phase components. The only phase which is required to be specified is that of 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 PPD filter wavelet extraction (PPD-WE) 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 of calculating the total number of overlapping source wavelet combinations (i.e., (59)). Appendix G provides for more detail on the utilization of 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: x l £ + l 1 A 0 0 0 x\k 0 0 0 0 0 " 0 x2k+l 0 a2 0 0 0 x2k 0 b2 0 0 0 u2 k x2k+l 0 0 1 A 0 x3k + 0 0 0 0 0 0 (67) x4k+\ 0 0 0 0 x4k 0 0 0 V 0 u4k x5k+l_ 0 0 0 0 an_ _x5k_ 0 0 0 unk 69 The P P D - W E algorithm models the amplitude modulating terms ( A M T ) (states xlk and x3k) of the A M S source wavelet with a first order Taylor series approximation 6. The rate of change terms (states x2k and x4k) are then approximated as Gauss-Markov processes. This allows for considerable flexibility and controllability when assigning prior to the amplitude modulating terms (e.g., smoothness). The time constant terms of states x2k and x4k are important and fundamental parameters within the P P D - W E algorithm. It is desired that states x2k and x4k result in a smooth trajectory of the amplitude modulation terms of the 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 of the user specified sinusoid. The P P D - W E applies a linear range of possible time constant terms equal to the number of specified particles (i.e., Kalman Filters) for states x2k and x4k. In the subsequent test bed simulations, the range of time constant terms varies from 0.6 ms (i.e., ai = a2 = 0.92, highly maneuverable) to 6.22 ms (i.e., aj = a2 = 0.992, sluggish). The SIR (i.e., degeneracy check) portion of the P P D - W E estimates the optimal value of the time constant for states x2k and x4k within the range specified. 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 of 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 of two source wavelets are overlapped at time index, k, for the event condition. zk = x ^ k sin{eokA + qAk) + x3^ sin(a>kA + <j?3 ^ ) + x5 ^ (68) In (68), co again represents the dominant frequency of 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 6 Please note that the amplitude modulating terms defined by (62) and (67) are forced to be positive. 70 phase of 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 of phases is generally refined based upon the initial values estimated from the first pass of data. For example, i f after the first pass of processing the seismogram it is found that the H M M estimates of /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 of the P P D - W E algorithm are the simplicity of implementation, minimal parameter specification and there is no theoretical limit on the number of overlapping source wavelets. Table 4.4 outlines the P P D - W E filter formulation. A disadvantage of the P P D - W E algorithm is that any errors generate during the wavelet extraction process wi l l propagate as the seismogram is sequentially and chronologically processed. 71 T A B L E 4.4 P P D - W E F I L T E R F O R M U L A T I O N Step Description Mathematical Representation 1 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 x2k and x4k- Based upon these parameter formulate the J M L G S system and measurement equations. System Dynamics See (62). Measurement Equations See (63) for the case of 3 overlapping source wavelets. 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. A t each time index, draw Ns samples for 3 state F S M C D . Update Kalman filter measurement equations based upon y\ as outlined in Appendix F ((A21) and (A22). If y\ =y\, then states xlk and x3k are jittered [25] (i.e., add on U[-1.5,+1.5]) to facilitate a diversity of particles. P(y\). P(vl), P(y\h & P(A I y{-\)> P(<p\ )&p(<p\ \ <p{ ),i = l,2 andy = P ...,Npi, where NPJ = number of fixed-grid value (i.e., phase interval resolution). y'k Propagate sinusoids based current time index (t =k?) and phase estimates. Sinusoid(l ,3) = sin(cot+f (l,3)k ) 72 T A B L E 4.4 ( C O N T I N U E D ) P P D - W E F I L T E R F O R M U L A T I O N Step Description Mathematical Representation Util izing (13)-(17) outlined in Table 2.1 and = wik_]N(zk \zik,Sik),, i updated sinusoids of Step 5, propagate the j ^ system and measurement equations, ' " " s ' N calculate importance weights1 for particles w\ = w | / f w | and then update and normalize the weights. i=l Obtain asymptotically optimal estimate of Ns . . the state vector (i.e., estimate source wavelet xk ~ S ^kx\ to be extracted). 1 = 1 Sampling Importance Re-sampling (SIR). - 1 Re-sample i f Neff<NT ^eff ~ T V i=i 9 Use xk and H M M filter equations (Table 2.2) to estimate/(l,3)k. 10 Utilize (18)-(20) to update the bank of 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(zk \ zX:k_x,y\.k) = N(zk;z'k,S'k), whereNdenotes a Gaussian 73 The implementation of the P P D - W E algorithm is outlined in more detail by considering the following analysis of synthetic seismograms. In the first test case, the mixed phased Berlage wavelet illustrated in Fig. 4.15 is convolved with the reflection coefficients outlined in F ig 4.16 to give the output illustrated in Fig. 4.34. The synthetic seismogram shown in Fig . 4.34 is similar to that illustrated in Fig. 4.17, but in this case the additive Gauss-Markov noise has a significantly higher variance of 400 units 2 and a lower time constant of 0.01 ms (i.e., significantly higher correlation in measurement noise) which results in 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 units2 and a time constant of 0.01 ms. A 100 Hz , eighth order, zero-phase shift, low-pass Butterworth frequency filter is then applied to the synthetic seismogram shown in F ig 4.34 to give the output illustrated in F ig 4.35. Figure 4.35 shows the filtered time series data superimposed upon the raw synthetic seismogram illustrated in F ig 4.34. In Fig 4.34 the initial 20 ms of 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 of 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 - Filtered Data -AMS-ESinusoid(105°Phase) \ Figure 4.35. The output after applying a 100 Hz, eighth order, zero-phase shift, low-pass Butterworth frequency filter to 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 x2k and x4k was 90000 (units/s)2. There were a total of Ns = 500 particles (Kalman filters) specified. The first pass of the P P D - W E algorithm implemented on the synthetic seismogram shown in Fig. 4.34 resulted in the ini t ia l /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 3k; therefore, the P P D - W E was re-executed with this phase restriction specified for/3 k . 30-" | | i i • I | I I ' 1 . | I . 1 I • | I . . i I | i 1 . \ I | I I I I I j . i ) , I , I , , 1 , j I I I I I | 0 10 20 30 40 50 60 70 80 90 100 Time [ms] Figure 4.36. The initial HMM/3 k phase estimates for first pass of PPD-WE algorithm on data shown in Fig. 4.35. 75 The second pass of 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 in F ig 4.16. Figure 4.38 shows the estimated and true Berlage source wavelets for the second reflection coefficient illustrated in F ig 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 ig 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 of the corresponding A M S s . The P P D - W E algorithm obtained an optimal estimate time constant estimate of 0.96 ms (i.e., a/ = a2 = 0.949) for states x2/( and x4k. One can readily synthesize reflection coefficients from the P P D - W E estimated source wavelets. A simplistic approached consist of 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 of the reflection coefficients. For example, Fig. 4.40 illustrates the normalized cross-correlation function for estimated source wavelets shown in Figs. 4.37 and 4.38. The peak value of the cross-correlation function in Fig . 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 of the reflection coefficient is estimated by averaging the dominant amplitude ratios of the two estimated A M S s . F ig 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 in Figs. 4.37 and 4.38. 76 100 I j . I I I • | i I I i I | i I I i I | 1 I I ! I | ! I I i I , I I , , , j : I i 1 I [ I I ! I I ,' 0 10 20 30 40 50 60 70 80 90 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 i . , r , , , . , i , , , , , , , i , / 0 5 10 15 20 25 30 35 40 45 50 Time [ms] 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 1 0.8 E < 0.6 0.4 0.2A 0 M T T T T T T T T T i l i i l i i l ! I ! i i l i l | ! i ! i l j I i I ! ! | 1 i ! 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 Time [ms] 40 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 of 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 ' k = y \ , y ' k = y \ and y l k = y \ , respectively, as outlined in Appendix G . 79 In the second test bed analysis of the P P D - W E algorithm, the mixed phased Berlage wavelet illustrated in Fig. 4.15 is convolved with the reflection coefficients outlined in F ig 4.43 to give the output illustrated in Fig. 4.44. The synthetic seismogram shown in Fig. 4.44 has additive Gauss-Markov noise with a variance of400 units2 and a time constant of 0.01 ms. 1 0.8-I 0.6-0.4-^ 0.2: — 0 Q. E -0.2 < -0.4: -0.6--0.8--1- I I I I I | i I ! I i | i i I I i : i 1 i i i 0 10 20 30 40 50 60 70 80 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 Time (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 units2 and a time constant of 0.01 ms. 80 A 100 Hz , eighth order, zero-phase shift, low-pass Butterworth frequency filter is then applied to the synthetic seismogram shown in F ig 4.44 to give the output illustrated in F ig 4.45. Figure 4.45 shows the filtered time series data superimposed upon the raw synthetic seismogram illustrated in F ig 4.44. In F ig 4.45 the initial 20 ms of 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 of 105° (±3°) of the sinusoid of 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, low-pass 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 NS = 600 particles (Kalman filters) specified. The first pass of the P P D - W E algorithm implemented on the synthetic seismogram shown in Fig. 4.44 resulted in the initial /3 k phase estimates illustrated in Fig. 4.46. A s is evident from Fig. 4.46, the phases 50° to 80° are dominant for /3 k ; therefore, the P P D - W E was re-executed with this phase resolution specified for/3 k . 81 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 PPD-WE algorithm on data shown in Fig. 4.44. The second pass of the P P D - W E algorithm provided for the A M S wavelet estimates illustrated in Figs. 4.47 and 4.48. Figure 4.47 shows the estimated and true Berlage source wavelets for the first reflection coefficient illustrated in Fig 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 ig 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 ig 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 of 1.14 ms (i.e., aj = a2 = 0.957) for states x2k and x4k. Subsequent to the extraction of the first arriving Berlage source wavelet illustrated in Fig . 4.47, the P P D - W E algorithm is then applied to the estimated residual wavelet shown in Fig. 4.48. Figure 4.50 illustrate the estimated residual wavelet shown in Fig . 4.48 with an initial phase estimate of 44° (±3°) for the sinusoid of the source wavelet to be extracted and a time index t* o f 9 ms. The first pass of the P P D - W E algorithm implemented on the estimated residual wavelet shown in Fig. 4.50 resulted in the in i t ia l /3 k phase estimates illustrated in Fig. 4.51. A s is evident from Fig. 4.51, the phases 110° to 140° are dominant for f3k ; therefore, the P P D - W E was re-executed with this phase resolution specified for f3k. 82 100 0 10 20 30 40 50 60 70 80 90 100 Time [ms] ^ J E s t i m a t e d B ^ 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] ^ ^ s t i m a t e d l R e s ^ ^ 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 20 o_ £ -20: < -40' : / \ \ Vs. /J ' ^N V — t* = 9 ms 10 20 30 40 50 Time [ms] 70 80 90 100 • 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 ' 1 1 1 1 1 1 1 1 ' 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 Time [ms] Figure 4.51. 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.50. The second pass of 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 ig 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 in Fig 4.43 The P P D - W E algorithm did a good job in obtaining accurate amplitude oscillations but it had difficulty in obtaining the correct amplitude scaling. This test bed outlines the current limitation of the P P D - W E algorithm in that any errors generate during the wavelet extraction process wi 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 in Figs. 4.52 and 4.54. The P P D - W E algorithm obtained an optimal estimate time constant estimate of 1.14 ms (i.e., ai = a2 = 0.957) for states x2k and x4k. 85 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 Time [ms] 70 80 90 100 -AMT of AMS-E • AMT of Overlapping AMS Figure 4.54. The AMTs 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 Fig. 4.52. The peak value of the cross-correlation 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 of the cross-correlation function in Fig. 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 30 35 40 45 50 25 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 of the reflection coefficients are estimated by averaging the four dominant amplitude ratios of 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 ig . 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 _ 0.2-I o . 0-=3 E < -0.2 ~| -0.4' -0.6: -0.8: 18 19 20 21 22 23 24 25 26 27 Time [ms] 28 29 30 31 32 33 34 -~ Estimated Reflection Coefficients - True Reflection Coefficients Figure 4.57. The estimated and true reflection coefficients for synthetic seismogram illustrated in Fig. 4.44. 88 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) of the source wavelet to be extracted (SWE) 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 DF. 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 of possible DFs (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 of the dominant frequency. The P P D - W E filter is re-executed with the estimate of the H M M -F E DF specified. The P P D - W E measurement equation for the case when there is only the S W E present is given as zk =x^k sin{o)t + (pl k ), where t = kA (69) In (69), co represents the D F of the S W E and / h 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 ig 4.35 the D F of 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. At the first zero crossing of 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' ( 7 1 ) B y utilizing (71) one can easily calculate a corresponding phase for any D F specified. For example, in the case outlined in Fig. 4.58 a D F o f 40 H z would result in a corresponding phase of 119.5° and a D F of 70 H z result in a corresponding phase of 74.2°. 90 The implementation of the H M M - F E filter within the P P D - W E is outlined by considering the synthetic seismogram illustrated in Fig. 4.34. In this case, the D F is assumed unknown and an initial frequency range of 30 H z to 70 H z is specified with a increment resolution of 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 of 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 Fig. 4.44. 91 A s is evident from Fig 4.59, the H M M - F E filter did an impressive job on locking onto the correct D F value of 50 H z within a relatively short time frame. The D F value of 50 H z is then specified within the P P D - W E algorithm and the synthetic seismogram illustrated in 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 of 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) = 60° specified. 92 I I i I i I 1 ' ' I I I I I I I I I I ' I P I P I I I i I i 1 I I • ! | I • i I | • i i i I i I . I | I I I | . i . I | i . . ; i . ! , • . . i • . • , | 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 Time [ms] 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 and <t> = 60° specified. 93 -60-25 T T T T T i 0 5 10 15 20 30 35 40 45 50 55 60 65 Time [ms] Figure 4.63. Illustration of BSW4 with parameters f=40 Hz, n=2,a= 170 and <j> = 60° specified. ' i i • I > M 35 40 T r r r H i 70 75 80 The synthetic seismogram illustrated in F ig 4.64 was generated by 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 and summing the results. The time series data illustrated in Fig. 4.64 has additive Gauss-Markov measurement noise with variance of 100 units 2 and a 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 i i i ! 10 20 40 50 i I ' •• > • i I I i i I I I 100 110 120 130 140 70 80 Time (ms) 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 I | I I I I | I I I I | i I I I I I ! i ! j I M ; | I ! II | I II I J I I I I j I I I ! | ! f I i i ! i j ! | I 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 Time [ms] 11111111111 90 95 100 •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 Hz , eighth order, zero-phase shift, low-pass Butterworth frequency filter is then applied to the synthetic seismogram shown in F ig 4.64 to give the output illustrated in F ig 4.67. Figure 4.67 shows the filtered time series data superimposed upon the raw synthetic seismogram illustrated in F ig 4.64 (with the initial 20 ms of 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 of 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 of 0.1 Hz . 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 in Fig . 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 wi 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 of the H M M - F E filter and an initial value of m = 35 Hz is specified. A s is evident from F ig 4.68, the H M M - F E filter did an impressive job on locking onto the correct D F value of 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, low-pass 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 in Figs. 4.69 and 4.71 7. Figure 4.69 shows the estimated and true BSW1 source wavelet. Figure 4.70 illustrates the estimated residual wavelet and the actual residual wavelet (i.e., BSW1+BSW3+BSW4). 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. 7 Please note that the P P D - W E output has a 100 Hz , eighth order, zero-phase shift, low-pass Butterworth frequency filter applied. 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. The estimated residual wavelet and the actual residual wavelet (i.e., B S W 1 + B S W 3 + BSW4). Subsequent to the extraction of the first arriving Berlage source wavelet illustrated in Fig . 4.64, the P P D - W E algorithm is then applied to the estimated residual wavelet shown in Fig. 4.70. Figure 4.71 illustrates the estimated residual wavelet shown in Fig. 4.70 with the time parameters t* and t' set to 11 ms and 8.96 ms, respectively. In the implementation of 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 of 55 Hz) was specified with an increment resolution of 0.1 Hz . Figure 4.72 shows the output of the H M M - F E filter after the first pass of the H M M - F E filter and an initial value of co = 40 Hz specified. A s is evident from Fig 4.72, the H M M - F E filter locked onto the correct D F value of 50 H z within 23 ms. 98 Bayesian Recursive Estimation 60^ 401 <u -a 20^ 3 : = 0= a . E -20-<£ •40| f = 8.9 ms -60; -80: t* = 11 ms ^ 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 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 of the P P D - W E algorithm on the data shown in Fig. 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., BSW3+BSW4). 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 Figure 4.73. The estimated and true BSW2 source wavelets. | I I l i i | I I I • I | I 1 • I I | : •• i ( I | I I : I I I . . I I 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 1 1 I . . | I I , I i | | | 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 Time [ms] Estimated Residual (BSW3+BSW4) True Residiual (BSW3+BSW4) | Figure 4.74. Estimated residual wavelet and the actual residual wavelet (i.e., BSW3+BSW4). 100 Subsequent to the extraction of Berlage source wavelets BSW1 and B S W 2 , the P P D - W E algorithm is then applied to the estimated residual wavelet shown in Fig. 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 of 30 H z to 47 H z was specified with an increment resolution of 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 of approximately 50 H z (peak 1 to peak 3). Figure 4.76 shows the output of the H M M - F E filter after the first pass of the H M M - F E filter and an initial value of co = 40 Hz specified. A s is evident from Fig 4.77, the H M M - F E filter locked onto a D F value of approximately 44.5 H z within 30 ms. | . ! I i I | I i I i . j | i I I . I ; I > j I . | I I I i I | I M I I | I I I I 0 5 10 15 20 25 30 35 40 45 50 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 of the P P D - W E algorithm on the data shown in Fig . 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 w i l l propagate as the seismogram is sequentially and chronologically processed. To the best of my knowledge, the previously processed synthetic seismogram is the first example o f the blind extraction of 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 of 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 of stone columns utilizing vibro-compaction). This makes the post analysis of 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 in separating the closely spaced overlapping source wavelets. I I I i i I i I | H I , , | I I , , , | • I , I I j , , I , . 30 35 40 45 50 55 60 65 Time [ms] 70 75 80 — Estimated BSW3 • True BSW3 Figure 4.77. The estimated and true B S W 3 source wavelets. 60-40; 20; T3 =3 _ fl-"a. E < -20; -40--60-/ / V \\ 10 15 20 25 30 35 40 45 Time [ms] i " 1 50 55 60 65 70 75 80 • Estimated BSW4 — True BSW4 Figure 4.78. The estimated and true B S W 4 source wavelets. 103 h) Utilization of the P P D - W E Algorithm within Standard Frequency Domain Deconvolution Techniques A standard frequency domain methodology in estimating the reflection series, juk, is the water level technique (WLT) [51]. The W L T is outlined by considering (7), representing the Z transform by the Fourier transform (i.e., z = e10" ) and rearranging terms as follows: * r « ; = ^ (72) S(co) where x¥(co) denotes the Fourier transform of the reflection series /J(t). Theoretically, for a known source wavelet, one could simply implement (72) and calculate x¥(co). The reflection series, nk, is then estimated by taking the inverse Fourier transform of ^(co). Unfortunately, due to inaccuracies in the specification of the source wavelet, the bandlimited nature of the source wavelet and additive measurement noise, the implementation of (72) is highly unstable and inaccurate. To mitigate the previously outlined limitations, (72) is modified by firstly multiplying the numerator and denominator by the complex conjugate of S(co) (denoted as S*(co)) and secondly, by introducing an additive scalar value to the denominator which is referred to as the water level (A) [51]. Implementation of these two modifications to (72) gives Z(co)S*(co) ^Z(co)S*(co) ( ? 3 ) S(co)S*(co) + A Ps(co) + A where Ps(co) denotes the power spectrum o f the source wavelet (i.e., the Fourier transform of the autocorrelation of S(t)). In general terms, the setting of 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 and the source wavelet (i.e., Z(co)S*(co)). 104 The implementation of the W L T in conjunction with the P P D - W E technique is illustrated by considering the synthetic time series illustrated in Fig. 4.79. The simulated time series shown in this figure is a typical seismogram which one may encounter and it was generated without my prior knowledge of the source wavelet(s) or reflection series. 0 50 100 150 200 250 300 350 400 450 500 Time [ms] Figure 4.79. A typical synthetic seismogram. The amplitude spectrum of 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 Hz . 0.081 — Frequency[Hz] Figure 4.80. The amplitude spectrum of the seismogram illustrated in Fig. 4.79. 105 A 100 H z low-pass Butterworth frequency filter was then applied to the seismogram of Fig. 4.79. Figure 4.81 illustrates the filtered seismogram (black line) superimposed upon the noisy seismogram (light grey plot) of Fig. 4.79. 0.8 o.6; 0.4; • 0 2 3 0 §- -0.2; •* -0.4; -0.6! -0.8 -1-1 50 100 150 200 300 350 400 450 250 Time [ms] Figure 4.81. The superposition of filtered seismogram onto seismogram shown in Fig. 500 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 of the filtered time series was ignored due to the fact no signal was present. A s a result of the estimated frequency bandwidth of the seismogram (i.e., Fig. 4.80), the H M M - F E component of the P P D - W E algorithm had a frequency range of 40 H z to 60 H z specified. Figure 4.82 illustrates the estimated dominant frequency output of the H M M - F E . A s is evident from Fig. 4.82, the H M M - F E locks onto a dominant frequency range of approximately 49 H z to 51 H z within 5 ms from the onset of the first arriving source wavelet. 106 45 : . . • . . ) I . , , , | , . , . , , , , , , , , | , , , I , , : I , , , 0 5 10 15 20 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 Fig. 4.81. Figure 4.83 shows the P P D - W E estimated first arriving source wavelet superimposed upon the true source wavelet. As is evident from Fig. 4.83, the P P D - W E algorithm estimated the correct dominant frequency of 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 of view, decay to zero and start up again, the estimated time series beyond 105 ms is set to zero. F ig . 4.84 illustrates the estimated source wavelet of 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 of separating all possible overlapping source wavelets, the first deconvolution attempt was performed using the estimated source wavelet of Fig. 4.84 and the W L T . If there is minimal source wavelet variation within the time series, the estimated reflection coefficients wi l l be very similar in shape. Figure 4.85 illustrates the superimposition of the true reflection series onto the estimated reflection series for the synthetic seismogram shown in Fig. 4.79. For these results, an 8 t h order Butterworth low pass filter was applied to the noisy seismogram and the water level was set to 0.2% of the maximum value of the power spectrum of 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 0 -0.2 -0.4 '• if! f 1 U H I ti *i ' A nil i i i ) I - H f _1 I I I • ' _l I I l _ -1— I — I — 1 — I I I I I I I I I 1 I I I 0 50 100 150 200 250 300 350 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 (MCS) considerations in the implementation of the P P D - W E algorithm. 1) Step 4 of Table 4.4, where, i f y\ =y\ , then states xlk and x3k are jittered by means of a random number generator. 2) The response of 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 by 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. I l l The ability of P P D - W E algorithm to respond to varying noise realizations is also of concern. Figure 4.87 illustrates an example of a 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. 1:1 i , . . . , 8 i i 0 10 20 30 40 50 60 70 80 90 100 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 units2. Figure 4.88 illustrates the output of the P P D - W E algorithm for ten additive Gauss-Markov noise (time constant of 0.001 ms and a variance of 0.005 units2) realizations. A s is evident from Fig . 4.88, the P P D - W E had high repeatability when processing the ten synthetic seismograms. 112 0.8:i Time [ms] Figure 4.88. PPD - W E source wavelet estimates for ten additive Gauss-Markov noise (time constant of 0.001 ms and a variance of 0.005 units 2) realizations. The bold light 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 of 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 of 99 ms and the second Berlage wavelet has an arrival time of 108 ms. These two source wavelets have associated phases o f 190° and 20°, respectively. This results in a phase difference of 170° that is 10° from the indiscernible value o f 180°. Figure 4.90 shows 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. This is highlighted in 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 of 99 ms is shown. Figure 4.92 illustrates the error residual between the estimated PPD Beralge wavelet and the true Berlage wavelet for the source wavelet with an arrival time of 108 ms. Due to the nature of these and similar results, it is believed that the current formulation of 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 of ±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 of ±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 PPD-WE 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 of 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 of 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 of low strain (<10~6) in-situ compression (P) and shear (S) wave velocities. 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 Young's modulus. The estimation of elastic constants is essential in geotechnical foundation designs [30]. Another important use of estimated shear wave velocities in geotechnical design is in the liquefaction assessment of soils [5]. Details of 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 of 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 in 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 difference and v is the corresponding SCPT interval velocity). It is of 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 of the SC data acquisition system and seismic sensors utilized are given in [18]. For the S H wave analysis, the preprocessing of the seismic time series data captured on the X and Y axes is three-fold. 1) Apply bandpass frequency filter. 2) Rotate the X and Y axes responses onto the full waveform axis utilizing hodograms and polarization analysis [10], [17]. 3) Apply 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 SCPT configuration where a S H source is located at the outriggers of the in-situ testing vehicle. A triaxial SC 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 Fig. 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 New Zealand on February 13 t h, 2006 near the Tauranga Eastern Motorway by Perry Dril l ing Ltd. o f Tauranga New 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 of the rod strings. The source consists of a sledge hammer horizontally impacting point source steel beams located underneath the outriggers and an electrical contact trigger is utilized. Two stacked S H sources were generated for each depth increment. Fig. 4.95 illustrates raw X axis time series data captured from SCPT test site SC 64. Fig . 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. Fig. 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. Typ ica l S C P T configuration. 120 Figure 4.95. Raw X axis V S P from S C P T hole SC 64. 121 1 , , , . , , . ; , , ! 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 Time [ms] Figure 4.96. Raw Y axis VSP from SCPT hole SC 64. 122 18 20 40 BO BO 100 120 140 160 180 200 220 240 260 280 300 Time [ms] Figure 4.97. V S P for the processed seismic data captured at test site SC 64 for depths 2 m to 17 m. There is clear evidence of overlapping source wavelets. 123 In the implementation of the P P D - W E technique there are slightly differing possible source wavelet estimation realizations depending upon the specification of the initial P P D - W E filter parameters of / (zero crossing), /* (time to source wavelet overlap), and Rk, (variance o f measurement noise) and the utilization of Monte Carlo techniques to obtain realizations o f yk ~ I y{-i) an^ J i t t e r m e particles (i.e., Step 4 of 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 of 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, t2, is introduced. Unfortunately, there are some situations where the A M T of the source wavelet w i l l slowly start to diverge from the true value due to the degeneracy check. To mitigate this effect, the degeneracy check is turned off after time t2 which results in a noisier estimated source wavelet but the diversity of the particles is maintained. In the P P D - W E M C algorithm the investigator initially specifies minimum t^min (first seismogram zero crossing) and t2min parameters. These parameters are modified within each iteration of 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 t2 ~ N(0,225). In addition, the measurement noise variance is increased from an initial user specified minimum, Rmin, according to a specified increment value, Rinc (i.e., R = Rmin + Rinc). 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 of the acquired seismograms. For example, Fig. 4.98 shows the amplitude spectrum of 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 Hz . 124 0.04 T3 =5.0.02-E < 40 50 60 Frequency [Hz] Figure 4.98. The amplitude spectrum of the seismogram captured at 2 m. In the implementation of the P P D - W E algorithm on the SC 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 of the PPD-W E M C algorithm and the minimum noise variance, Rmin , was set to 0.1 with an increment increase of 0.1. This results in a measurement noise variance of 3.6 for the final run of the P P D - W E M C algorithm. 25 30 35 40 Time [ms] 45 50 55 60 Figure 4.99. Estimating tmin for the seismogram acquired at 12 m. 125 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 SCPT hole SC 64. The time series recorded at 2.0 m was assumed to contain no source multiples. Appendix H provides illustrations of 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 Fig . 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 in Table 4.5. The traces shown in Fig. 4.100 have very high calculated interval correlation coefficients as shown in Table 4.5. Table 4.6 outlines the estimated interval velocities i f the processed seismograms of F ig 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 of 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 of the elastic constant. The crosscorrelation interval arrival times derived from the processed seismograms of Fig. 4.97 are also sensitive to the time window specified for the crosscorrelation function. Fig. 4.101 illustrates the residual wavelets derived by subtracting the estimated primary source wavelets outlined in Fig. 4.100 from the processed seismograms shown in Fig . 4.97. A s is evident from Fig. 4.101 and the residual wavelet plots shown in Appendix H , the majority of 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, Fig. 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. Fig. 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 PPD-WEMC 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 I N T E R V A L 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) Correlation Coefficient 2 . 0 - 3 . 0 1 0 6 0 . 9 0 3 .0 -4 .0 1 2 0 0 . 9 6 4 . 0 - 5 . 0 9 6 0 . 9 7 5 .0 -6 .0 1 1 8 0 . 9 8 6 .0 -7 .0 141 0 . 9 6 7 .0 -8 .0 1 5 6 0 . 9 9 8 .0-9 .0 1 5 0 0 . 9 5 9 . 0 - 1 0 . 0 1 3 9 0 . 9 5 1 0 . 0 - 1 1 . 0 143 0 . 9 6 1 1 . 0 - 1 2 . 0 1 2 2 0 . 98 1 2 . 0 - 1 3 . 0 1 1 0 0 . 9 9 1 3 . 0 - 1 4 . 0 1 6 8 0 . 9 6 1 4 . 0 - 1 5 . 0 148 0 . 9 9 6 1 5 . 0 - 1 6 . 0 1 6 6 0 . 9 9 1 6 . 0 - 1 7 . 0 2 1 6 0 . 9 9 5 T A B L E 4 . 6 I N T E R V A L V E L O C I T I E S ( P R O C E S S E D S E I S M O G R A M 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 S C 6 4 Interval Depth (m) Crosscorrelation Interval Velocity Estimate (m/s) Correlation Coefficient 2 . 0 - 3 . 0 9 3 0 . 9 0 3 .0 -4 .0 1 0 4 0 . 8 6 4 . 0 - 5 . 0 9 2 0 . 9 3 5 .0 -6 .0 1 1 4 0 . 9 7 6 .0 -7 .0 121 0 .81 7 .0 -8 .0 1 6 4 0 . 9 7 8 .0 -9 .0 1 5 2 0 . 8 4 9 . 0 - 1 0 . 0 1 6 7 0 . 9 4 1 0 . 0 - 1 1 . 0 1 4 4 0 . 9 7 1 1 . 0 - 1 2 . 0 1 1 8 0 . 7 7 1 2 . 0 - 1 3 . 0 1 0 8 0 .91 1 3 . 0 - 1 4 . 0 1 5 4 0 . 8 0 1 4 . 0 - 1 5 . 0 1 6 7 0 . 9 6 1 5 . 0 - 1 6 . 0 148 0 . 8 5 1 6 . 0 - 1 7 . 0 2 8 7 0 . 7 9 1 2 8 "I • •• ! ! 1 i ! i 1 ! •• I ! 1 1 1 1 : 1 : i ! ! ; 1 , , , 1 , , , , 1 , , : , 1 25 50 75 100 125 150 175 200 225 Time [ms] Figure 4.101. Residual wavelets VSP for depths 3 m to 17 m of SCPT hole SC 64. 129 150 200 Time [ms] 250 300 350 Figure 4.102. Estimated reflection coefficients for depth 16 m utilizing the W L T . 150 200 Time [ms] 250 300 350 Figure 4.103. Estimated reflection coefficients for depth 17 m utilizing the W L T . 130 Fig. 4.104 shows the superposition of the P P D - W E estimated source wavelets outlined in Fig. 4.100 (depths 4 m to 17 m) and corrected zero time offset. Fig. 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 of the central trough. 0 25 50 Time [ms] Figure 4.104. Superposition of estimated primary 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 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 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 of dense material and the S H source for S C P T SC 64 was positioned approximately 4m from the bank of 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 SCPT SC 64 Depth (m) Estimated Offset of Secondary Arriving Source Wavelet (ms) 3.0 22 4.0 25 5.0 25 6.0 26 7.0 28 8.0 28 9.0 28 10.0 28 11.0 27 12.0 29 13.0 27 14.0 31 15.0 30 16.0 30 17.0 29 The estimated time offsets outlined in Table 4.7 are consistent with the physics of a primary source wavelet traveling directly to the SC and a secondary source wavelet reflecting off the bund and then traveling directly to the SC. For example, the estimated arrival time of the primary arriving source wavelet at depth 17 m is 132 ms. This results in an estimated average velocity of 130 m/s (for a radial source offset of 2.3 m). The extra travel distance for the source wavelet to travel to the bund, reflect and then travel to the SC is approximately 3.8 m. This calculation is based upon a S H source bund offset of 4 m, average in-situ velocity of 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 B R E 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 of 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 of physical systems that are non-linear and/or non-stationary. In B R E , the posterior density function is estimated, so that the conditional mean estimates of 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 of 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 of B R E , where the requirement of 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 of this thesis addressed the two important seismic signal processing problems of 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, in real time, the acoustic signals collected by an array of appropriate seismic transducers. Seismic activity is often observed in the vicinity of 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 of large underground explosions. Section IV(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 of identifying events during passive seismic monitoring. The event- detection algorithm builds upon previous designs, where the major improvements consist of modeling the problem as a J M L G S and the implementation of a 134 H M M filter for quantifying the phase of the seismic events. Simulation results were presented where source wavelets with additive Gauss-Markov background noise were analyzed and the results evaluated in terms of the improvement in signal-to-noise ratio. There is up to an 80-fold improvement in signal-to-noise ratio after implementation of the P S - S E E D algorithm. Although the described filter has been utilized for automating the identification of 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 IV(C) of this thesis outlined two novel PF 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 (PPD-WE), respectively. The PPD and P P D - W E algorithms are similar to that of the P S - S E E D algorithm with the added complexity of modeling and reconstructing several overlapping source wavelets. The PPD and P P D - W E algorithms utilize R B P F , J M L G S , and H M M filter formulations in order to separate the overlapping wavelets according to their distinct phase components. Some of the demonstrated advantages of 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 of minimum phase source wavelet is not required. • Avoids problems associated with band-limited source wavelets as in the case of frequency domain deconvolution. This is due to the fact the PPD 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 of 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 of 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 of 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 wi l l propagate as the seismogram is sequentially and chronologically processed. 136 References 1) Aldrige, D.F . (1990) The Berlage wavelet. Geophysics. 55, No . 11,1508-1511. 2) Al len , 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, No . 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 velocity measurements and simplified procedures. NISTIR 6277. National Institute of Standards and 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 New Concept in Bl ind Seismic Deconvolution. IEEE Transactions on Geoscience and Remote Sensing, 2, No . 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, No . 4, 418-422. 137 9) Baziw, E . (2004) State-space seismic cone minimum variance deconvolution. In Proceedings of the 2nd 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 of 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 of 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 8th Annual Consortium for the Development of Specialized Seismic Techniques Technical Meeting: V o l . 1. (pp. 135-164) Vancouver, Canada: C D S S T Publishing. 13) Baziw, E . (2007). Implementation of 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 of Seismic Cone Data Using Digital Filtering Techniques. In Proceedings of the 12th International Conference on Soil Mechanics and Foundation Engineering, Rio de Janeiro, 13-18 Aug. A . A . Balkema, Rotterdam. 16) Baziw, E . (2002). Derivation o f Seismic Cone Interval Velocities Util izing Forward Modeling and the Downhil l Simplex Method. Can. Geotech. J. 39, 1-12. 138 17) Baziw, E . (2004). Two 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 in Seismic Cone Penetration Testing. In Proceedings of the 3rd International Symposium on Integrated Technical Approaches to Site Characterization (ITASCE), Argonne, IL, 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 (3rd ed.). New Jersey: Prentice-Hall. 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. No. 6, pp. 116-130. 21) Casella, G , & Robert, C P . (1996). Rao-Blackwellisation of sampling schemes. Biometrika, 83, No. 1, 81-94. 22) Crump, M . D . (1974) A Kalman filter approach to the deconvolution of seismic signals. Geophysics, 39, 1-14. 23) Cheng, Q., Rong C. & L i T.(1996). Simultaneous Wavelet Estimation and Deconvolution of Reflection Seismic Signals. IEEE Transactions on Geoscience and Remote Sensing Letter, 34, No. 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 Sky Montana, U S A : IEEE Publishers. 25) Del 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). On sequential Monte Carlo sampling methods for Bayesian filtering. Statistics and Computing, 10, No. 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 of jump Markov linear systems. IEEE Transactions on Signal Processing, 49, No. 3, 613-624. 29) Doucet, A . , & Andrieu, C. (2001) Iterative algorithms for state estimation of jump Markov linear systems. IEEE Transactions on Signal Processing, 49, No . 6, 1216-1227. 30) Finn, W . D . L . (1984). Dynamic response analysis of soils in engineering practice. In Mechanics of engineering materials. John Wi ley & Sons Ltd., New York . Chapter 13. 31) Frieden, B .R. (1983) Probability, Statistical Optics, and Data Testing (1st ed). Berlin, Germany: Springer-Verlag Press. 32) Ge, M . , & Kaiser, P .K . (1992) Interpretation of physical status of arrival picks for microseismic source location, Bulletin Seismic Society of America, 80, 1643-1660. 33) Gelb, A . (\91 A) Applied Optimal Estimation (4th ed.). Cambridge, Mass: M I T Press. 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. & Amin i , A . (2005). Numerical Simulation of Seismic Cone Signals. Canadian Geotechnical Journal. V . 42, No. 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, No . 6. 797-808. 38) Kanasewich, E.R. (1981) Time Sequence Analysis (3 r d ed.) Edmonton, Alberta (Canada): University of Alberta Press. 39) Kaaresen, K . F . & Taxt, T. (1998). Multichannel Bl ind Deconvolution of Seismic Signals. Geophysics, 63, No. 6, 2093-2107. 40) L iu , J.S. & Chen, R. Blind (1995) Deconvolution via sequential imputations. Journal of the American Statistical Association, 90, No. 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 s t ed.) New Jersey: Prentice-Hall. 45) Papoulis, A . (1965) Probability, Random Variables, and Stochastic Process (1 s t ed.) New York, N Y : McGraw-Hi l l . 46) Rosec O., Boucher J., Nsir i B . & Chonavel T. (2003). Bl ind Seismic Deconvolution Using Statistical M C M C Methods. IEEE Journal of Oceanic Engineering, 2, No. 4, 418-422. 47) Sheriff, R .E . & Geldart, L .P . (1982) Exploration Seismology, V o l . 1, (2 n d ed.) Cambridge, U K : Cambridge University Press. 48) Sheriff, R . E . & Geldart, L .P . (1983) Exploration Seismology. V o l . 2. (2 n d ed.) Cambridge, U K : Cambridge University Press. 49) Talebi, S. & Boone, T.J. (1998) Source parameters of injection-induced microseismicity. Pure Applied Geophyics, 153, 113-130. 50) Talebi, S., Ge, M . , Rochon, P. & Mottahed, P. (1994) Analysis of induced seismicity in 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 of Texas at Austin, Rotterdam: Balkema Publishers. 51) Ulrych, T.J., and Sacchi, M . D . (2005) Information-Based Inversion and Processing with Applications (1st ed.) , Amsterdam, The Netherlands: Elsevier B . V . 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 of (8) is outlined as follows [45]: Proof For any three random variables (r.v.) x i , X2, x 3 , we have the following: p(xl\x2) = j™x p(xx \x2,x3 )p(x2 I x2 )dx2 ( A l ) for a Markov sequence of order one we have the following: p(xx | x2,x3) = p(xx | x2) (A2) Substituting xj = xk , x2 = xk-i, and x3 = zi:k-i into (A2) and then substituting this result in ( A l ) gives the (8). B. Derivation of the BRE Update Equation The derivation of (9) is outlined as follows: Derivation P(z\:k \xk)P(xk) P(xk \z\:k) = P(Zl;k) (A3> Equation (A3) is from Bayes' rule. Separating the pdf p(zi:k) into p(zk, zi.-k-i) and utilizing 143 factorized joint probability 2 results in the following: P("k \zi:k)= ' » » ^ = <A4) P(Zk>Z\:k-\) P(Zk I Zl.-t-l. xk )P(z\:k-\ I xk )P(xk ) p(tk \Z\:k-\)p(Z\:k-\) From Bayes' rule (A4) becomes as follows: P(xk \ z\k)= ^ Z k * Z l k ~ l ' X k ^ X k ^ Z l : k ~ l ^ Z l ± ~ 1 ^ X k ^ (A5) P(Zk | Z\:k-l )p(Zl.yt-l )P(Xk ) Canceling out terms in (A5) and using the assumption that we have independent measurements (white measurement noise) 3 results in the update equation (9). C. Derivation of the Recursive Weights Update Equation The derivation of (30) is outlined as follows: Derivation From (28), we have the following: ^P*(xQ:k\zi:k) i=l q(xl0:k\z\:k) 2p(a,b|c) = p(a|b,c)-p(b|c) and p(a,b) = p(a|b)-p(b) 3 P( Zk\Z\:k-\ >xk) = P(Zk\xk) Wk , ~i P (X0k \z\-k >r K , where w\= — — — (Af,\ N, . k j j ( A 6 ) Ywk 4(x0:k\zlk) i=l 144 Note: ~ i P*(xl-kK-k)„P(xO:k\zl:k) Wk= , t , t . x , t . t , (A7) If the importance density function is chosen to factorized such that: ^X0:k I Z\:k ) = 4(xk\ X0:k-\ > Z\:k Mx0:k-l I Z\:k > (A8) then one can augment old particles xlQ.k_i by x k ~a(xk \xok-i'z\k^ t 0 §et n e w particles xlQ.k. From (A5) we have the following: p(x0.k\z1.k)= P ( Z k 1 Zl:k~l' X°'k X ° : h * Z l k ~ l ) P ( Z l : k ~ 1 ) P ( X ° : k ^ (A9> p(Zk | Zl:k-\ )p(Zi;k-l )P(x0:k ) Canceling out terms in and using the assumption that we have independent measurements results in the following: p ( x Q . k \ z x . k ) = P ( Z k U Q : k ) p ( X k ' X ° : k - x U x : k - l ) (A10) P(Zk\Z\:k-\) Utilizing factorized joint probability on (A 10) gives the following: P(x0-k \z\k)= P ( Z k * X k ' X 0 : k ~ l ^ * k * X°:k~1' Z l : k ~ l ) p ( X°:k~l * Zl:k~1) ( A l 1) P(Zk I Zi.*-i ) 145 Using the fact that we have a first order Markov process results in the following: P(x0:k I Z l . * ) = P(zk \ xk )p(xk I xk-l) P( zk i ; P(xQ:k-\\Z\:k-\) (A12) or P(x0:k I Zl:k ) x(Zk\xk )P(xk I xk-\ )P(x0:k-\ I Zl;jfc-1) (A13) Substituting (A 13) and (A8) into (A7) gives the following: ~ i ~i w k = w k -P(zk \X1)P(X <*(xk\xO:k-Vzl:k-0 (AU) Equation (A14) becomes (30) i f it is assumed that q(xlk \ xlQ.k_i>z\.k) = q(x\ I xk-\yZlk) (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 of s(t) and n(t) d(t) = s(t) + n(t) (A15) 146 The joint probability density of d(t) and s(t) is defined as follows: pds(d,s) = p(d\s)ps(s) (A16) 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 of s we have the following [31]: p(d\s) = p(n\s) = pn(n) = pn(d-s) (A18) 147 E . PS-SEED F S M C D and H M M Filter Implementation A s outlined in Table 4.1, the first step of the P S - S E E D algorithm is to initialize a bank o f K F s equal in number (Ns) to the particles pre-specified. The P S - S E E D algorithm utilizes a F S M C D (within the J M G L S formulation) to sequentially update the set of K F s for each time increment as was outlined in Table 3.1. The two states of the F S M C D ylk ~ p{ylk | y(_i) l ? are event (p(yk)) and no-event (p(yk)). A Markov chain is uniquely defined by the initial distribution and the transitional probabilities at time k = 0. For the simulations presented in Section I V . A . l . c , the initial distribution for the no-event and event states was 1 2 set to p(y0 ) = 0.9 and p(yQ) = 0.1, respectively. The initial transitional probabilities were specified as p(y\ \y\_l) = 0.8, p(y\ \y\_x) = 0.8, p(y\ \y\_x) = b.2, and 2 1 p(yk I yk_x) = 0.2. In general terms, the probability o f a source wavelet being present ("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 of the initial F S M C D probabilities as long as a sufficient number of particles are implemented. A t each time increment k the F S M C D is updated according to the following equation: p(y\) p(y\). In addition to updating the pdf of the F S M C D at each time increment k, a Monte Carlo technique is utilized to obtain realizations of the F S M C D equal to the number of event/no-event particles (Ns). In this step, a random number generator is utilized to obtain Ns samples of the uniform distribution If [0,1]k (i = 1 to If If [0,1]k < p(y\ ) then an y\\y\-\ y\\y2k-x yl\y\-x y\\y\-i P(A-0 p(yli) (A19) 148 event has occurred otherwise no-event is present. In terms of 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^ +vk) as outlined in Step 1, Case 1 of Table 4.1. For the event condition, the measurement equation given as the addition of the background noise and source wavelet (e.g., = x l k + xlk sin(cok&. + (pk) + vk) as outlined in Step 1, Case 2 of Table 4.1. is defined by the overlapping source wavelet combination (e.g., (22) for the case of three overlapping source wavelets). Based upon the output of the bank of K F s , the PS-SEED algorithm utilizes a R B P F to obtain asymptotically optimal estimate of the modulating amplitude term of the A M S source wavelet defined by the state vector (Steps 6 and 7 of Table 4.1). From the output of the R B P F filter, a H M M filter is utilized to refine the initially specified phase component of the sinusoid within the user specified phase resolution window (Step 9 of Table 4.1). Once the PS-SEED algorithm obtains an optimal phase estimate at time index k, (18)-(20) of Table 3.1 are utilized to update the bank of K F s . The time index k is then incremented and the P S - S E E D algorithm iterates to Step 4 of 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 of the PPD algorithm is to initialize a bank of K F s equal in number (Ns) to that defined by (65). The PPD algorithm utilizes a F S M C D (within the J M G L S formulation) to sequentially update the set of K F s for each time increment as was outlined in Table 3.1. The two states of the F S M C D yj^ ~ p[yk \ y{_\) a r e event 1 2 (p(yk)) and no-event (p(yk))- A Markov chain is uniquely defined by the initial distribution and the transitional probabilities at time k = 0. For the simulations presented in 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, respectively. The initial transitional probabilities were specified as p(y\ \y{_l) = 0A, p(y\ \ yk_x) = OA, p(y\ \ y\_x) = 0.6, and 2 1 P(yk I yk-l) = 0-6- I n general terms, the probability of 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 PPD algorithm behaves robustly to the specification of the initial F S M C D probabilities as long as a sufficient number of particles are implemented. At 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 of the F S M C D equal to the number of event/no-event particles ( A ^ ) . In this step, a random number generator is utilized to obtain N$ samples of the uniform distribution If[0,1]k (i = 1 to Nfs). If If[0,1]\ < p(y\^then an event has occurred otherwise no-event is present. For the event condition the number of K F s sequentially updated is equal to the calculated number of combinations defined by (64). In terms of 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 measurement equation is equal to the 150 background noise (zk -xAk). For the event condition, the measurement equation is defined by the overlapping source wavelet combination (e.g., (63) for the case of three overlapping source wavelets). Based upon the output of the bank of K F s , the P P D algorithm utilizes a R B P F to obtain asymptotically optimal estimates of the modulating amplitude terms of the A M S source wavelets defined by the state vector (Steps 6 and 7 of Table 4.3). From the output of the R B P F filter, H M M filters are utilized to refine the initially specified phase components of the sinusoids within the user specified phase resolution window (Step 9 of 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 of the pdf of the A M S source wavelet phases was set to the uniform distribution, while the fixed-grid transitional pdf p(<plk \<plk-\) was set quite high (e.g., 0.996) and the remaining values 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 (Np) phases specified for each phase resolution window to reflect the possible phase shifts of ±3° in 1° increments. Once the P P D obtains optimal phase estimates at time index k, (18)-(20) of Table 3.1 are utilized to update the bank of KFs . The time index k is then incremented and the P P D algorithm iterates to Step 4 of Table 4.1. 151 G. Implementation of the Three State F S M C D within the PPD-WE Algorithm present at time index k (p(yk)), the source wavelet understudy plus measurement noise present at time index k(p(yk)) and the source wavelet understudy is overlapped with other time series data plus measurement noise at time index k(p(yk)). 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 in Section IVB.2.f, the initial distribution for the three state F S M C D was set to p(ylQ) = 0.5, p(y2 ) = 0.4, and p(y^) = 0.1, respectively. This initial distribution is based upon the fact that at the start of the time series under analysis only measurement noise or the source wavelet understudy wi l l be present (ability to window in on initial seismic wavelet recordings). The initial transitional probabilities were specified as 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 of (A20) add up to 1.0 according to the law of total probability. A t each time increment k the pdf of the three state F S M C D is updated according to the following equation: The three states of the P P D - W E F S M C D measurement noise y\\yl0 y\\yi y\\y\ [0.1429 0.0357 0.0909' yf\ylo y \ \ y \ y\\y\ = 0.5714 0.3214 0.0909 y\\y\ y\\yl y\\y\ L 0 - 2 8 5 7 ° - 6 4 2 9 ° - 8 1 8 2 (A20) 152 P(y\) ~y\ \yU y\ \yU y\ \yU p(y\-0 p(y2k) = yl \yU yl \yU yl \yU p(y\-\) (A21) p(y\) y\ yU yl \yU yl y\-\ p(y\-x) The calculated probabilities p(y\) o f (A21) are sorted from lowest to highest. The Inverse Transform Method is then utilized to obtain realizations of the three state F S M C D equal to the number of particles, Ns, specified by the investigator. In this step, a random number generator is utilized to obtain Ns samples of the uniform distribution If [0,1] k (i - 1 to Ns ). The realizations of the three state F S M C D are then calculated as follows: yk y\ 2 yk yl if Ui[QX]k<p{yk) if P(yl)< uUo.iJk < p{y\)+p{yl) otherwise (A22) The realizations outlined in (A22) are based upon the assumption 1 2 3 that p(yk)<p(yk)<p(yk). In terms of the J M G L S , the system equation (67 does not change but the measurement equation is updated based upon the estimated value of y'k. For the case yjj. =y\ the measurement equation is set equal to the background noise ( zk = x5k ). For the case ylk = v | the measurement equation is set as zk =xlk sin((okA + (p\ ^ ) + x 5 k (A23) For the case y'k = the measurement equation is set to (68 and states xlk and x3k are jittered (i.e., add on U[-1.5,+1.5]) to facilitate a diversity of particles and minimize sample impoverishment [35]. 153 H . PPD-WE Estimated Source and Residual Wavelets for SCPT 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 60 70 SO 90 "10O T i m e [ m s ] Figure H.2. The Averaged Estimated First Arriving Source Wavelet at 3.0m. 1 ; o.a; -1 "l • • • , r-> . . . . . . , I O 10 20 30 40 50 60 TO SO 90 10O Time [ms] Figure H.3. The estimated first arriving source wavelet superimposed upon residual wavelet at 3.0 m. 154 2) Estimated Source Wavelet at 4m 1:1 Time [ms] Figure H.4. Thirty five P P D - W E primary source wavelet estimations at depth of 4.0 m 1:: Time [ms] Figure H.5. The Averaged Estimated First Arriv ing Source Wavelet at 4.0m. o.s; Time [ms] 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 3) Estimated Source Wavelet at 5m 100 T i m e [ms] Figure H.8. The Averaged Estimated First Arriving Source Wavelet at 5.0m. 1 T i m e [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:: Time [ms] Figure H . l l . The Averaged Estimated First Arriving Source Wavelet at 6.0m. 1-. — — °-8: t • * • < • ; < < < : > i ^ ; j , , ! ! , , • , f O 25 50 75 100 Time [ms] Figure H.12. The estimated first arriving source wavelet superimposed upon residual wavelet at 6.0 m. 157 5) Estimated Source Wavelet at 7m 50 T i m e [ m s ] Figure H.13. Thirty five P P D - W E primary source wavelet estimations at depth of 7.0 m. T i m e [ m s ] Figure H.14. The Averaged Estimated First Arriving Source Wavelet at 7.0m. T i m e [ m s ] Figure H.15. The estimated first arriving source wavelet superimposed upon residual wavelet at 7.0 m. 158 6) Estimated Source Wavelet at 8m 100 Time [ms] Figure H.16. Thirty five P P D - W E primary source wavelet estimations at depth of 8.0 m. -in o--1 1 1 ' ; 1 i •• •• ' •• 1 1 • •• 1 ' 1 1 1 • • • : 1 O 25 50 75 IOO Time [ms] 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 7) Estimated Source Wavelet at 9m '1 : : . . . 1 . . . . , i , , , , , , , , , , , ( O 25 50 75 100 Time [ms] 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 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 75 100 Time [ms] Figure H.25. Thirty five P P D - W E primary source wavelet estimations at depth of 11.0m. 1 ;; — .8 Time [ms] Figure H.26. The Averaged Estimated First Arr iv ing Source Wavelet at 11.0m. yesian Recursive Estimation 1. _ _ 0.8--1 \ , , , . , , ^ , , , , , , , , , . . , 1 O 25 50 75 100 Time [ms] Figure H.27. The estimated first arriving source wavelet superimposed upon residual wavelet at 11.0 m. 162 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 5 0 T i m e [ m s ] Figure H.31. Thirty five P P D - W E primary source wavelet estimations at depth of 13.0m. 5 0 T i m e [ m s ] Figure H.32. The Averaged Estimated First Arriving Source Wavelet at 13.0m. 5 0 T i m e [ m s ] 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 T i m e [ m s ] Figure H.34. Thirty five P P D - W E primary source wavelet estimations at depth of 14.0m O 25 50 75 100 T i m e [ m s ] Figure H.35. The Averaged Estimated First Arriving Source Wavelet at 14.0m. 1:1 T i m e [ m s ] Figure H.36. The estimated first arriving source wavelet superimposed upon residual wavelet at 14.0 m. 165 T i m e [ m s ] Figure H.38. The Averaged Estimated First Arriving Source Wavelet at 15.0m. O 2 5 5 0 7 5 1 0 0 T i m e [ m s ] 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 75 100 T i m e [ m s ] Figure H.40. Thirty five P P D - W E primary source wavelet estimations at depth of 16.0m. o3 O 25 50 75 100 T i m e [ m s ] Figure H.41. The Averaged Estimated First Arriving Source Wavelet at 16.0m. 1 o o -• -i •• • 1 1 1 1 > . . . > 1 . , , : ! 1 > 1 ! . 1 O 25 50 75 100 Tim© [ms] Figure H.42. The estimated first arriving source wavelet superimposed upon residual wavelet at 16.0 m. 167 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
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
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 |
AggregatedSourceRepository | 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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0052619/manifest