Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Application of Bayesian recursive estimation for seismic signal processing Baziw, Erick 2007

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
ubc_2007-266762.pdf [ 7.03MB ]
Metadata
JSON: 1.0052619.json
JSON-LD: 1.0052619+ld.json
RDF/XML (Pretty): 1.0052619.xml
RDF/JSON: 1.0052619+rdf.json
Turtle: 1.0052619+rdf-turtle.txt
N-Triples: 1.0052619+rdf-ntriples.txt
Original Record: 1.0052619 +original-record.json
Full Text
1.0052619.txt
Citation
1.0052619.ris

Full Text

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

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 14 4
United States 10 5
Poland 7 0
Japan 4 0
Indonesia 4 15
India 3 0
Canada 3 0
Russia 1 0
Ukraine 1 0
United Kingdom 1 0
Republic of Korea 1 1
Philippines 1 0
Norway 1 0
City Views Downloads
Unknown 14 4
Shenzhen 7 4
Ashburn 6 0
Beijing 4 0
Tokyo 4 0
Vancouver 3 0
Mumbai 2 0
Chengdu 2 0
Sunnyvale 2 0
Wuhan 1 0
Makati City 1 0
Newcastle upon Tyne 1 0
Albuquerque 1 1

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052619/manifest

Comment

Related Items