Bayesian Deconvolution of Seismic Array Data

for Ripple-Fired Explosions Eric A. Suess Advisor: Professor Robert Shumway Outline:

1. The Problem

2. Model

3. Bayesian Statistics

4. Bayesian Deconvolution of Seismic Array Data 5. Results

• Simulated Data

• Real Data

6. Conclusions and Future Work Treaties:

1. 1963 Limited Nuclear Test Ban Treaty (LTBT)

2. 1968 Non-Proliferation of Nuclear Weapons Treaty (NPT)

3. 1974 Threshold Test Ban Treaty (TTBT)

4. 1976 Peaceful Nuclear Explosions Treaty (PNET)

5. 1996 Comprehensive Test Ban Treaty (CTBT)

Background: Much of the focus in the past has been on distinguishing possible \underline{nuclear explosions} from \underline{earthquakes}. Currently, since the testing treaties have put limitations on the permissible sizes of the nuclear explosions, other smaller seismic events such as \underline{industrial mining explosions} have become of interest in the discrimination problem. The work presented here is related to distinguishing low-level nuclear explosions from ripple-fired mining explosions that are on the same seismic level.

The Problem:

• Monitoring seismic events at Regional Distances for low-level nuclear tests.

• Other seismic sources need to be ruled out.

– Ripple-Fired Mining Explosions Ripple-Fired Explosions: This is a mining technique in which explosions of single devices (or groups of devices) are detonated in succession.

Monitoring: Arrays of receivers are put in place at Regional Distances and seismic data is continually collected. Seismic disturbances that are above the baseline noise of the area are investigated.

Model:

y_{k}(t) = s_{k}(t) + ∑_{j=1}^{m} a_{j}s_{k}(t − j) + ε_{k}(t)

Amplitudes are distributed according to a random

Bernoulli-Gaussian model. (ref. Cheng, Chen and Li 1996)

p(a_{j}|η) ∼ (1 − η)I(a_{j} = 0) + ηT N(μ_{α}, σ_{α}^{2})I(a_{j} > 0)

Signal and path effects follow an AR(3) model. (ref.

Dargahi-Noubary 1995, Tjøstheim 1975)

s_{k}(t) + ϕ_{1}s_{k}(t − 1) + ... + ϕ_{p}s_{k}(t − p) = e_{k}(t)

e_{k}(t) i.i.d N(0, σ^{2}) and define the precision τ = \frac{1}{σ^{2}} . ε_{k}(t) i.i.d. N(0, cσ^{2}), c = 1/SNR, c ≥ 0 is fixed.

Truncated Normal:

p(a_{j}) = c^{−1}(μ_{α}, σ_{α}^{2})\frac{1}{\sqrt{2π}σ_{α}} exp \{−\frac{(a_{j} − μ_{α})^{2}}{σ_{α}^{2}} \} I(a_{j} > 0)

where

c^{−1}(μ_{α}, σ_{α}^{2}) = ∫_{0}^{∞} \frac{1}{\sqrt{2π}σ_{α}} exp \{−\frac{(x − μ_{α})^{2}}{σ_{α}^{2}} \}

Bayesian Statistics: Model:

p(Y|Θ)

Prior distribution:

p(Θ)

Joint distribution:

p(Y, Θ) = p(Y|Θ)p(Θ)

Posterior distribution:

p(Θ|Y) = \frac{p(Y|Θ)p(Θ)}{∫ p(Y|Θ)p(Θ)dΘ}
∝ p(Y|Θ)p(Θ)

Gibbs Sampler: (ref. Gelfand and Smith 1990)

Θ = (θ_{1}, θ_{2}) p(Θ|Y) = p((θ_{1}, θ_{2})|Y) Given θ_{1}^{(0)} and θ_{2}^{(0)} for h = 1, ..., Reps

1. Sample θ_{1}^{(h)} from p(θ_{1}|Y, θ_{2}^{(h−1)})

2. Sample θ_{2}^{(h)} from p(θ_{2}|Y, θ_{1}^{(h)})

3. Set h = h + 1 and go to 1.

(θ_{1}^{(BurnIn+1)}, θ_{2}^{(BurnIn+1)}), ..., (θ_{1}^{(Reps)}, θ_{2}^{(Reps)})

Θ^{(1)}, ..., Θ^{(Reps)} are realizations of a stationary Markov Chain, with transition probability from Θ^{(h−1)} to Θ^{(h)},

T (Θ^{(h−1)}, Θ^{(h)}) = p(θ_{1}|Y, θ_{2}^{(h−1)})p(θ_{2}|Y, θ_{1}^{(h)})

By Ergodic theory, we can calculate estimates of say θ_{1} by

\frac{1}{Reps} ∑_{h} θ_{1}^{(i)} ⟶^{a.s.} E[θ_{1}|Y],

Reps → ∞.

Priors: ; η ∼ BET A(β_{1}, β_{2}) ; ϕ ∼ N_{p}(ϕ_{0}, Σ_{0}) ; τ ∼ GAMMA(γ_{1}, γ_{2})

The Model:

y_{k}(t) = s_{k}(t) + ∑_{j=1}^{m} a_{j}s_{k}(t − j) + ε_{k}(t)

The parameter set: Θ = {η, a, ϕ, τ, S} Hyperparamters: β_{1}, β_{2}, μ_{α}, σ_{α}^{2}, ϕ_{0}, Σ_{0}, γ_{1}, γ_{2} Fixed parameters: c, m Likelihood: Y = [y_{1}, ..., y_{q}]′ y_{k} = [y_{k}(1), ..., y_{k}(n)]′

