-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathget_coseis.m
executable file
·209 lines (192 loc) · 4.77 KB
/
get_coseis.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
function get_coseis
%DMM 10/2011
%
%Simulate metrics to extract real-time coseismic offsets
dir='/diego-local/Research/Data/El Mayor'
cd(dir)
load d_all_reduced_300.mat
% dir='/diego-local/Research/Data/Toki/';
% cd(dir);
% load toki.mat
%Run parameters
Tavg=20; %Average for final offset
thresh=0.015; %minimum horizontal offset
p=0.25; %Percent of peak variance for detection
Tvar=20; %LAgging variance window
%
% N=coseis.N(:,22:422);
% E=coseis.E(:,22:422);
% U=coseis.U(:,22:422);
% T=coseis.T(:,22:422);
N=coseis.N;
E=coseis.E;
U=coseis.U;
T=coseis.T;
H=sqrt(E.^2+N.^2);
mH=max(H');
i=find(mH>thresh);
H=H(i,:);
E=E(i,:);
N=N(i,:);
U=U(i,:);
T=T(i,:);
%Update old data
coseis.lat=coseis.lat(i);
coseis.lon=coseis.lon(i);
coseis.stdn=coseis.stdn(i);
coseis.stde=coseis.stde(i);
coseis.stdu=coseis.stdu(i);
Hvar=movingSD(H,Tvar); %Get trailing variance
iter=size(Hvar,2);
nsta=size(Hvar,1);
%get maximum var and horiz. displacement at each epoch
mH=zeros(size(Hvar));
mH(:,1)=H(:,1);
mHvar=zeros(size(Hvar));
mHvar(:,1)=Hvar(:,1);
for k=2:iter
%Get maximum offsets
mHvar_temp=Hvar(:,k);
mH_temp=H(:,k);
%Get maximum offsets
gt=mH_temp>mH(:,k-1); %If new epoch offset is alrger
i=find(gt); %Which stations need updates
mH(:,k)=mH(:,k-1); %these are the old maximum values
mH(i,k)=mH_temp(i); %And these are the one being updated
%Get maximum variances
gt=mHvar_temp>mHvar(:,k-1); %If new epoch variance is alrger
i=find(gt); %Which stations need updates
mHvar(:,k)=mHvar(:,k-1); %these are the old maximum values
mHvar(i,k)=mHvar_temp(i); %And these are the one being updated
end
%Find when it drops below p% of Variance but
%respect horizontal offset threshold
I=zeros(size(Hvar)); %This flags when the threshold is breached
I1=ones(size(H,1),1);
for k=2:iter
if k==250
a=0;
end
lt=Hvar(:,k)<(p*mHvar(:,k));
i=find(lt); %the ones that breach the variance threshold
gt=H(:,k)>thresh;
j=find(gt); %The ones that are over the threshold
ij=find(lt==1 & gt==1); %The ones that satisfy both
I(ij,k)=I1(ij,1); %These are the old values
%Now validate that this hasn't happened BEFORE maximum offset is
%detected
maxH=max(mH(:,k));
i=find(mH(:,k)==maxH);
if I(i,k)==0 %Maximum offset hasn't stabilized
I(:,k)=zeros(nsta,1);
end
end
%Mark first occurence of variance threshold breahcing plus T seconds
Efinal=NaN(size(Hvar));
Nfinal=Efinal;
Ufinal=Efinal;
Tfinal=Efinal;
for k=1:nsta
Itemp=I(k,:);
i=find(Itemp);
if ~isempty(i)
i1=i(1);
if i1+Tavg>iter
i2=iter;
else
i2=i1+Tavg;
end
%Get averages over selected interval
Efinal(k,i1:end)=mean(E(k,i1:i2));
Nfinal(k,i1:end)=mean(N(k,i1:i2));
Ufinal(k,i1:end)=mean(U(k,i1:i2));
Tfinal(k,i1:end)=T(k,i1:end);
end
end
%Update trigger counter matrix
for k=1:nsta
i=find(I(k,:));
if ~isempty(i)
I(k,i(1):end)=1;
end
end
%Plot data
subplot(2,2,1)
Hfinal=sqrt(Nfinal.^2+Efinal.^2);
plot(T',H')
grid on
xlabel('Seconds after origin time')
ylabel('Horizontal Disp. (m)')
hold on
plot(Tfinal',Hfinal','LineWidth',4)
stem(206,1,'LineWidth',2)
xlim([0 300])
ylim([0 0.6])
%Plot station availability
t=1:1:iter;
nthresh=10*ones(1,iter);
NS=sum(I);
subplot(2,2,[2 4])
plot(t,NS,t,nthresh,'LineWidth',2);
grid on
xlabel('Seconds after origin','FontSize',20)
ylabel('Stations detecting anomalous motion','FontSize',20)
legend('#stations','threshold')
subplot(2,2,3)
plot(T',Hvar')
grid on
xlabel('Seconds after origin time')
ylabel('Trailing Variance (m^2)')
xlim([0 300])
figure
subplot(1,3,1)
plot(T',N')
grid on
xlabel('Seconds after origin time')
ylabel('North (m)')
hold on
plot(Tfinal',Nfinal','LineWidth',4)
xlim([0 300])
subplot(1,3,2)
plot(T',E')
grid on
xlabel('Seconds after origin time')
ylabel('East (m)')
hold on
plot(Tfinal',Efinal','LineWidth',4)
xlim([0 300])
subplot(1,3,3)
plot(T',U')
grid on
xlabel('Seconds after origin time')
ylabel('Up (m)')
hold on
plot(Tfinal',Ufinal','LineWidth',4)
xlim([0 300])
%Cut data to a particular epoch
%tcut=185; %Tokachi
tcut=202; %El Mayor
i=find(I(:,tcut));
display(['Rapid solution when ' num2str(max(size(i))) ' stations detect anomalous motion.']);
cd(dir)
coseis.N=Nfinal(i,tcut)';
coseis.E=Efinal(i,tcut)';
coseis.U=Ufinal(i,tcut)';
coseis.T=Tfinal(i,tcut)';
coseis.lat=coseis.lat(i);
coseis.lon=coseis.lon(i);
coseis.stdn=coseis.stdn(i);
coseis.stde=coseis.stde(i);
coseis.stdu=coseis.stdu(i);
%Further culling of dataset
i=find(coseis.lat<33.7);
coseis.N=coseis.N(i);
coseis.E=coseis.E(i);
coseis.U=coseis.U(i);
coseis.T=coseis.T(i);
coseis.lat=coseis.lat(i);
coseis.lon=coseis.lon(i);
coseis.stdn=coseis.stdn(i);
coseis.stde=coseis.stde(i);
coseis.stdu=coseis.stdu(i);
save('rapid_coseis14.mat','coseis')