forked from uafgeotools/array_processing
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathdetection.py
More file actions
91 lines (65 loc) · 2.63 KB
/
Copy pathdetection.py
File metadata and controls
91 lines (65 loc) · 2.63 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
import numpy as np
from obspy.core import Stream
from .generic import phaseAlignIdx, phaseAlignData
def fstatbland(dtmp, fs, tau):
r"""
Calculates the F-statistic and SNR based on Blandford's method (see Notes).
Args:
dtmp: ``(m, n)`` array; time series with ``m`` samples from ``n``
traces as columns
fs (int or float): Sample rate [Hz]
tau: ``(n(n-1)//2)`` array; time delays of relative signal arrivals
(TDOA) for all unique sensor pairings
Returns:
tuple: Tuple containing:
- **fstat** – F-statistic
- **snr** – SNR
References:
Blandford, R. R., 1974. An automatic event detector at the Tonto
Forest Seismic Observatory. Geophysics, vol. 39, no. 5,
p. 633–643. `https://library.seg.org/doi/abs/10.1190/1.1440453
<https://library.seg.org/doi/abs/10.1190/1.1440453>`__
"""
m, n = dtmp.shape
wgt = np.ones(n)
#individual trace offsets from arrival model shifts. Zeros here
m_offset = [0 for i in range(n)]
# calculate beam delays
beam_delays = phaseAlignIdx(tau, fs, wgt, 0)
# apply shifts, resulting in a zero-padded array
beam = phaseAlignData(dtmp, beam_delays, wgt, 0, m, m_offset)
fnum = np.sum(np.sum(beam, axis=1)**2)
term1 = np.sum(beam, axis=1)/n
term1_0 = term1
for i in range(1, n):
term1 = np.vstack((term1, term1_0))
fden = np.sum(np.sum((beam.T - term1)**2))
fstat = (n-1) * fnum/(n * fden)
#calculate snr based on fstat
snr = np.sqrt((fstat-1)/n)
return fstat, snr
def calculate_semblance(data_in):
r"""
Calculates the semblance, a measure of multi-channel coherence, following
the definition of Neidell & Taner (1971). Assumes data are already
time-shifted to construct the beam.
Args:
data_in: Time-shifted ObsPy Stream or time-shifted NumPy array
Returns:
Multi-channel coherence defined on :math:`[0, 1]`
"""
if isinstance(data_in, Stream):
# check that all traces have the same length
if len(set([len(tr) for tr in data_in])) != 1:
raise ValueError('Traces in stream must have same length!')
n = len(data_in)
beam = np.sum([tr.data for tr in data_in], axis=0) / n
beampower = n * np.sum(beam**2)
avg_power = np.sum(np.sum([tr.data**2 for tr in data_in], axis=0))
elif isinstance(data_in, np.ndarray):
n = data_in.shape[0]
beam = np.sum(data_in, axis=0) / n
beampower = n * np.sum(beam**2)
avg_power = np.sum(np.sum([data_in**2], axis=0))
semblance = beampower / avg_power
return semblance