-
Notifications
You must be signed in to change notification settings - Fork 1
/
waveform.py
231 lines (184 loc) · 8.68 KB
/
waveform.py
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
import numpy as np
from numpy.lib.stride_tricks import as_strided
from os import path
import sys
from date2dec import *
from flatten import *
from pyrocko.gf import SatelliteTarget,Target,LocalEngine
from pyrocko import gf, moment_tensor as mtm, trace
from pyrocko import util, pile, model, config, trace, io, pile
class waveforms:
def __init__(self,network,reduction,wdir,event,weight=1.,phase='P',component='Z',
filter_corner=0.055,filter_order=4,filter_type='low',misfit_norm=2,taper_fade=2.0,
base=0,sig_base=0,extension='',dist='Unif'):
self.network=network
self.reduction=reduction
self.wdir=wdir
self.event = event
self.phase=phase
self.component=component
self.filter_corner=filter_corner
self.filter_order=filter_order
self.filter_type=filter_type
self.misfit_norm=misfit_norm
self.taper_fade=taper_fade
self.sigmad=1./weight
self.base=base
self.sig_base=sig_base
self.Mbase = 1
self.extension=extension
self.dist=dist
self.taper = trace.CosFader(xfade=self.taper_fade) # Cosine taper with fade in and out of 2s.
self.bw_filter = trace.ButterworthResponse(corner=self.filter_corner, # in Hz
order=self.filter_order,
type=self.filter_type) # "low"pass or "high"pass
self.setup = trace.MisfitSetup(description='Misfit Setup',
norm=2, # L1 or L2 norm
taper=self.taper,
filter=self.bw_filter,
domain='time_domain') # Possible domains are:
# time_domain, cc_max_norm (correlation)
# and frequency_domain
self.events = []
self.events.extend(model.load_events(filename=wdir+self.event))
origin = gf.Source(
lat=np.float(self.events[0].lat),
lon=np.float(self.events[0].lon))
# print util.time_to_str(events[0].time)
self.base_source = gf.MTSource.from_pyrocko_event(self.events[0])
self.base_source.set_origin(origin.lat, origin.lon)
# print util.time_to_str(self.base_source.time), events[0].lon, events[0].lat
# sys.exit()
self.type = 'Waveform'
def load(self,inv):
# load the data as a pyrocko pile and reform them into an array of traces
data = pile.make_pile([self.wdir+self.reduction])
self.traces = data.all()
# load station file
fname = self.wdir + self.network
stations_list = model.load_stations(fname)
for s in stations_list:
s.set_channels_by_name(*self.component.split())
self.targets = []
self.tmin, self.tmax = [], []
self.arrivals = []
self.names = []
for station,tr in zip(stations_list,self.traces): # iterate over all stations
# print station.lat, station.lon
target = Target(
lat = np.float(station.lat), # station lat.
lon = np.float(station.lon), # station lon.
store_id = inv.store, # The gf-store to be used for this target,
# we can also employ different gf-stores for different targets.
interpolation = 'multilinear', # interp. method between gf cells
quantity = 'displacement', # wanted retrieved quantity
codes = station.nsl() + ('BH'+self.component,)) # Station and network code
# Next we extract the expected arrival time for this station from the the store,
# so we can use this later to define a cut-out window for the optimization:
self.targets.append(target)
self.names.append(station.nsl()[1])
# print len(self.traces), len(self.targets)
for station,tr,target in zip(stations_list,self.traces,self.targets):
engine = LocalEngine(store_superdirs=inv.store_path)
store = engine.get_store(inv.store)
# trace.snuffle(tr, events=self.events)
arrival = store.t(self.phase, self.base_source, target) # expected P-wave arrival
# print arrival
tmin = self.base_source.time+arrival-15 # start 15s before theor. arrival
tmax = self.base_source.time+arrival+15 # end 15s after theor. arrival
# # print self.tmin,self.tmax
tr.chop(tmin=tmin, tmax=tmax)
self.tmin.append(tmin)
self.tmax.append(tmax)
self.arrivals.append(self.base_source.time+arrival)
self.Npoints = len(self.targets)
# data vector
self.d = []
self.d.append(map((lambda x: getattr(x,'ydata')),self.traces))
self.d=flatten(self.d)
# time vector
t = []
for i in xrange(self.Npoints):
t.append(self.traces[i].get_xdata())
# self.t.append(map((lambda x: getattr(x,'get_xdata()')),self.traces))
# convert time
self.t = time2dec(map(util.time_to_str, flatten(t)))
# print self.t
self.N = len(self.d)
# print self.t
# sys.exit()
def printbase(self):
print 'Time shift from GCMT:', self.base
return
def info(self):
print
print 'Waveforms data:',self.network
print 'Number of stations:', self.Npoints
print 'Lenght data vector:', self.N
def g(self,inv,m):
m = np.asarray(m)
# forward vector
self.gm=np.zeros((self.N))
# create synth traces for plot and misfit
self.syn = self.traces
# set trace to zero
for tr in self.syn:
tr.ydata.fill(0.)
start = 0
# one seismic trace can be associated to several kernel (ie coseismic+posteseismic)
for k in xrange(len(inv.kernels)):
kernel = inv.kernels[k]
# not sure: slow slips cannot produce seismic waves? so not necessary too loop over it
# if kernel.seismo is True:
# one seimic trace can be the result of slip of several patches
for j in xrange(kernel.Mseg):
seg = kernel.segments[j]
mp = as_strided(m[start:start+seg.Mpatch])
# update patch parameter
seg.ss,seg.ds,seg.x1,seg.x2,seg.x3,seg.l,seg.w,seg.strike,seg.dip = mp
# construct connecivities
if seg.connectivity is not False:
seg.connect(inv.segments[seg.connectindex])
seg.m = seg.tolist()
# update time event
seg.time += self.base
synt_traces = seg.engine(self.targets, inv.ref).\
pyrocko_traces()
# for syn,tr in zip(synt_traces,self.traces):
# print len(syn.ydata), len(tr.ydata)
# sys.exit()
temp = 0
for i in xrange(self.Npoints):
# print temp
syn = synt_traces[i]
# print len(syn.ydata)
# chop synt trace
try:
syn.chop(tmin=self.tmin[i], tmax=self.tmax[i])
gt = as_strided(self.gm[temp:temp+len(syn.ydata)])
time = as_strided(self.t[temp:temp+len(syn.ydata)])
gt += syn.ydata*inv.kernels[k].g(time)
# print len(gt), gt
# print len(self.gm), self.gm
# update synt trace: not sure that all synt. trace have same
# reference time...
# print util.time_to_str(syn.tmin), util.time_to_str(self.traces[i].tmin)
# ref time synt trace = time patch
self.syn[i].ydata += syn.ydata*inv.kernels[k].g(time)
except:
pass
temp += len(syn.ydata)
start += seg.Mpatch
return self.gm
def residual(self,inv,m):
misfit_list = [] # init a list for a all the singular misfits
norm_list = [] # init a list for a all the singular normalizations
for tr, syn in zip(self.traces, self.syn):
misfit, norm = tr.misfit(candidate=syn, setup=self.setup)
misfit_list.append(misfit), norm_list.append(norm) # append the misfit into a list
# self.res = np.asarray(misfit_list)/(self.sigmad*np.asarray(norm_list))
# print self.res
g=np.asarray(self.g(inv,m))
self.res = (self.d-g)/self.sigmad
# print self.res
return self.res