marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

initial import of classes for automatic picking purposes [just imported by me; module has originally been written by Ludger Küperkoch]

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2014-11-14 07:40:00 +01:00
parent c40aec192c
commit fbce83293d
1 changed files with 342 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

342
pylot/core/pick/CharFuns.py Normal file
View File

@ -0,0 +1,342 @@
# -*- coding: utf-8 -*-
"""
Created Oct/Nov 2014
Implementation of the Characteristic Functions (CF) published and described in:
Kueperkoch, L., Meier, T., Lee, J., Friederich, W., & EGELADOS Working Group, 2010:
Automated determination of P-phase arrival times at regional and local distances
using higher order statistics, Geophys. J. Int., 181, 1159-1170
Kueperkoch, L., Meier, T., Bruestle, A., Lee, J., Friederich, W., & EGELADOS
Working Group, 2012: Automated determination of S-phase arrival times using
autoregressive prediction: application ot local and regional distances, Geophys. J. Int.,
188, 687-702.
:author: MAGS2 EP3 working group
"""
import numpy as np
from obspy.core import Stream
class CharacteristicFunction(object):
'''
SuperClass for different types of characteristic functions.
'''
def __init__(self, data, cut, t2, order, t1=None, fnoise=0.001):
'''
Initialize data type object with information from the original
Seismogram.
:param: data
:type: `~obspy.core.stream.Stream`
:param: cut
:type: tuple
:param: t2
:type: float
:param: order
:type: int
:param: t1
:type: float (optional, only for AR)
:param: fnoise
:type: float (optional, only for AR)
'''
assert isinstance(data, Stream), "%s is not a Stream object" % str(data)
self.orig_data = data[0]
self.dt = self.orig_data.stats.delta
self.setCut(cut)
self.setTime1(t1)
self.setTime2(t2)
self.setOrder(order)
self.setFnoise(fnoise)
self.calcCF(self.getDataArray())
self.arpara = np.array([])
self.xpred = np.array([])
def __str__(self):
return '''\n\t{name} object:\n
Cut:\t\t{cut}\n
t1:\t{t1}\n
t2:\t{t2}\n
Order:\t\t{order}\n
Fnoise:\t{fnoise}\n
'''.format(name=type(self).__name__,
cut=self.getCut(),
t1=self.getTime1(),
t2=self.getTime2(),
order=self.getOrder(),
fnoise=self.getFnoise())
def getCut(self):
return self.cut
def setCut(self, cut):
self.cut = cut
def getTime1(self):
return self.t1
def setTime1(self, t1):
self.t1 = t1
def getTime2(self):
return self.t2
def setTime2(self, t2):
self.t2 = t2
def getOrder(self):
return self.order
def setOrder(self, order):
self.order = order
def getIncrement(self):
return self.dt
def getFnoise(self):
return self.fnoise
def setFnoise(self, fnoise):
self.fnoise = fnoise
def getCF(self):
return self.cf
def getDataArray(self, cut=None):
'''
If cut times are given, time series is cut from cut[0] (start time)
till cut[1] (stop time) in order to calculate CF for certain part
only where you expect the signal!
input: cut (tuple) ()
cutting window
'''
if cut is not None:
if self.cut[0] == 0:
start = 0
else:
start = self.cut[0] / self.dt
stop = self.cut[1] / self.dt
data = self.orig_data.data[start:stop]
return data
return self.orig_data.data
def calcCF(self, data=None):
self.cf = data
def arDet(self, data, order, rind, ldet):
pass
def arPred(self, data, arpara, rind, lpred):
pass
class AICcf(CharacteristicFunction):
'''
Function to calculate the Akaike Information Criterion (AIC) after
Maeda (1985).
:param: data, time series (whether seismogram or CF)
:type: tuple
Output: AIC function
'''
def calcCF(self, data):
print 'Calculating AIC ...'
xnp = self.getDataArray()
datlen = len(xnp)
k = np.arange(1, datlen)
cf = np.zeros(datlen)
cumsumcf = np.cumsum(np.power(xnp, 2))
i = np.where(cumsumcf == 0)
cumsumcf[i] = np.finfo(np.float64).eps
cf[k] = ((k - 1) * np.log(cumsumcf[k] / k) + (datlen - k + 1) *
np.log((cumsumcf[datlen - 1] -
cumsumcf[k - 1]) / (datlen - k + 1)))
cf[0] = cf[1]
inf = np.isinf(cf)
ff = np.where(inf == 'True')
if len(ff) >= 1:
cf[ff] = 0
self.cf = cf - np.mean(cf)
class HOScf(CharacteristicFunction):
'''
Function to calculate skewness (statistics of order 3) or kurtosis
(statistics of order 4), using one long moving window, as published
in Kueperkoch et al. (2010).
'''
def calcCF(self, data):
xnp = self.getDataArray(self.getCut())
if self.getOrder() == 3: # this is skewness
print 'Calculating skewness ...'
y = np.power(xnp, 3)
y1 = np.power(xnp, 2)
elif self.getOrder() == 4: # this is kurtosis
print 'Calculating kurtosis ...'
y = np.power(xnp, 4)
y1 = np.power(xnp, 2)
# Initialisation
# t2: long term moving window
ilta = round(self.getTime2() / self.getIncrement())
lta = y[0]
lta1 = y1[0]
# moving windows
LTA = np.zeros(len(xnp))
for j in range(3, len(xnp)):
if j <= ilta:
lta = (y[j] + lta * (j - 1)) / j
lta1 = (y1[j] + lta1 * (j - 1)) / j
else:
lta = (y[j] - y[j - ilta]) / ilta + lta
lta1 = (y1[j] - y1[j - ilta]) / ilta + lta1
# define LTA
if self.getOrder() == 3:
LTA[j] = lta / np.power(lta1, 1.5)
elif self.getOrder() == 4:
LTA[j] = lta / np.power(lta1, 2)
LTA[0:3] = 0
self.cf = LTA
class ARZcf(CharacteristicFunction):
def calcCF(self, data):
print 'Calculating AR-prediction error from single trace ...'
xnp = self.getDataArray(self.getCut())
# some parameters needed
# add noise to time series
xnoise = xnp + np.random.normal(0.0, 1.0, len(xnp)) * self.getFnoise() * max(abs(xnp))
tend = len(xnp)
# Time1: length of AR-determination window [sec]
# Time2: length of AR-prediction window [sec]
ldet = int(round(self.getTime1() / self.getIncrement())) # length of AR-determination window [samples]
lpred = int(np.ceil(self.getTime2() / self.getIncrement())) # length of AR-prediction window [samples]
cf = []
step = ldet + self.getOrder() - 1
for i in range(ldet + self.getOrder() - 1, tend - lpred + 1):
if i == step:
'''
In order to speed up the algorithm AR parameters are kept for time
intervals of length lpred
'''
# determination of AR coefficients
self.arDet(xnoise, self.getOrder(), i - ldet, i)
step = step + lpred
# AR prediction of waveform using calculated AR coefficients
self.arPred(xnp, self.arpara, i + 1, lpred)
# prediction error = CF
err = np.sqrt(np.sum(np.power(self.xpred[i:i + lpred] - xnp[i:i + lpred], 2)) / lpred)
cf.append(err)
# convert list to numpy array
cf = np.asarray(cf)
self.cf = cf
def arDet(self, data, order, rind, ldet):
'''
Function to calculate AR parameters arpara after Thomas Meier (CAU), published
in Kueperkoch et al. (2012). This function solves SLE using the Moore-
Penrose inverse, i.e. the least-squares approach.
:param: data, time series to calculate AR parameters from
:type: array
:param: order, order of AR process
:type: int
:param: rind, first running summation index
:type: int
:param: ldet, length of AR-determination window (=end of summation index)
:type: int
Output: AR parameters arpara
'''
# recursive calculation of data vector (right part of eq. 6.5 in Kueperkoch et al. (2012)
rhs = np.zeros(self.getOrder())
for k in range(0, self.getOrder()):
for i in range(rind, ldet):
rhs[k] = rhs[k] + data[i] * data[i - k]
# recursive calculation of data array (second sum at left part of eq. 6.5 in Kueperkoch et al. 2012)
A = np.array([[0, 0], [0, 0]])
for k in range(1, self.getOrder() + 1):
for j in range(1, k + 1):
for i in range(rind, ldet):
ki = k - 1
ji = j - 1
A[ki, ji] = A[ki, ji] + data[i - ji] * data[i - ki]
A[ji, ki] = A[ki, ji]
# apply Moore-Penrose inverse for SVD yielding the AR-parameters
self.arpara = np.dot(np.linalg.pinv(A), rhs)
def arPred(self, data, arpara, rind, lpred):
'''
Function to predict waveform, assuming an autoregressive process of order
p (=size(arpara)), with AR parameters arpara calculated in arDet. After
Thomas Meier (CAU), published in Kueperkoch et al. (2012).
:param: data, time series to be predicted
:type: array
:param: arpara, AR parameters
:type: float
:param: rind, first running summation index
:type: int
:param: lpred, length of prediction window (=end of summation index)
:type: int
Output: predicted waveform z
'''
# be sure of the summation indeces
if rind < len(arpara) + 1:
rind = len(arpara) + 1
if rind > len(data) - lpred + 1:
rind = len(data) - lpred + 1
if lpred < 1:
lpred = 1
if lpred > len(data) - 1:
lpred = len(data) - 1
z = np.append(data[0:rind], np.zeros(lpred))
for i in range(rind, rind + lpred):
for j in range(1, len(arpara) + 1):
ji = j - 1
z[i] = z[i] + arpara[ji] * z[i - ji]
self.xpred = z
class ARHcf(CharacteristicFunction):
pass
class AR3Ccf(CharacteristicFunction):
pass