In [1]:
%pylab inline
from __future__ import print_function
import matplotlib.pylab as plt
plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = 12, 8
import warnings
warnings.filterwarnings("ignore")
Budeme potřebovat:
- stažení seznamu s eventy
- stažení staničních metadat
- stažení seismických signálů
- odstranění přístrojové odezvy
- TauP
- trigger
- xcorrel
Načtení funkcí a modulů, clienty si pro lepší přehlednost přejmenujeme:
In [2]:
from obspy.clients.fdsn import Client as Client_FDSN
from obspy.clients.syngine import Client as Client_Syngine
from obspy.clients.fdsn.mass_downloader import RectangularDomain, Restrictions, MassDownloader
from obspy import UTCDateTime, read_inventory, read
from obspy.taup import TauPyModel
from obspy.geodetics.base import calc_vincenty_inverse, kilometer2degrees
from obspy.signal.trigger import recursive_sta_lta, plot_trigger, trigger_onset
from obspy.signal.cross_correlation import xcorr_pick_correction
from obspy.imaging.cm import pqlx
import numpy as np
import os, glob, shutil
Vstupní parametry¶
In [3]:
# oblast se stanicema
lon1 = 11; lon2 = 15
lat1 = 48.5; lat2 = 51
# časové období pro eventy
t1 = UTCDateTime('2015-01-01')
t2 = UTCDateTime('2016-01-01')
# rádius pro eventy (ve stupních)
rad1 = 20; rad2 = 100
# minimální magnitudo
minmagn = 5.0
# preferovaná frekvence signálů (v Hz)
sample_rate = 100
# simulační filtr WWSSN-SP
paz_wwssn_sp = {'poles': [-3.367788+3.731514j, -3.367788-3.731514j,
-7.037168+4.545562j, -7.037168-4.545562j],
'zeros': [0j, 0j],
'gain': 70.18386,
'sensitivity': 1.}
# trigger
sta = 1 # v sec
lta = 10 # v sec
trig1 = 4
trig2 = 3
# inicializace TAUP
model = TauPyModel(model="iasp91")
Budem potřebovat pomocné funkce pro sesazení signálů a pro vykreslování:
In [4]:
def median_picked(st, key, offset):
tr_mean = st[0].copy()
for i,tr in enumerate(st):
tr1 = tr.copy().trim(tr.stats[key]+offset[0], tr.stats[key]+offset[1])
tr1_data = tr1.data / np.linalg.norm(tr1.data)
if i == 0:
signals = tr1.data / np.linalg.norm(tr1.data)
else:
signals = np.vstack((signals, tr1.data / np.linalg.norm(tr1.data)))
tr_mean.data = np.median(signals, axis=0)
# rekursivní sta-lta trigger
df = tr_mean.stats.sampling_rate
cft = recursive_sta_lta(tr_mean.data, int(sta * df), int(lta * df))
t_teor_offs = trigger_onset(cft, trig1, trig2)[0][0] / df
tr_mean.stats.t_correl = tr_mean.stats.starttime + t_teor_offs
return tr_mean
def plot_all_picked(st, key, offset, show_mean=True):
fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(111)
for i,tr in enumerate(st):
tr1 = tr.copy().trim(tr.stats[key]+offset[0], tr.stats[key]+offset[1])
tm1 = tr1.times()
tr1_data = tr1.data / np.linalg.norm(tr1.data)
ax.plot(tm1, tr1_data)
if show_mean:
tr_mean = median_picked(st, key, offset)
ax.plot(tm1, tr_mean.data, 'k-')
ax.set_xlim((0, max(tm1)))
ax.set_title("Sorted: " + key)
return
def kresli_residua(st, baz, ax=None):
_lat, _lon, _prel, _stats = [],[],[],[]
for tr in st:
if tr.stats.station in ['AAE21']: continue
_lat.append(tr.stats.latitude)
_lon.append(tr.stats.longitude)
_prel.append(tr.stats.t_correl - tr.stats.t_teor)
_stats.append(tr.stats.station)
# korekce
if tr.stats.station in ['HKWD']:
_prel[-1] -= 1.
numcols, numrows = 100, 100
xi = np.linspace(np.amin(_lon), np.amax(_lon), numcols)
yi = np.linspace(np.amin(_lat), np.amax(_lat), numrows)
xi, yi = np.meshgrid(xi, yi)
zi = griddata(_lon, _lat, _prel, xi, yi, interp='linear')
if not ax:
fig = plt.figure(figsize=(15,10))
ax = fig.add_subplot(111)
im = ax.contourf(xi, yi, zi, 20, cmap=pqlx)
ax.scatter(_lon, _lat, c=_prel, s=50, vmin=zi.min(), vmax=zi.max(), marker='v', edgecolors='k',
linewidths=1, cmap=pqlx)
offs = 0.03
for i,stat in enumerate(_stats):
ax.text(_lon[i]+offs, _lat[i]+offs, stat)
ax.set_aspect('equal')
plt.colorbar(im, ax=ax)
# arrow
x1,x2 = ax.get_xlim()
y1,y2 = ax.get_ylim()
x0 = x1 + 0.9*(x2-x1)
y0 = y1 + 0.9*(y2-y1)
d = (x2-x1) / 7
phi = (baz - 180) * np.pi/180
ax.arrow(x0, y0, d*np.sin(phi), d*np.cos(phi), width=0.02, head_width=0.1, length_includes_head=True, color='r')
plt.draw()
Seznam eventů¶
In [5]:
# referenční bod (střed oblasti)
lat0 = 0.5 * (lat1 + lat2)
lon0 = 0.5 * (lon1 + lon2)
# eventy z finálního ISC katalogu
client = Client_FDSN('ISC')
catalog = client.get_events(starttime=t1, endtime=t2,
latitude=lat0, longitude=lon0, minradius=rad1, maxradius=rad2,
minmagnitude=minmagn, orderby="time-asc", contributor="ISC")
# vykreslení
plt.rcParams['figure.figsize'] = 12, 8
catalog.plot(),
Out[5]:
Pro zjednodušení vybereme jen jeden event:
In [6]:
# for event in catalog: # normálně by byl cyklus přes všechny eventy
event = catalog.filter("time > 2015-08-06T09:00", "time < 2015-08-06T09:30")[0]
print(event.__str__().split('\n')[0])
# vezmi nejlépe určené ohnisko u daného eventu
origin = event.preferred_origin() or event.origins[0]
print(origin.__str__().split('quality')[0])
Stažení staničních signálů a metadat¶
za pomoci MassDownloader
In [7]:
# !!!(stahování chvíli trvá, takže nejdřív spustit a pak vysvětlovat)
# stahovaný úsek: start = origin_time, délka = 30 min
starttime = origin.time
endtime = starttime + 30*60
# oblast
domain = RectangularDomain(minlatitude=lat1, maxlatitude=lat2, minlongitude=lon1, maxlongitude=lon2)
# restrikce
restrict = Restrictions(
starttime=starttime, endtime=endtime,
sanitize=True, # pouze stanice s metadatama
minimum_length=0.8,
reject_channels_with_gaps=True,
minimum_interstation_distance_in_m=1000,
exclude_networks = ["Z3"], # jsou na heslo
channel_priorities=["HHZ", "BHZ", "EHZ", "SHZ"],
location_priorities=["", "00", "01"])
# kam uložit mseedy a stationxml
mseed_storage = origin.time.strftime("%Y%m%d%H%M%S")
stationxml_storage = mseed_storage
# stahování dat a metadat
mdl = MassDownloader()
mdl.download(domain, restrict, download_chunk_size_in_mb=50, threads_per_client=3,
mseed_storage=mseed_storage, stationxml_storage=stationxml_storage),
Out[7]:
Protože nám v EIDĚ zmršili metadata, zkopírujeme si naši verzi:
In [8]:
for f in glob.glob(os.path.join('DLESS_XT_CZ', '*.xml')):
shutil.copy(f, stationxml_storage)
Načtení signálů a metadat, připnutí přístrojových odezev¶
In [9]:
# staniční metadata
files = sorted(glob.glob(os.path.join(stationxml_storage, "*.xml")))
inv = read_inventory(files[0])
for f in files[1:]:
inv += read_inventory(f)
inv.plot(projection='local', color_per_network=True),
# signály
st = read(os.path.join(mseed_storage, "*.mseed"))
print(len(st))
# přiřaď odezvu a lokaci
st.attach_response(inv)
for tr in st:
coor = inv.get_coordinates(tr.id, datetime=tr.stats.starttime)
tr.stats.latitude = coor['latitude']
tr.stats.longitude = coor['longitude']
tr.stats.elevation = coor['elevation']
Předpříprava signálů¶
Výřez signálů kolem začátku teoretické P vlny
In [10]:
for tr in st:
# teoretická P vlna
[dist_m, az, baz] = calc_vincenty_inverse(origin.latitude, origin.longitude,
tr.stats.latitude, tr.stats.longitude)
dist = kilometer2degrees(dist_m/1000.)
P_arrival = model.get_travel_times(source_depth_in_km=origin.depth/1000.,
distance_in_degree=dist, phase_list=['ttp+'])[0]
tr.stats.t_teor = origin.time + P_arrival.time
tr.trim(tr.stats.t_teor-300, tr.stats.t_teor+300)
# průměrný backazimuth
[dist_m, az, baz0] = calc_vincenty_inverse(origin.latitude, origin.longitude, lat0, lon0)
dist0 = kilometer2degrees(dist_m/1000.)
Odstranění přístrojové odezvy a simulace WWSSN-SP
In [11]:
st.detrend('linear')
st.remove_response(pre_filt=[0.1, 0.5, 2, 4], water_level=1e6, output="DISP")
st.simulate(paz_remove=None, paz_simulate=paz_wwssn_sp)
for tr in st:
tr.trim(tr.stats.t_teor-60, tr.stats.t_teor+60)
st.detrend('linear')
st.filter('bandpass', freqmin=0.2, freqmax=5., corners=2, zerophase=True)
Out[11]:
Přeškálování na 100 Hz
In [12]:
for tr in st:
if tr.stats.sampling_rate != sample_rate:
if tr.stats.sampling_rate % sample_rate == 0.:
factor = int(round(tr.stats.sampling_rate / sample_rate))
tr.decimate(factor)
else:
tr.interpolate(sample_rate)
1. Sesazení signálů podle teoretického picku¶
Průměr přes všechny stopy sesazené podle teoretického picku:
In [13]:
tr_mean = median_picked(st, 't_teor', (-10,10))
plot_all_picked(st, 't_teor', (-10,10))
Trigger na průměrnou stopu, odhalí případný větší offset všech signálů od teoretického picku:
In [14]:
# rekursivní sta-lta trigger
df = tr_mean.stats.sampling_rate
cft = recursive_sta_lta(tr_mean.data, int(sta * df), int(lta * df))
t_teor_offs = trigger_onset(cft, trig1, trig2)[0][0] / df - 10
print(t_teor_offs)
plot_trigger(tr_mean, cft, trig1, trig2)
for tr in st:
tr.stats.t_teor += t_teor_offs
tr_mean.stats.t_teor = tr_mean.stats.starttime + 10 + t_teor_offs
2. Sesazení signálů podle korelace s průměrným signálem¶
In [15]:
for tr in st:
dt, coeff = xcorr_pick_correction(tr_mean.stats.t_teor, tr_mean, tr.stats.t_teor, tr, 0.5, 2., 1., plot=False)
tr.stats.t_correl = tr.stats.t_teor + dt
# nové sesazení
tr_mean = median_picked(st, 't_correl', (-10,10))
plot_all_picked(st, 't_correl', (-10,10))
3. Poslední sesazení signálů podle korelace s průměrným signálem¶
tentokrát hledáme vlnu až 5 sec daleko od teoretického picku
In [16]:
for tr in st:
dt, coeff = xcorr_pick_correction(tr_mean.stats.t_correl, tr_mean, tr.stats.t_correl, tr, 0.5, 2., 5., plot=False)
tr.stats.t_correl += dt
tr.stats.correl = coeff
# nové sesazení
tr_mean = median_picked(st, 't_correl', (-10,10))
plot_all_picked(st, 't_correl', (-10,10))
VÝPIS: stanice ...... korelace ................ P_theo ................... P_correl-P_theo
In [17]:
for tr in st:
print("{:<15s} {:3.0f}% {} {:6.2f}".format(tr.id, 100.*tr.stats.correl,
tr.stats.t_teor, tr.stats.t_correl - tr.stats.t_teor))
Vykreslíme jen signály s horší korelací:
In [18]:
fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(111)
for tr in st:
if tr.stats.correl < 0.8:
tr1 = tr.copy().trim(tr.stats.t_correl-10, tr.stats.t_correl+10)
tm1 = tr1.times()
tr1_data = tr1.data / np.linalg.norm(tr1.data)
ax.plot(tm1, tr1_data, label="{} {:4.2f}".format(tr1.id, tr1.stats.correl))
tm1 = tr_mean.times()
ax.plot(tm1, tr_mean.data, 'k-', label='MEAN')
ax.set_xlim((0, max(tm1)))
ax.legend(loc='upper left')
ax.set_title("Signals with worse correlations")
Out[18]:
Jednoduché vykreslení odchylek od teoretických časů příchodů P vln:
In [19]:
kresli_residua(st, baz0)
Syntetické seismogramy - budoucnost pikování?¶
Idea: místo korelací se signálem zprůmětovaným přes všechny stanice zkusit korelace se syntetickými signály napočtenými pro každou stanici zvlášť
Syngine = IRIS Synthetics Engine¶
http://ds.iris.edu/ds/products/syngine/
http://nbviewer.jupyter.org/gist/krischer/3e655576e4d17e6c95f2
In [20]:
client = Client_Syngine()
Syntetický signál pro KHC a dané zemětřesení:
In [21]:
st_synt = client.get_waveforms(model="ak135f_1s", network="CZ", station="KHC", units="displacement",
eventid="GCMT:C201508060922A", starttime="P-60", endtime="P+60")
print(st_synt)
st_synt.plot()
Použijeme stejné filtrování jako u reálných signálů:
In [22]:
st_synt.simulate(paz_remove=None, paz_simulate=paz_wwssn_sp)
st_synt.filter('bandpass', freqmin=0.2, freqmax=5., corners=2, zerophase=True)
st_synt.plot()
Porovnáme reálnou a syntetickou P vlnu:
In [23]:
st_KHC = st.select(station="KHC")
st_KHC += st_synt
t1 = UTCDateTime('2015-08-06T09:34:20')
t2 = t1 + 70
print('RELATIVNI NORMOVANI')
st_KHC.plot(equal_scale=False, starttime=t1, endtime=t2)
In [24]:
print('ABSOLUTNI NORMOVANI')
st_KHC.plot(equal_scale=True, starttime=t1, endtime=t2)