p(Y|Θ) = ∏_{k=1}^{q} p(y_{k}|Θ)
= ∏_{k=1}^{q}(2π)^{−n/2} (\frac{τ}{c} )^{n/2} exp \{−\frac{τ}{2c} ∑_{t} ε_{k}^{2}(t)\}

Overall Prior:

p(Θ) = p(η, a, ϕ, τ, S)

= p(η)p(a|η)p(ϕ)p(τ)p(S|ϕ, τ)

Joint Density:

p(Y, Θ) = p(Y|Θ)p(Θ)

Joint Posterior:

p(Θ|Y) ∝ p(Y|Θ)p(Θ)

Conditional Marginal Posterior Distributions:

p(η|Y, rest) ∼ beta(β_{1}^{∗}, β_{2}^{∗})

For fixed j = 1, ..., m

p(a_{j}|Y, rest) ∼ (1 − η_{j})I(a_{j} = 0) + η_{j}T N(μ_{a_{j}} , σ_{a_{j}}^{2} )I(a_{j} > 0)

p(ϕ|Y, rest) ∼ N_{p}(ϕ_{∗}, Σ_{∗})
p(τ|Y, rest) ∼ gamma(γ_{1}^{∗}, γ_{2}^{∗})

For fixed i = 1, ..., n and k = 1, ..., q

p(s_{k}(i)|Y, rest) ∼ N(μ_{s_{k}(i)}, σ_{s_{k}(i)}^{2})

\begin{matrix}η|Y, rest ∼ beta(β_{1}^{∗}, β_{2}^{∗}) \\ β_{1}^{∗} = m − n_{a} + β_{1}, \\ β_{2}^{∗} = n_{a} + β_{2} \end{matrix}

a_{j}|Y, rest ∼ Bernoulli-Gaussian

μ_{a_{j}} = σ_{a_{j}}^{2} [\frac{μ_{α}}{σ_{α}^{2}} + (\frac{τ}{c} ) ∑_{k} ∑_{t} ε_{k}^{∗}(t[−j])s_{k}(t − j)]

σ_{a_{j}}^{−2} = \frac{1}{σ_{α}} + (\frac{τ}{c} ) ∑_{k} ∑_{t} s_{k}^{2}(t − j)

η_{j} = \frac{η}{η + (1 − η) [\frac{c(μ_{a_{j}} ,σ_{a_{j}}^{2} )}{c(μ_{α},σ_{α}^{2})}] (\frac{σ_{a_{j}}}{σ_{α}} ) exp \{\frac{1}{2} [\frac{μ_{a_{j}}^{2}}{σ_{a_{j}}^{2}} − \frac{μ_{α}^{2}}{σ_{α}^{2}} ]\}}

ϕ|Y, rest ∼ N_{p}(ϕ_{∗}, Σ_{∗})

ϕ_{∗} = Σ_{∗} [−τ ∑_{k} ∑_{t} s_{k}(t)s̃_{k}(t) + Σ_{0}^{−1}ϕ_{0}]

Σ_{∗}^{−1} = τ ∑_{k} ∑_{t} s̃_{k}(t)s̃′_{k}(t) + Σ_{0}^{−1}

where

s̃_{k}(t) = [s_{k}(t − 1), ..., s_{k}(t − p)]′

\begin{matrix}τ|Y, rest ∼ gamma(γ_{1}^{∗}, γ_{2}^{∗}) \\ γ_{1}^{∗} = q(2n − l − p)/2 + γ_{1} \\ γ_{2}^{∗} = \frac{1}{2c} ∑_{k} ∑_{t} ε_{k}^{2}(t) + \frac{1}{2} ∑_{k} ∑_{t} e_{k}^{2}(t) \end{matrix}

s_{k}(t)|Y, rest ∼ Normal ; μ_{s_{k}(i)} = σ_{s_{k}(i)}^{2} [(\frac{τ}{c} ) ∑_{t} a_{t−i}ε′(t[−i]) − τ ∑_{t} ϕ_{t−i}e′_{k}(t[−i])] ; σ_{s_{k}(i)}^{−2} = (\frac{τ}{c} ) ∑_{t} a_{t−i}^{2} + τ ∑_{t} ϕ_{t−i}^{2}

Steps To Perform The Gibbs Sampler: Given the initial values:

{η^{(0)}, a^{(0)}, ϕ^{(0)}, τ^{(0)}, S^{(0)}}

for h = 1 to Reps:

1. Sample η^{(h)} from a beta.

2. Sample a_{j}^{(h)}, j = 1, ..., m, from a Bernoulli-Gaussian.

3. Sample ϕ^{(h)} from a p-variate normal.

4. Sample τ^{(h)} from a gamma.

5. Sample s_{k}(i)^{(h)}, i = 1, ..., n, k = 1, ..., q, from a normal.

6. Set h = h + 1 and go to Step 1.

Conclusions and Further Work:

1. one

2. two References: Cheng, Q., Chen, R., and Li,T.-H. (1996) Simultaneous Wavelet Estimation and Deconvolution of Reflection Seismic Signals IEEE Transactions On Geoscience and Remote Sensing 34, 377-384. Dargahi-Noubary. (1992) Stochastic Modeling and Identification of Seismic Records Based on Established Deterministic Formulations. Journal of Time Series Analysis 16, 201-220. Tjøstheim, D. (1975) Autoregressive Representation of Seismic P-wave Signals with an Application to the Problem of Short-Period Discriminants Geophys. J. R. astr. Soc. 43, 269-291.

Gelfand, A.E. and Smith, A.F.M. (1990) Sampling-Based Approaches to Calculating Marginal Densities Journal of the American Statistical Association 85, 398-409. Pearson, D.C., Strump, B.W., and Anderson, D.P. (1996) Physical Constraints on Mining Explosions http://www.geology.smu.edu/ dpa-www/papers/srl/SRL.html