Spektrum signálu a filtrování


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")

Populating the interactive namespace from numpy and matplotlib
In [2]:
from obspy.core import read, UTCDateTime

Načteme denní signál ze stanice KHC:

In [3]:
stZNE  = read('KHC*')
print(stZNE)

3 Trace(s) in Stream:
CZ.KHC..BHE | 2012-04-17T00:00:00.019500Z - 2012-04-17T23:59:59.969500Z | 20.0 Hz, 1728000 samples
CZ.KHC..BHN | 2012-04-17T00:00:00.019500Z - 2012-04-17T23:59:59.969500Z | 20.0 Hz, 1728000 samples
CZ.KHC..BHZ | 2012-04-17T00:00:00.019500Z - 2012-04-17T23:59:59.969500Z | 20.0 Hz, 1728000 samples

Jeden z možných druhů vykreslení je denní graf:

In [4]:
stZNE.select(channel="*Z").plot(type="dayplot")

Způsob vykreslení se dá modifikovat, pro lepší přehlednost si signál předfiltrujeme pomocí filter. Denní plot umožňuje zobrazit počátky i zemětřesení, katalog si stahuje sám, pokud je k dispozici.

In [5]:
stZ = stZNE.select(channel="*Z").copy()
stZ.filter("lowpass", freq=0.1, corners=2)
stZ.plot(type='dayplot', interval=120, one_tick_per_line=True, color='k', number_of_ticks=7,
        vertical_scaling_range=2e3, size=(1100,600), bgcolor='lightgray', show_y_UTC_label=False, 
        tick_format="%H:%M", x_labels_size=10, y_labels_size=10, events={'min_magnitude': 6.0})

Časově frekvenční reprezentace povrchových vln


Vybereme 90 minutové okno a vlnu převedeme na 2 Hz

In [6]:
st = stZNE.copy()
t1 = UTCDateTime('2012-04-17T04:00:00')
t2 = t1 + 90*60
st.trim(t1,t2).decimate(factor=10).detrend('linear')
for tr in st:
    tr.stats.channel = 'L' + tr.stats.channel[1:]
print(st)
st[2].plot()

3 Trace(s) in Stream:
CZ.KHC..LHE | 2012-04-17T04:00:00.019500Z - 2012-04-17T05:30:00.019500Z | 2.0 Hz, 10801 samples
CZ.KHC..LHN | 2012-04-17T04:00:00.019500Z - 2012-04-17T05:30:00.019500Z | 2.0 Hz, 10801 samples
CZ.KHC..LHZ | 2012-04-17T04:00:00.019500Z - 2012-04-17T05:30:00.019500Z | 2.0 Hz, 10801 samples

V ObsPy existuje několik způsobů, jak počítat časově frekvenční spektrum signálu. Místo Stream.spectrogram využijeme propracovanějšího modulu tf_misfit.

In [7]:
from obspy.signal.tf_misfit import plot_tfr
from obspy.imaging.cm import pqlx, viridis_r   # různé barevné škály
plt.rcParams['figure.figsize'] = 12, 8
In [8]:
for tr in st:
    if tr.stats.channel[-1] == 'N':
        continue
    print(tr.id)
    plot_tfr(tr.data, dt=tr.stats.delta, w0=8, fmin=0.01, fmax=0.5, fft_zero_pad_fac=1, cmap=pqlx)

CZ.KHC..LHE

CZ.KHC..LHZ

Odfiltrujeme šum na 4 s. Použit je Butterworth bandpass filtr 2. řádu tam a zpátky (neposunuje vlnu), celkem tedy 4. řádu.

In [9]:
stZ_filt = st.select(channel="*Z").copy()
stZ_filt.filter('bandpass', freqmin=0.0124, freqmax=0.1, corners=2, zerophase=True)
plot_tfr(stZ_filt[0].data, dt=stZ_filt[0].stats.delta, w0=8, fmin=0.01, fmax=0.5, fft_zero_pad_fac=1, cmap=pqlx)

Spektrum ještě jednou, tentokrát přímo pro povrchové vlny a korektní barevnou škálu:

In [10]:
t = stZ_filt[0].stats.starttime
tr = stZ_filt.slice(t+2400, t+3900)[0]
plot_tfr(tr.data, dt=tr.stats.delta, w0=8, fmin=0.02, fmax=0.1, fft_zero_pad_fac=2, cmap=viridis_r)

Příklad, jak pomocí envelope dostat a následně i vykreslit časovou obálku signálu.

In [11]:
from obspy.signal.filter import envelope

fig = plt.figure(figsize=(15,4))
ax = fig.add_subplot(111)

tm = tr.times()
data_envelope = envelope(tr.data)

ax.plot(tm, tr.data, 'k-')
ax.plot(tm, data_envelope, 'b:')
ax.plot(tm, -data_envelope, 'b:')
ax.set_xlim(tm[0],tm[-1])
Out[11]:
(0.0, 1500.0)