-
Notifications
You must be signed in to change notification settings - Fork 16
/
a1_ccf_ambnoise_P.m
538 lines (451 loc) · 18.8 KB
/
a1_ccf_ambnoise_P.m
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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
% Calculate ambient noise cross correlation record from multiple stationpairs
% for Vertical (Z) only using the methods from Bensen et al. (2007) GJI
% DOI:10.1111/j.1365-246X.2007.03374.x
%
% Expects files organized like so:
% {datadirectory}/{station}/{station}.{yyyy}.{jday}.{hh}.{mm}.{SS}.{COMP}.sac
% e.g.: mydata/CC05/CC05.2018.112.00.00.00.BDH.sac
%
% JBR, Jan 2020: Implemented frequency-time normalization after
% Shen et al. (2012) BSSA; DOI:10.1785/0120120023. This greatly improves signal
% extraction compared to typical one-bit noralization and whitening of Bensen et
% al. (2007) GJI. Faster FiltFiltM() can be replaced with MATLAB's slower
% built-in filtfilt().
%
% JBR, update: We have found that doing no time or frequency normalization at all
% can produce higher SNR traces than doing one-bit or time-frequency normalization.
% Therefore, the default is to use the raw seismograms as is without any preprocessing.
%
% (NOTE: FUNCTIONIZE IN THE FUTURE)
% Patty Lin -- 10/2014
% Natalie Accardo
% Josh Russell
% https://github.com/jbrussell
clear;
setup_parameters;
strSACcomp = 'BDH'; % 'Z'
strNAMEcomp = 'PP'; % 'ZZ'
IsFigure1 = 1;
IsFigure2 = 0;
% OUTPUT SETTINGS
IsOutputFullstack = 1; % Save full year ccf stacks
IsOutputMonthstack = 0; % save month ccf stacks
IsOutputDaystack = 0; % save day ccf stacks
IsOutputSinglestack = 0; % save single ccf before stacking
IsOutputSeismograms = 0; % save raw seismograms before cross-correlating
% GENERAL PROCESSING
IsRemoveIR = 0; % remove instrument response
units_RemoveIR = 'M'; % 'M' displacement | 'M/S' velocity
IsDetrend = 1; % detrend the data
IsTaper = 1; % Apply cosine taper to data chunks
%%%%%%%%%%% OPTIONS FOR PREPROCESSING %%%%%%%%%%%%
% (1) ONE-BIT NORMALIZATION & SPECTRAL WHITENING? (Bensen et al. 2007)
IsSpecWhiten = 0; % Whiten spectrum
IsOBN = 0; % One-bit normalization
% (2) TIME-FREQUENCY NORMALIZATION (Ekstrom et al. 2009; Shen et al. 2011)
IsFTN = 0; % Frequency-time normalization? (If 1, applied instead of whitening and one-bit normalization)
frange_FTN = [1/60 1/10]; % frequency range over which to construct FTN seismograms
% (3) BASIC PREFILTER (Ekstrom 2011)
IsPrefilter = 0; % apply butterworth bandpass filter before cross-correlation?
frange_prefilt = [1/100 1/10];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Setup parallel pool
% Nworkers = 4; % number of workers in pool for parallel processing
delete(gcp('nocreate'));
% parpool(Nworkers);
% input path
datadir = parameters.datapath;
PZpath = parameters.PZpath;
figpath = parameters.figpath;
seis_path = parameters.seis_path;
orientation_path = parameters.orientation_path;
dt = parameters.dt;
winlength = parameters.winlength;
year = ''; %'2012';
Nstart_sec = parameters.Nstart_sec; % (seconds) offset start of file
Nstart = Nstart_sec/dt; % Number of samples
comp = parameters.comp;
%dist_min = 20;
dist_min = parameters.mindist;
% Build File Structure: cross-correlations
ccf_path = parameters.ccfpath;
ccf_winlength_path = [ccf_path,'window',num2str(winlength),'hr/'];
ccf_singlestack_path = [ccf_winlength_path,'single/'];
ccf_daystack_path = [ccf_winlength_path,'dayStack/'];
ccf_monthstack_path = [ccf_winlength_path,'monthStack/'];
ccf_fullstack_path = [ccf_winlength_path,'fullStack/'];
if ~exist(ccf_path)
mkdir(ccf_path)
end
if ~exist(ccf_winlength_path)
mkdir(ccf_winlength_path)
end
if ~exist(ccf_singlestack_path)
mkdir(ccf_singlestack_path)
end
if ~exist(ccf_daystack_path)
mkdir(ccf_daystack_path)
end
if ~exist(ccf_monthstack_path)
mkdir(ccf_monthstack_path)
end
if ~exist(ccf_fullstack_path)
mkdir(ccf_fullstack_path)
end
PATHS = {ccf_singlestack_path; ccf_daystack_path; ccf_monthstack_path; ccf_fullstack_path};
for ipath = 1:length(PATHS)
ccfZ_path = [PATHS{ipath},'ccf',strNAMEcomp,'/'];
if ~exist(ccfZ_path)
mkdir(ccfZ_path);
end
end
% Build File Structure: figures
fig_winlength_path = [figpath,'window',num2str(winlength),'hr/'];
if ~exist(figpath)
mkdir(figpath);
end
if ~exist(fig_winlength_path)
mkdir(fig_winlength_path)
end
% Build File Structure: windowed seismograms
seis_winlength_path = [seis_path,'window',num2str(winlength),'hr/'];
if ~exist(seis_path)
mkdir(seis_path);
end
if ~exist(seis_winlength_path)
mkdir(seis_winlength_path)
end
warning('off','MATLAB:nargchk:deprecated')
%% ------------------- loop through center station station-------------------
stalist = parameters.stalist;
nsta=parameters.nsta; % number of target stations to calculate for
% READ OBS ORIENTATIONS
% [slist, orientations] = textread(orientation_path,'%s%f\n');
% Calculate filter coefficients for FTN
if IsFTN
[ b, a ] = get_filter_TFcoeffs( frange_FTN, dt );
end
for ista1=1:nsta
sta1=char(stalist(ista1,:));
% Build station directories
for ipath = 1:length(PATHS)
ccfZ_path = [PATHS{ipath},'ccf',strNAMEcomp,'/'];
if ~exist([ccfZ_path,sta1])
mkdir([ccfZ_path,sta1]);
end
end
seisZ_path = [seis_winlength_path,strNAMEcomp(1),'/'];
if ~exist([seisZ_path,sta1])
mkdir([seisZ_path,sta1]);
end
list1 = dir([datadir,sta1,'/*',strSACcomp,'.sac']);
for ista2=1:nsta
clear lat1 lat2 lon1 lon2 dist az baz vec_tz2 Z2raw vec_tz Z1raw
sta2=char(stalist(ista2,:));
% if same station, skip
if(strcmp(sta1,sta2))
continue
end
% check to see if we've already done this ccf
if exist([ccfZ_path,sta1,'/',sta1,'_',sta2,'_f.mat'])
display('CCF already exist, skip this pair');
continue
end
display(['performing cross-correlation for staion pair : ',sta1,' ', sta2]);
% -------------loop through each half day--------------------
nday_stack=0;
coh_sumZ = 0;
coh_num = 0;
% Get a list of all available data
ihday = 0;
month_counter = 0;
imonth = 0;
for ifil = 1:length(list1)
file1cZ = list1(ifil).name;
% Check that day file exists for station 2
Nchar = length(sta1);
file2cZ = dir([datadir,sta2,'/',sta2,file1cZ(Nchar+1:end)]);
str = strsplit(file1cZ,'.');
hdayid = [str{2},'.',str{3},'.',str{4},'.',str{5},'.',str{6}];
if isempty(file2cZ)
disp(['No data for ',sta2,' on day ',hdayid,'... skipping'])
continue
end
file2cZ = file2cZ.name;
if month_counter == 0
coh_sumZ_month = 0;
coh_num_month = 0;
end
clear data1cZ data2cZ
ihday = ihday +1;
month_counter = month_counter + 1;
clear temp
%temp = strsplit(daylist1(ihday).name,'.');
disp(['Looking at ',hdayid,' ',sta2]);
data1cZ= dir([datadir,sta1,'/',year,'/',sta1,'.',hdayid,'.*',strSACcomp,'.sac']);
data2cZ= dir([datadir,sta2,'/',year,'/',sta2,'.',hdayid,'.*',strSACcomp,'.sac']);
data1cZ = [datadir,sta1,'/',year,'/',data1cZ.name];
data2cZ = [datadir,sta2,'/',year,'/',data2cZ.name];
%------------------- TEST IF DATA EXIST------------------------
[S1Zt,S1Zraw,S1,S1Ztstart] = load_sac(data1cZ);
[S2Zt,S2Zraw,S2,S2Ztstart] = load_sac(data2cZ);
% Check that sample rates are the same
if S1.DELTA ~= S2.DELTA
error('S1 and S2 sample rates don''t match!');
end
% Make sure all times are relative to same reference point
starttime = S1Ztstart;
S1Zt = S1Zt + seconds(S1Ztstart-starttime);
S2Zt = S2Zt + seconds(S2Ztstart-starttime);
% Ensure that files have same start time to within 1 sample
if abs(seconds(S1Ztstart-S2Ztstart)) > S1.DELTA
error('Station files do not have same start time');
end
% Make sure sample rates all match
if (abs(S1.DELTA-dt) >= 0.01*dt ) || (abs(S2.DELTA-dt) >= 0.01*dt )
error('sampling interval does not match data! check dt');
end
%------------------- Remove instrument response ------------------------
if IsRemoveIR
pzfile1 = dir([PZpath,'/RESP.*.',sta1,'.*.*Z']); % PZ for H1 and H2 are identical
pzfile2 = dir([PZpath,'/RESP.*.',sta2,'.*.*Z']);
% Read RESP file for station 1
[z,p,c,units] = read_sac_RESP([PZpath,pzfile1.name],units_RemoveIR);
dt1 = abs(S1Zt(1)-S1Zt(2));
dt2 = abs(S2Zt(1)-S2Zt(2));
% Remove instrument response for station 1 Z
S1Zraw = rm_resp(S1Zraw,z,p,c,dt1);
% Read RESP file for station 2
[z,p,c,units] = read_sac_RESP([PZpath,pzfile2.name],units_RemoveIR);
% Remove instrument response for station 2 Z
S2Zraw = rm_resp(S2Zraw,z,p,c,dt2);
end
% Check to make sure there's actual data
zeroind1 = find(S1Zraw == 0);
zeroind2 = find(S2Zraw == 0);
if length(zeroind1) == length(S1Zraw) || length(zeroind2) == length(S2Zraw)
disp('All zeros!');
continue
end
if(length(S1Zt)*length(S2Zt)==0)
display(['no data for ! station ',sta2]);
continue
end
% Determine the time span to cut to ... this will change with
% different segments
clear tcut
minT1Z = min(S1Zt);
minT2Z = min(S2Zt);
if length(S1Zraw) < 30000
disp(['Sta1 ',sta1,' : ',num2str(length(S1Zraw)),' is too short!'])
continue
elseif length(S2Zraw) < 30000
disp(['Sta2 ',sta2,' : ',num2str(length(S2Zraw)),' is too short!'])
continue
end
if(~exist('lat2','var'));
lat1=S1.STLA;
lon1=S1.STLO;
dep1=S1.STEL; % depth is negative for OBS and positive for land stations
lat2=S2.STLA;
lon2=S2.STLO;
dep2=S2.STEL; % depth is negative for OBS and positive for land stations
% Get the interstation distance and azimuth
[delta,S1az]=distance(lat1,lon1,lat2,lon2);
[delta,S2az]=distance(lat2,lon2,lat1,lon1);
dist=deg2km(delta);
if(dist < dist_min)
display(['distance shorter than ',num2str(dist_min),' km, skip']);
break
end
end % if lat variabls
stapairsinfo.stanames = {sta1,sta2};
stapairsinfo.lats = [lat1,lat2];
stapairsinfo.lons = [lon1,lon2];
stapairsinfo.dt = dt;
stapairsinfo.r = dist;
% % Frequency-time normalization
% if IsFTN
% [ S1Zraw ] = FTN( S1Zraw, frange_FTN, dt );
% [ S2Zraw ] = FTN( S2Zraw, frange_FTN, dt );
% [ S1H1raw ] = FTN( S1H1raw, frange_FTN, dt );
% [ S2H1raw ] = FTN( S2H1raw, frange_FTN, dt );
% [ S1H2raw ] = FTN( S1H2raw, frange_FTN, dt );
% [ S2H2raw ] = FTN( S2H2raw, frange_FTN, dt );
% end
% START WINDOWING
hour_length = winlength;
nwin = floor(24/hour_length)*2-1; %
win_length = hour_length*60*60/dt; % length of individual windows.
win_start = 1;
coh_sumZ_day = 0;
coh_num_day = 0;
last_pt = win_length*.5*(nwin-1)+1+Nstart+win_length;
if last_pt < length(S1Zraw)
nwin = nwin + 1;
end
% tic
parfor iwin = 1:nwin
% clear tcut S1R S2R S1T S2T S1Z S2Z fftS1R fftS2R fftS1T fftS2T fftS1Z fftS2Z
% cut in time
if hour_length == 24
pts_begin = Nstart;
pts_end = length(S1Zraw)-Nstart;
else
pts_begin = win_length*.5*(iwin-1)+1+Nstart;
pts_end = pts_begin+win_length;
end
if pts_begin > length(S1Zraw) || pts_begin > length(S2Zraw) || pts_end > length(S1Zraw) || pts_end > length(S2Zraw)
% disp('(Z) Points greater than the data... fixing window');
pts_begin = length(S1Zraw)-win_length-Nstart;
pts_end = pts_begin+win_length;
%continue
end
tcut = [pts_begin:pts_end] * dt;
% cut in tim Z for STA1
S1Z=interp1(S1Zt,S1Zraw,tcut);
S1Z(isnan(S1Z))=0;
% cut in tim Z for STA2
S2Z=interp1(S2Zt,S2Zraw,tcut);
S2Z(isnan(S2Z))=0;
%detrend
if IsDetrend
S1Z = detrend(S1Z);
S2Z = detrend(S2Z);
end
% Apply cosine taper
if IsTaper
S1Z = cos_taper(S1Z);
S2Z = cos_taper(S2Z);
end
% Apply prefilter
if IsPrefilter
[bb,aa] = butter(2,frange_prefilt*2*dt);
S1Z = FiltFiltM(bb,aa,S1Z);
S2Z = FiltFiltM(bb,aa,S2Z);
end
if IsFigure2
figure(49)
clf
%Z
subplot(5,1,1)
plot(tcut,S1Z,'-k')
%ylim([-0.15e-5 0.15e-5])
xlim([tcut(1) tcut(end)])
title(strNAMEcomp(1));
hold on
pause;
%return
end
%-------------------- Vertical Component --------------
if IsFTN
% Frequency-time normalization
[ S1Z ] = FTN( S1Z, b, a );
[ S2Z ] = FTN( S2Z, b, a );
fftS1Z = fft(S1Z);
fftS2Z = fft(S2Z);
else
% One-bit normalization
if IsOBN
S1Z = runwin_norm(S1Z);
S2Z = runwin_norm(S2Z);
end
%fft
fftS1Z = fft(S1Z);
fftS2Z = fft(S2Z);
%Whiten
if IsSpecWhiten
fftS1Z = spectrumwhiten_smooth(fftS1Z,0.001);
fftS2Z = spectrumwhiten_smooth(fftS2Z,0.001);
end
end
% calcaulate daily CCF and stack for radial (1-->2 Causal; 2-->1 Acausal)
coh_trace = fftS2Z .* conj(fftS1Z);
coh_trace = coh_trace ./ abs(fftS1Z) ./ abs(fftS2Z);
coh_trace(isnan(coh_trace)) = 0;
coh_sumZ = coh_sumZ + coh_trace;
coh_trace_Z = coh_trace;
coh_sumZ_day = coh_sumZ_day + coh_trace;
coh_sumZ_month = coh_sumZ_month + coh_trace;
% coh_num = coh_num + 1;
coh_num_day = coh_num_day + 1;
coh_num_month = coh_num_month + 1;
% toc
if IsOutputSinglestack % save individual xcor
ccfZ_singlestack_path = [ccf_singlestack_path,'ccf',strNAMEcomp,'/'];
save(sprintf('%s%s/%s_%s_%d_f.mat',ccfZ_singlestack_path,sta1,sta1,sta2,coh_num),'coh_trace_Z','stapairsinfo');
end
if IsOutputSeismograms % save seismograms before xcor
seisZ_path = [seis_winlength_path,strNAMEcomp(1),'/'];
save(sprintf('%s%s/%s_%d_f.mat',seisZ_path,sta1,sta1,coh_num),'S1Z','stapairsinfo');
end
end % end window
% toc
coh_num = coh_num + coh_num_day;
if IsOutputDaystack
% Save day stack
daystr = datestr(starttime,'YYYYmmddHHMMSS');
ccfZ_daystack_path = [ccf_daystack_path,'ccf',strNAMEcomp,'/'];
clear coh_sum
coh_sum = coh_sumZ_day;
save(sprintf('%s%s/%s_%s_%s_f.mat',ccfZ_daystack_path,sta1,sta1,sta2,daystr),'coh_sum','coh_num_day','stapairsinfo','starttime');
end
if IsOutputMonthstack
% Save 30 day (month) stack
if month_counter == 30
imonth = imonth + 1;
ccfZ_monthstack_path = [ccf_monthstack_path,'ccf',strNAMEcomp,'/'];
clear coh_sum
coh_sum = coh_sumZ_month;
save(sprintf('%s%s/%s_%s_month%d_f.mat',ccfZ_monthstack_path,sta1,sta1,sta2,imonth),'coh_sum','coh_num_month','stapairsinfo');
month_counter = 0; % start over
end
end
end % end hday
if coh_num > 1
if IsFigure1
f101 = figure(101);clf;
% set(gcf,'position',[400 400 600 300]);
subplot(3,1,2)
T = length(coh_sumZ);
faxis = [0:(T-mod(T-1,2))/2 , -(T-mod(T,2))/2:-1]/dt/T;
ind = find(faxis>0);
plot(faxis(ind),smooth(real(coh_sumZ(ind)/coh_num),100));
title(sprintf('%s %s coherency %s ,station distance: %f km',sta1,sta2,strNAMEcomp(1),dist));
% xlim([0.01 1/(dt*2.5)])
xlim([0.01 0.5])
%xlim([0.04 0.16])
xlabel('Frequency')
subplot(3,1,3)
costap_wid = 0.2;
coperiod = 1./[0.5 0.01];
ccf = coh_sumZ ./ coh_num;
N = length(ccf);
ccf_ifft = real(ifft(ccf,N)); % inverse FFT to get time domain
ccf_ifft = fftshift(ccf_ifft); % rearrange values as [-lag lag]
ccf_ifft = detrend(ccf_ifft);
ccf_ifft = cos_taper(ccf_ifft);
[ ccf_filtered ] = tukey_filt( fft(fftshift(ccf_ifft)),coperiod,dt,costap_wid );
ccf_ifft = fftshift(real(ifft(ccf_filtered)));
time = ([0:N-1]-floor(N/2))*dt;
time = [time(time<0), time(time>=0)];
plot(time,ccf_ifft,'-r');
phv_min_win = 0.8; % km/s
% xlim([-1 1]*max([stapairsinfo.r/phv_min_win,100]));
xlim([-250 250])
%xlim([0.04 0.16])
xlabel('Lag Time')
drawnow
print(f101,'-dpng',[fig_winlength_path,sta1,'_',sta2,'_',strNAMEcomp,'.png']);
%pause;
end
if IsOutputFullstack
ccfZ_fullstack_path = [ccf_fullstack_path,'ccf',strNAMEcomp,'/'];
clear coh_sum
coh_sum = coh_sumZ;
save(sprintf('%s%s/%s_%s_f.mat',ccfZ_fullstack_path,sta1,sta1,sta2),'coh_sum','coh_num','stapairsinfo');
end
end
end % ista2
end % ista1
delete(gcp('nocreate')); % remove parallel pools