From 728cf15f19ac9cc70729a4dc0bd441127cd1c602 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 24 Jun 2015 09:29:59 +0200 Subject: [PATCH 1/6] Demean: Mean to be subtracted from time series is now calculated from noise part only. --- pylot/core/pick/utils.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 58fd99a2..971ab1fb 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -47,14 +47,11 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): x = X[0].data t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, X[0].stats.delta) - # get latest possible pick - # get noise window inoise = getnoisewin(t, Pick1, TSNR[0], TSNR[1]) # get signal window isignal = getsignalwin(t, Pick1, TSNR[2]) # remove mean - meanwin = np.hstack((inoise, isignal)) - x = x - np.mean(x[meanwin]) + x = x - np.mean(x[inoise]) # calculate noise level nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac # get time where signal exceeds nlevel @@ -337,7 +334,7 @@ def getSNR(X, TSNR, t1): return # demean over entire snr window - x -= x[inoise[0]:isignal[-1]].mean() + x = x - np.mean(x[np.hstack([inoise, isignal])]) # calculate ratios noiselevel = np.sqrt(np.mean(np.square(x[inoise]))) From 7ede5a7859f1f96a7b3d5b8888e5ae24e19a194d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 24 Jun 2015 09:31:01 +0200 Subject: [PATCH 2/6] Modified parameters. --- autoPyLoT_local.in | 1 + 1 file changed, 1 insertion(+) diff --git a/autoPyLoT_local.in b/autoPyLoT_local.in index dd6d9815..68ded8a6 100644 --- a/autoPyLoT_local.in +++ b/autoPyLoT_local.in @@ -1,6 +1,7 @@ %This is a parameter input file for autoPyLoT. %All main and special settings regarding data handling %and picking are to be set here! +%Parameters are optimized for local data sets! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #main settings# From 61c745cd79a1ed35b6f5c62baf7934010b7b866f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 24 Jun 2015 09:31:18 +0200 Subject: [PATCH 3/6] Modified parameters. --- autoPyLoT_regional.in | 33 ++++++++++++++++----------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/autoPyLoT_regional.in b/autoPyLoT_regional.in index 3bc1bc9b..7d66c537 100644 --- a/autoPyLoT_regional.in +++ b/autoPyLoT_regional.in @@ -1,13 +1,14 @@ %This is a parameter input file for autoPyLoT. %All main and special settings regarding data handling %and picking are to be set here! +%Parameters are optimized for regional data sets! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #main settings# /DATA/Egelados #rootpath# %project path EVENT_DATA/LOCAL #datapath# %data path -2006.02_Nisyros #database# %name of data base -e0032.033.06 #eventID# %event ID for single event processing +2006.01_Nisyros #database# %name of data base +e1412.008.06 #eventID# %event ID for single event processing PILOT #datastructure# %choose data structure 2 #iplot# %flag for plotting: 0 none, 1, partly, >1 everything AUTOPHASES_AIC_HOS4_ARH #phasefile# %name of autoPILOT output phase file @@ -30,8 +31,8 @@ HYPOSAT #locrt# %location routine used ("HYPO #common settings picker# 20 #pstart# %start time [s] for calculating CF for P-picking 160 #pstop# %end time [s] for calculating CF for P-picking -1.0 #sstart# %start time [s] after or before(-) P-onset for calculating CF for S-picking -50 #sstop# %end time [s] after P-onset for calculating CF for S-picking +3.0 #sstart# %start time [s] after or before(-) P-onset for calculating CF for S-picking +100 #sstop# %end time [s] after P-onset for calculating CF for S-picking 3 10 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] 3 12 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] 3 8 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] @@ -48,13 +49,12 @@ HOS #algoP# %choose algorithm for P-onset 0.6 #tdet2z# %for AR-picker, length of AR determination window [s] for Z-component, 2nd pick 0.2 #tpred2z# %for AR-picker, length of AR prediction window [s] for Z-component, 2nd pick 0.001 #addnoise# %add noise to seismogram for stable AR prediction -4 0.1 1.0 0.5 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s] -4 #pickwinP# %for initial AIC pick, length of P-pick window [s] +4 0.2 2.0 1.5 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s] +4 #pickwinP# %for initial AIC and refined pick, length of P-pick window [s] 8 #Precalcwin# %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick) -0 #peps4aic# %for HOS/AR, artificial uplift of samples of AIC-function (P) -0.2 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s] -0.1 #tsmoothP# %for HOS/AR, take average of samples for smoothing CF [s] -0.001 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P) +3.0 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s] +0.3 #tsmoothP# %for HOS/AR, take average of samples for smoothing CF [s] +0.3 #ausP# %for HOS/AR, artificial uplift of samples (aus) of CF (P) 1.3 #nfacP# %for HOS/AR, noise factor for noise level determination (P) #H-components# ARH #algoS# %choose algorithm for S-onset determination (ARH or AR3) @@ -63,18 +63,17 @@ ARH #algoS# %choose algorithm for S-onset 0.6 #tdet2h# %for HOS/AR, length of AR-determinaton window [s], H-components, 2nd pick 0.3 #tpred2h# %for HOS/AR, length of AR-prediction window [s], H-components, 2nd pick 4 #Sarorder# %for AR-picker, order of AR process of H-components -20 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H) -5 #pickwinS# %for initial AIC pick, length of S-pick window [s] -6 0.2 2.0 1.5 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s] -0.05 #aictsmoothS# %for AIC-picker, take average of samples for smoothing of AIC-function [s] -0.02 #tsmoothS# %for AR-picker, take average of samples for smoothing CF [s] (S) -0.2 #pepsS# %for AR-picker, artificial uplift of samples of CF (S) +10 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H) +6 #pickwinS# %for initial AIC and refined pick, length of S-pick window [s] +5 0.2 3.0 3.0 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s] +3.0 #aictsmoothS# %for AIC-picker, take average of samples for smoothing of AIC-function [s] +1.0 #tsmoothS# %for AR-picker, take average of samples for smoothing CF [s] (S) 0.4 #ausS# %for HOS/AR, artificial uplift of samples (aus) of CF (S) 1.5 #nfacS# %for AR-picker, noise factor for noise level determination (S) %first-motion picker% 1 #minfmweight# %minimum required p weight for first-motion determination 2 #minFMSNR# %miniumum required SNR for first-motion determination -5.0 #fmpickwin# %pick window around P onset for calculating zero crossings +6.0 #fmpickwin# %pick window around P onset for calculating zero crossings %quality assessment% #inital AIC onset# 0.04 0.08 0.16 0.32 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P From b44fbe0b02060493b422b5d0708bc17281f03dfc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 24 Jun 2015 10:23:24 +0200 Subject: [PATCH 4/6] Some Cosmetics. --- autoPyLoT.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index 167fa8ae..c14514ee 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -13,6 +13,7 @@ from pylot.core.read import Data, AutoPickParameter from pylot.core.pick.run_autopicking import run_autopicking from pylot.core.util.structure import DATASTRUCTURE from pylot.core.pick.utils import wadaticheck +import pdb __version__ = _getVersionString() @@ -30,6 +31,18 @@ def autoPyLoT(inputfile): .. rubric:: Example ''' + print '************************************' + print '*********autoPyLoT starting*********' + print 'The Python picking and Location Tool' + print ' Version ', _getVersionString(), '2015' + print '**Authors:' + print '**S. Wehling-Benatelli' + print '** Ruhr-University Bochum' + print '**L. Kueperkoch' + print '** BESTEC GmbH' + print '**K. Olbert' + print '** Christian-Albrechts University Kiel' + print '************************************' # reading parameter file @@ -115,7 +128,6 @@ def autoPyLoT(inputfile): station = wfdat[0].stats.station allonsets = {station: picks} for i in range(len(wfdat)): - #for i in range(0,5): stationID = wfdat[i].stats.station #check if station has already been processed if stationID not in procstats: @@ -139,6 +151,10 @@ def autoPyLoT(inputfile): print '-------Finished event %s!-------' % parameter.getParam('eventID') print '------------------------------------------' + print '************************************' + print '*********autoPyLoT terminates*******' + print 'The Python picking and Location Tool' + print '************************************' if __name__ == "__main__": # parse arguments From 68bbea98545c66ad706eba350bd7180d3337ff5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 24 Jun 2015 14:15:54 +0200 Subject: [PATCH 5/6] Implemented new function for quality control: checksignallength, checks signal length in order to detect spuriously picked noise peaks. --- pylot/core/pick/utils.py | 86 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 971ab1fb..c869ef69 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -9,6 +9,7 @@ """ import numpy as np +import scipy as sc import matplotlib.pyplot as plt from obspy.core import Stream, UTCDateTime import warnings @@ -511,3 +512,88 @@ def wadaticheck(pickdic, dttolerance, iplot): plt.close(iplot) return checkedonsets + + +def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): + ''' + Function to detect spuriously picked noise peaks. + Uses envelope to determine, how many samples [per cent] after + P onset are below certain threshold, calculated from noise + level times noise factor. + + : param: X, time series (seismogram) + : type: `~obspy.core.stream.Stream` + + : param: pick, initial (AIC) P onset time + : type: float + + : param: TSNR, length of time windows around initial pick [s] + : type: tuple (T_noise, T_gap, T_signal) + + : param: minsiglength, minium required signal length [s] to + declare pick as P onset + : type: float + + : param: nfac, noise factor (nfac * noise level = threshold) + : type: float + + : param: minpercent, minimum required percentage of samples + above calculated threshold + : type: float + + : param: iplot, if iplot > 1, results are shown in figure + : type: int + ''' + + assert isinstance(X, Stream), "%s is not a stream object" % str(X) + + print 'Checking signal length ...' + + x = X[0].data + t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, + X[0].stats.delta) + + # generate envelope function from Hilbert transform + y = np.imag(sc.signal.hilbert(x)) + e = np.sqrt(np.power(x, 2) + np.power(y, 2)) + # get noise window + inoise = getnoisewin(t, pick, TSNR[0], TSNR[1]) + # get signal window + isignal = getsignalwin(t, pick, TSNR[2]) + # calculate minimum adjusted signal level + minsiglevel = max(e[inoise]) * nfac + # minimum adjusted number of samples over minimum signal level + minnum = len(isignal) * minpercent/100 + # get number of samples above minimum adjusted signal level + numoverthr = len(np.where(e[isignal] >= minsiglevel)[0]) + + if numoverthr >= minnum: + print 'checksignallength: Signal reached required length.' + returnflag = 1 + else: + print 'checksignallength: Signal shorter than required minimum signal length!' + print 'Presumably picked picked noise peak, pick is rejected!' + returnflag = 0 + + if iplot == 2: + plt.figure(iplot) + p1, = plt.plot(t,x, 'k') + p2, = plt.plot(t[inoise], e[inoise]) + p3, = plt.plot(t[isignal],e[isignal], 'r') + p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], \ + [minsiglevel, minsiglevel], 'g') + p5, = plt.plot([pick, pick], [min(x), max(x)], 'c') + plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \ + 'Envelope Signal Window', 'Minimum Signal Level', \ + 'Onset'], loc='best') + plt.xlabel('Time [s] since %s' % X[0].stats.starttime) + plt.ylabel('Counts') + plt.title('Check for Signal Length, Station %s' % X[0].stats.station) + plt.yticks([]) + plt.show() + raw_input() + plt.close(iplot) + + return returnflag + + From 017532272f6d13ba6f939c1e2cec342e8386a4e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 24 Jun 2015 14:17:43 +0200 Subject: [PATCH 6/6] Modified: new function checksignallength added. --- pylot/core/pick/run_autopicking.py | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/pylot/core/pick/run_autopicking.py b/pylot/core/pick/run_autopicking.py index 7bf63cef..65bf0514 100755 --- a/pylot/core/pick/run_autopicking.py +++ b/pylot/core/pick/run_autopicking.py @@ -74,6 +74,10 @@ def run_autopicking(wfstream, pickparam): minFMSNR = pickparam.getParam('minFMSNR') fmpickwin = pickparam.getParam('fmpickwin') minfmweight = pickparam.getParam('minfmweight') + # parameters for checking signal length + minsiglength = pickparam.getParam('minsiglength') + minpercent = pickparam.getParam('minpercent') + nfacsl = pickparam.getParam('noisefactor') # initialize output Pweight = 4 # weight for P onset @@ -94,6 +98,7 @@ def run_autopicking(wfstream, pickparam): aicSflag = 0 aicPflag = 0 + Pflag = 0 Sflag = 0 # split components @@ -152,9 +157,15 @@ def run_autopicking(wfstream, pickparam): # of class AutoPicking aicpick = AICPicker(aiccf, tsnrz, pickwinP, iplot, None, tsmoothP) ############################################################## + # check signal length to detect spuriously picked noise peaks + z_copy[0].data = tr_filt.data + Pflag = checksignallength(z_copy, aicpick.getpick(), tsnrz, minsiglength, \ + nfacsl, minpercent, iplot) + ############################################################## # go on with processing if AIC onset passes quality control if (aicpick.getSlope() >= minAICPslope and - aicpick.getSNR() >= minAICPSNR): + aicpick.getSNR() >= minAICPSNR and + Pflag == 1): aicPflag = 1 print 'AIC P-pick passes quality control: Slope: %f, SNR: %f' % \ (aicpick.getSlope(), aicpick.getSNR()) @@ -190,7 +201,7 @@ def run_autopicking(wfstream, pickparam): mpickP = refPpick.getpick() ############################################################# if mpickP is not None: - # quality assessment + # quality assessment # get earliest and latest possible pick and symmetrized uncertainty [lpickP, epickP, Perror] = earllatepicker(z_copy, nfacP, tsnrz, mpickP, iplot) @@ -227,8 +238,8 @@ def run_autopicking(wfstream, pickparam): Sflag = 0 else: - print 'run_autopicking: No vertical component data available, ' \ - 'skipping station!' + print 'run_autopicking: No vertical component data availabler!, ' \ + 'Skip station!' if edat is not None and ndat is not None and len(edat) > 0 and len( ndat) > 0 and Pweight < 4: