tcs-ah/epos-ip-scripts
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added TemplateMatchingDetection application

Browse Source
This commit is contained in:
Joanna Kocot 2019-04-05 13:31:09 +02:00
parent 0b494bba22
commit 75e691de5f
9 changed files with 669 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/StressMaps/LICENCE.txt
View File

@ -1,3 +1,3 @@
This application is a part of the IS-EPOS e-PLATFORM This application is a part of the IS-EPOS e-PLATFORM
Creative Commons Attribution-ShareAlike 4.0 International License: https://creativecommons.org/licenses/by-sa/4.0/ This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License: https://creativecommons.org/licenses/by-sa/4.0/

15
src/TemplateMatchingDetection/LICENCE.txt Normal file
View File

@ -0,0 +1,15 @@
This application is a part of the IS-EPOS e-PLATFORM
Copyright 2019 Olivier Lengliné
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

83
src/TemplateMatchingDetection/TM_EPOS.m Normal file
View File

@ -0,0 +1,83 @@
function [T, yc] = TM_EPOS(pickinfo,event_wfm,continuousfile,nwin,npre,fmin,fmax, ...
min_dist,min_cor)
% Check for integers
if(floor(nwin) ~= nwin); error('nwin is not an integer'); end
if(floor(npre) ~= npre); error('npre is not an integer'); end
if(floor(min_dist) ~= min_dist); error('min_dist is not an integer'); end
% Check frequencies
if(fmin>fmax); error('fmin > fmax'); end
%% Parameters definition (example)
% % Filenames
% pickinfo = 'picks.xml'; % Picking files
% event_wfm= 'VN.TBVB..HHE.300000.SAC'; % Event waveform
% continuousfile= 'VN.TBVB..HHE.620000.SAC'; % Continuous waveform
%
% % Set Template window parameters
% nwin = 512; % Total number of samples in the template window
% npre = 50; % Number of points before the arrival time
%
% % Set filter parameters
% fmin = 1; % Minimum frequency
% fmax = 15; % Maximum frequency
%
% % Set parameters for declaring a detection
% min_dist = 200; % Minimum required distance (in samples) between 2 detections
% min_cor = 0.7; % Correlation threshold for declaring a detection (optional)
% Get the P-wave pick info
[tp,sta,ntwk,channel]=read_xml(pickinfo);
% Read the continuous sac file
Fc=readsac(continuousfile);
if(Fc.tau==-1); error('The continuous sac file does not exist'); end
% And we will filter it.
fs = round(1./Fc.delta) ; fn = fs/2;
[b,a] = butter(4,[fmin fmax]./fn);
yc = Fc.trace - mean(Fc.trace); yc = filtfilt(b,a,yc);
% Read the event trace file
Fe=readsac(event_wfm);
if(Fe.tau==-1); error('The event sac file does not exist'); end
% Check if the two waveforms correspond to the same station, channel, and
% have the same sampling frequency
if(Fe.delta ~= Fc.delta); error('Sampling frequencies are not the same'); end
if(Fe.kcmpnm ~= Fc.kcmpnm); error('Components are not the same'); end
if(Fe.kstnm ~= Fc.kstnm); error('Stations are not the same'); end
% Check if picks correspond to the supplied waveform event file
if(strcmp(Fe.kstnm,sta))
if(strcmp(Fe.kcmpnm,channel))
if(strcmp(Fe.knetwk,ntwk))
[mo,day]=jd2md(Fe.nzjday,Fe.nzyear);
t0 = datenum(Fe.nzyear,mo,day,Fe.nzhour,Fe.nzmin,Fe.sec);
t1 = t0 + (Fe.npts*Fe.delta)/(24*3600);
% Check if the picking time is within the event waveform
% file. It should be the case but just to be sure.
if( (tp > t0) && (tp < t1))
% Extract the template waveform window
y=mk_template(Fe.trace,fmin,fmax,tp,t0,1./Fe.delta,nwin,npre);
% Compute the correlation
R=compute_correlation(yc,y);
% Find detections
T= find_events(R,min_dist,min_cor);
else
error('Picking time not within the event waveform duration')
end
else
error('Not the same Network')
end
else
error('Not the same Channel')
end
else
error('Not the same station')
end

68
src/TemplateMatchingDetection/compute_correlation.m Normal file
View File

@ -0,0 +1,68 @@
function R=compute_correlation(v,H)
% where H is the template waveform
% and v is the continuous signal
% we require that both signals are filtered similarly
% Get fourier transform of the continuous signal
[e,nv] = size(v);
if(nv ==1); v = v'; nv = e; end % Make sure we have the good shape
f_v = fft(v);
% Transform template signal
[nh1,nh2] = size(H); % Get size
if(nh2 ==1); H = H'; nh2 = nh1; end % Make sure we have the good shape
H = H - mean(H); % remove mean
H = H.*tukeywin(nh2,0.1)'; % Tapering
H = H./std(H); % Normalize by standard deviation
npad = nv - nh2;
H = padarray(H,[0 npad],0,'post');
f_H = fft(H); % Get the Fourier transform of H (flipped and same size
% as v.
% Compute the standard deviation of continuous signal
s_v= std_v(v,nv,nh2);
% Compute correlation
R = fliplr(real(ifft(f_H.*conj(f_v))));
% Get normalized values
R = R./s_v;
end
function s_v= std_v(v,nv,nw)
% This function computes the stabdard deviation on the continuous signal
% on windows of size nw
% Compute std for the first window
s1 = sum(v(1:nw));
s2 = v(1:nw)*v(1:nw)';
s_v(1) = sqrt( (s2-s1*s1/nw)/(nw-1));
% Initialize all values with this one
s_v(1:nv) = s_v(1);
% Increment
for i = 2:nv-nw+1
s1 = s1 - v(i-1) + v(i+nw-1);
s2 = s2 - v(i-1)*v(i-1) + v(i+nw-1)*v(i+nw-1);
s_v(i ) = sqrt( (s2-s1*s1/nw)/(nw-1));
end
% For the end of file
for i = nv-nw+2:nv
s1 = s1 - v(i-1);
s2 = s2 - v(i-1)*v(i-1);
s_v(i ) = sqrt( (s2-s1*s1/(nv-i))/(nv-i-1));
end
s_v = s_v * nw;
end

19
src/TemplateMatchingDetection/find_events.m Normal file
View File

@ -0,0 +1,19 @@
function T= find_events(R,min_dist,min_cor)
% Take the enveloppe of the correlation signal
% (not very much used here but it s betetr in case we stack several
% correlation functions
R = abs(hilbert(R));
% Set the correrlation threshold based on the median absolute devaition if
% it has not been set
if nargin == 2
min_cor = mean(R) + 8 * mad(R);
end
% Build the time vector
t = 1:length(R);
% Find Peaks in the correlation function separated at least by min_dist and
% whose values are higher than min_cor
[~,T]=findpeaks(R,t,'MinpeakDistance',min_dist,'MinpeakHeight',min_cor);

38
src/TemplateMatchingDetection/jd2md.m Normal file
View File

@ -0,0 +1,38 @@
function [mo,day]=jd2md(jd,yr)
% function [mo,day]=jd2md(jd,<yr>)
%
% Function to transfer Julian day to date with the option to have leap years
% Input 'leap' or year number for <yr> if you want to use leap years
% Default is regular year
%
% Remark: a leap year is a year with 366 days: every 4th year is a leap year
% BUT every 100th is NOT _AND_ every 400th _IS_ again a leap year...
%
% György Hetényi
% 22 Nov 2004
% ------------------------------------------------------------------
for i=1:length(jd)
if nargin<2 yr(i)=1; end
if ischar(yr(i))==1 & yr(i)=='leap' type='leap';
elseif mod(yr(i),4)==0 & mod(yr(i),100)~=0 type='leap';
elseif mod(yr(i),400)==0 type='leap';
else type='regu';
end
clear mos
if type=='leap'
mos=[31,29,31,30,31,30,31,31,30,31,30,31];
end
if type=='regu'
mos=[31,28,31,30,31,30,31,31,30,31,30,31];
end
k=1;
while jd(i)>sum(mos(1:k))
k=k+1;
end
mo(i)=k;
day(i)=jd(i)-sum(mos(1:mo(i)-1));
end

17
src/TemplateMatchingDetection/mk_template.m Normal file
View File

@ -0,0 +1,17 @@
function y=mk_template(y,fmin,fmax,tp,t0,fs,nwin,npre)
% This function build the template waveform window
% Filter the whole seismogram (to avoid edge problem)
fn = fs/2;
[b,a] = butter(4,[fmin fmax]./fn);
y = y - mean(y);
y = filtfilt(b,a,y);
% Get the first point of the seismogram for windowing P-wave
I = round((tp-t0)*24*3600*fs);
I = I -npre;
y = y(I:I+nwin-1);
y = y .* tukeywin(nwin,0.05);

141
src/TemplateMatchingDetection/read_xml.m Normal file
View File

@ -0,0 +1,141 @@
function [tp,sta,ntwk,channel]=read_xml(xmlfilename)
% Load the xml File architecture
F = parseXML(xmlfilename);
k1=xml_field(F,'eventParameters'); H = F.Children(k1);
k2=xml_field(H,'event'); H = H.Children(k2);
k3=xml_field(H,'pick');
% Several picks are possible
for ik = 1:length(k3)
H_test = H.Children(k3(ik));
% Test if we have the P wave
k=xml_field(H_test,'phaseHint');
% If this is the P-wave read the informations
if(strcmp(H_test.Children(k).Children(1).Data,'P'))
k=xml_field(H_test,'waveformID'); A = H_test.Children(k);
[sta,ntwk,channel] = read_xml_attributes(A);
k=xml_field(H_test,'time'); G=H_test.Children(k);
k=xml_field(G,'value'); G = G.Children(k);
tp = datenum(G.Children(1).Data,'yyyy-mm-ddTHH:MM:SS.FFF');
end
end
end
% Read waveform attributes
function [sta,ntwk,channel] = read_xml_attributes(A)
na = length(A.Attributes);
for ia = 1:na
if(strcmp(A.Attributes(ia).Name,'channelCode'))
channel = A.Attributes(ia).Value;
end
if(strcmp(A.Attributes(ia).Name,'networkCode'))
ntwk = A.Attributes(ia).Value;
end
if(strcmp(A.Attributes(ia).Name,'stationCode'))
sta = A.Attributes(ia).Value;
end
end
end
% Get the current xml field with the name keystring
function [ikey]=xml_field(F,keystring)
nchilds = length(F.Children);
ikey = [];
for i = 1:nchilds
if(strcmp(F.Children(i).Name,keystring))
ikey = [ikey i];
end
end
nmatch = length(ikey);
if(nmatch == 0)
error('No field %s', keystring)
end
end
function theStruct = parseXML(filename)
% PARSEXML Convert XML file to a MATLAB structure.
% THis functions is from MATLAB xmlread doc
try
tree = xmlread(filename);
catch
error('Failed to read XML file %s.',filename);
end
% Recurse over child nodes. This could run into problems
% with very deeply nested trees.
try
theStruct = parseChildNodes(tree);
catch
error('Unable to parse XML file %s.',filename);
end
end
% ----- Local function PARSECHILDNODES -----
function children = parseChildNodes(theNode)
% Recurse over node children.
children = [];
if theNode.hasChildNodes
childNodes = theNode.getChildNodes;
numChildNodes = childNodes.getLength;
allocCell = cell(1, numChildNodes);
children = struct( ...
'Name', allocCell, 'Attributes', allocCell, ...
'Data', allocCell, 'Children', allocCell);
for count = 1:numChildNodes
theChild = childNodes.item(count-1);
children(count) = makeStructFromNode(theChild);
end
end
end
% ----- Local function MAKESTRUCTFROMNODE -----
function nodeStruct = makeStructFromNode(theNode)
% Create structure of node info.
nodeStruct = struct( ...
'Name', char(theNode.getNodeName), ...
'Attributes', parseAttributes(theNode), ...
'Data', '', ...
'Children', parseChildNodes(theNode));
if any(strcmp(methods(theNode), 'getData'))
nodeStruct.Data = char(theNode.getData);
else
nodeStruct.Data = '';
end
end
% ----- Local function PARSEATTRIBUTES -----
function attributes = parseAttributes(theNode)
% Create attributes structure.
attributes = [];
if theNode.hasAttributes
theAttributes = theNode.getAttributes;
numAttributes = theAttributes.getLength;
allocCell = cell(1, numAttributes);
attributes = struct('Name', allocCell, 'Value', ...
allocCell);
for count = 1:numAttributes
attrib = theAttributes.item(count-1);
attributes(count).Name = char(attrib.getName);
attributes(count).Value = char(attrib.getValue);
end
end
end

287
src/TemplateMatchingDetection/readsac.m Normal file
View File

@ -0,0 +1,287 @@
function F=readsac(files)
% F=readsac('files')
%
% Read a SAC-file or a set of SAC-files.
%
% Input :
% 'files' corresponds either to the complete name of one SAC-file,
% either to a set of SAC-files with the logical expression *.
% Output :
% F is a structure (length equal to the number of SAC-files read).
% It contains ALL the fields of the SAC's header and the trace in
% the field "trace".
%
% EXAMPLE :
% F=readsac('2004.02.23-17.31.00.ABH.SHZ.SAC')
% reads the file 2004.02.23-17.31.00.ABH.SHZ.SAC and saves it into
% the structure F (dimension: 1*1).
% F=readsac('Data/*BHZ*.SAC')
% reads all the files of the form *BHZ*.SAC in the directory "Data".
% The size of F equals the number of these files.
%
% From J. Vergne and G. Hetenyi
%------------------------------------------------------------------
% By default, the signals are read in mode 'r'
modlect='r';
% Determine the type of "files"
rep_files=fileparts(files);
list_files=dir(files);
if length(list_files)==0
% disp('"readsac": File(s) do not exist');
F.tau=-1;
else
for ifile=1:length(list_files)
nomfich=fullfile(rep_files,list_files(ifile).name);
clear dsac;
% Read signals
% ------------
% Following the signal type, reading procedure can be different:
% 'l' - IEEE floating point with little-endian byte ordering
% 'b' - IEEE floating point with big-endian byte ordering
[C,MAXSIZE,ENDIAN]=computer;
if ENDIAN=='L' bool_l=0;
else bool_l=1;
end
fidev=fopen(nomfich,modlect,'l');
if fidev > 0
h1=fread(fidev,70,'float'); % --------------
h2=fread(fidev,40,'long'); % reading header
h3=fread(fidev,192,'uchar'); % --------------
% testing byte-order, h2(7) must! be 6.
if h2(7)~=6
bool_l=1;
fclose(fidev);
fidev=fopen(nomfich,modlect,'b');
h1=fread(fidev,70,'float');
h2=fread(fidev,40,'long');
h3=fread(fidev,192,'uchar');
end
% reading trace
tamp=fread(fidev,inf,'float');
dsac.h1=h1;
dsac.h2=h2;
dsac.h3=h3;
% PART 1: reading float-type parameters
% ------------------------------------------
% %- comment are from original version,
% % comments ares from SAC-homepage
dsac.delta=h1(1); %- sampling time interval
dsac.depmin=h1(2); % minimum value of dependent variable
dsac.depmax=h1(3); % maximum value of dependent variable
dsac.scale=h1(4); % multiplying scale factor for dependent variable (not currently used)
dsac.odelta=h1(5); % observed increment if different from nominal value
dsac.b=h1(6); %- begin time (d<EFBFBD>calage du 1er <EFBFBD>chantillon) (beginning value of independent variable)
dsac.e=h1(7); % ending value of independent variable
dsac.o=h1(8); %- event origin time (seconds relative to reference time)
dsac.a=h1(9); % first arrival time (seconds relative to reference time)
dsac.internal10=h1(10); % INTERNAL
dsac.t0=h1(11); %- user defined time picks or markers
dsac.t1=h1(12); %- (seconds relative to reference time)
dsac.t2=h1(13); %-
dsac.t3=h1(14); %-
dsac.t4=h1(15); %-
dsac.t5=h1(16); %
dsac.t6=h1(17); %
dsac.t7=h1(18); %
dsac.t8=h1(19); %
dsac.t9=h1(20); %
dsac.f=h1(21); % fini or end of event time (seconds relative to reference time)
dsac.resp0=h1(22); % instrument response parameters (not currently used)
dsac.resp1=h1(23); %
dsac.resp2=h1(24); %
dsac.resp3=h1(25); %
dsac.resp4=h1(26); %
dsac.resp5=h1(27); %
dsac.resp6=h1(28); %
dsac.resp7=h1(29); %
dsac.resp8=h1(30); %
dsac.resp9=h1(31); %
dsac.stla=h1(32); %- station latitude (degrees, north positive)
dsac.stlo=h1(33); %- station longitude (degrees, east positive)
dsac.stel=h1(34); %- station elevation (meters)
dsac.stdp=h1(35); % station depth below surface (meters)(not currently used)
dsac.evla=h1(36); %- event latitude (degrees, north positive)
dsac.evlo=h1(37); %- event longitude (degrees, east positive)
dsac.evel=h1(38); % event elevation (meters)(not currently used)
dsac.evdp=h1(39); %- event depth below surface (meters)
dsac.mag=h1(40); % event magnitude
dsac.user0=h1(41); %- used defined variable storage area, floating!
dsac.user1=h1(42); %-
dsac.user2=h1(43); %-
dsac.user3=h1(44); %
dsac.user4=h1(45); %
dsac.user5=h1(46); %
dsac.user6=h1(47); %
dsac.user7=h1(48); %
dsac.user8=h1(49); %
dsac.user9=h1(50); %
dsac.dist=h1(51); %- station to event distance (km)
dsac.az=h1(52); %- event to station azimuth (degrees)
dsac.baz=h1(53); %- station to event azimuth (degrees)
dsac.gcarc=h1(54); %- station to event great circle arc length (degrees)
dsac.internal55=h1(55); % INTERNAL
dsac.internal56=h1(56); % INTERNAL
dsac.depmen=h1(57); % mean value of dependent variable
dsac.cmpaz=h1(58); %- component azimuth (degrees clockwise from north)
dsac.cmpinc=h1(59); %- component incidence angle (degrees from vertical)
dsac.xminimum=h1(60); % minimum value of X (spectral files only)
dsac.xmaximum=h1(61); % maximum value of X (spectral files only)
dsac.yminimum=h1(62); % minimum value of Y (spectral files only)
dsac.ymaximum=h1(63); % maximum value of Y (spectral files only)
dsac.unused64=h1(64); % UNUSED
dsac.unused65=h1(65); % UNUSED
dsac.unused66=h1(66); % UNUSED
dsac.unused67=h1(67); % UNUSED
dsac.unused68=h1(68); % UNUSED
dsac.unused69=h1(69); % UNUSED
dsac.unused70=h1(70); % UNUSED
% PART 2: reading long-type parameters
% ------------------------------------
% GMT time corresponding to reference (0) time in file
dsac.nzyear=h2(1); %- year
dsac.nzjday=h2(2); %- julian day
dsac.nzhour=h2(3); %- hour
dsac.nzmin=h2(4); %- minute
dsac.nzsec=h2(5); % second
dsac.nzmsec=h2(6); % millisecond
dsac.sec=h2(5)+h2(6)/1000; %- full second (NOT an original SAC-attribute!)
dsac.nvhdr=h2(7); % header version number: 6!
dsac.norid=h2(8); % origin ID (CSS 3.0)
dsac.nevid=h2(9); % event ID (CSS 3.0)
dsac.npts=h2(10); %- number of points per data component
dsac.internal81=h2(11);% INTERNAL
dsac.nwfid=h2(12); % waveform ID (CSS 3.0)
dsac.nxsize=h2(13); % spectral length (spectral files only)
dsac.nysize=h2(14); % spectral width (spectral files only)
dsac.unused85=h2(15); % UNUSED
dsac.iftype=h2(16); % type of file (required)(see SAC-page)
dsac.idep=h2(17); % type of dependent variable (see SAC-page)
dsac.iztype=h2(18); % reference time equivalence (see SAC-page)
dsac.unused89=h2(19); % UNUSED
dsac.iinst=h2(20); % type of recording instrument (not currently used)
dsac.istreg=h2(21); % station geogrpahic region (not currently used)
dsac.ievreg=h2(22); % event geographic location (not currently used)
dsac.ievtyp=h2(23); % type of event (see SAC-page)
dsac.iqual=h2(24); % quality of data (not currently used)(see SAC-page)
dsac.isynth=h2(25); % synthetic data flag (not currently used)
dsac.imagtyp=h2(26); % magnitude type
dsac.imagsrc=h2(27); % source of magnitude information
dsac.unused98=h2(28); % UNUSED
dsac.unused99=h2(29); % UNUSED
dsac.unused100=h2(30); % UNUSED
dsac.unused101=h2(31); % UNUSED
dsac.unused102=h2(32); % UNUSED
dsac.unused103=h2(33); % UNUSED
dsac.unused104=h2(34); % UNUSED
dsac.unused105=h2(35); % UNUSED
dsac.leven=h2(36); % TRUE if data is evenly spaced
dsac.lspol=h2(37); % TRUE if station components have positive polarity (left-hand rule)
dsac.lovrok=h2(38); % TRUE if it is okay to overwrite this file on disk
dsac.lcalda=h2(39); % TRUE if DIST,AZ,BAZ and GCARC are to be calculated form station and event coordinates
dsac.unused110=h2(40); % UNUSED
% PART 3: reading uchar-type parameters
% -------------------------------------
imin1=min(find(h3(1:8)==0 | h3(1:8)==32));
imin2=min(find(h3(9:24)==0 | h3(9:24)==32));
dsac.kstnm=rm_blanc(h3(1:1+imin1-1))'; %- station name
dsac.kevnm=rm_blanc(h3(9:9+imin2-1))'; %- event name
dsac.khole=rm_blanc(h3(25:32))'; % hole identification if nuclear event
dsac.ko=rm_blanc(h3(33:40))'; % event origin time identification
dsac.ka=rm_blanc(h3(41:48))'; % first arrival time identification
dsac.kt0=rm_blanc(h3(49:56))'; %- user defined time pick identifications, h1(11:20)
dsac.kt1=rm_blanc(h3(57:64))'; %-
dsac.kt2=rm_blanc(h3(65:72))'; %-
dsac.kt3=rm_blanc(h3(73:80))'; %-
dsac.kt4=rm_blanc(h3(81:88))'; %-
dsac.kt5=rm_blanc(h3(89:96))'; %
dsac.kt6=rm_blanc(h3(97:104))'; %
dsac.kt7=rm_blanc(h3(105:112))'; %
dsac.kt8=rm_blanc(h3(113:120))'; %
dsac.kt9=rm_blanc(h3(121:128))'; %
dsac.kf=rm_blanc(h3(129:136))'; % fini identification
dsac.kuser0=rm_blanc(h3(137:144))'; %- user defined variable storage area
dsac.kuser1=rm_blanc(h3(145:152))'; %-
dsac.kuser2=rm_blanc(h3(153:160))'; %-
dsac.kcmpnm=rm_blanc(h3(161:168))'; %- component name
dsac.knetwk=rm_blanc(h3(169:176))'; % name of seismic network
dsac.kdatrd=rm_blanc(h3(177:184))'; % date data was read onto computer
dsac.kinst=rm_blanc(h3(185:192))'; %- generic name of recording instrument
% reading trace
% -------------
dsac.trace=tamp;
dsac.tau = 1; % Added to check if the file exist
fclose(fidev);
else
%disp(['"readsac": file ' nomfich ' : reading error - file does not exist'])
dsac.tau=-1;
end
F(ifile)=dsac;
end
end
% -------------------------------------------------------------
function ch1=rm_blanc(ch)
% looks for whitespace in 'ch' and removes them
if ischar(ch)
ch1=ch(find(double(ch)~=32) & double(ch)~=0);
else
ch1=ch(find(ch~=32 & ch~=0));
ch1=char(ch1);
end