由于篇幅过长，一共分为八个文档，此为第二部分，内容如下目录：
Beamforming - FK Analysis （FK分析）
Downsampling Seismograms（下采样）
下面的脚本展示了如何对地震记录进行下采样。目前支持简单的整数下采样。如果未明确禁用，则在下采样之前应先进行低通滤波以防产生混叠。为了比较，也绘制了过滤过的未下采样数据。应用的处理步骤记录在每个Trace的trace.stats.processing中。请注意处理开始的相位，因为默认情况下应用的滤波器不是零相位类型。这可以通过手动应用零相滤波器并在下下采样期间停用自动滤波来避免（no_filter = True）。
import numpy as np
import matplotlib.pyplot as plt
import obspy
# Read the seismogram
st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new")
# There is only one trace in the Stream object, let’s work on that trace...
tr = st[0]
# Decimate the 200 Hz data by a factor of 4 to 50 Hz. Note that this
# automatically includes a lowpass filtering with corner frequency 20 Hz.
# We work on a copy of the original data just to demonstrate the effects of
# downsampling.
tr_new = tr.copy()
tr_new.decimate(factor=4, strict_length=False)
# For comparison also only filter the original data (same filter options as in
# automatically applied filtering during downsampling, corner frequency
# 0.4 * new sampling rate)
tr_filt = tr.copy()
tr_filt.filter(’lowpass’, freq=0.4 * tr.stats.sampling_rate / 4.0)
# Now let’s plot the raw and filtered data...
t = np.arange(0, tr.stats.npts / tr.stats.sampling_rate, tr.stats.delta)
t_new = np.arange(0, tr_new.stats.npts / tr_new.stats.sampling_rate,tr_new.stats.delta)
plt.plot(t, tr.data, ’k’, label=’Raw’, alpha=0.3)
plt.plot(t, tr_filt.data, ’b’, label=’Lowpassed’, alpha=0.7)
plt.plot(t_new, tr_new.data, ’r’, label=’Lowpassed/Downsampled’, alpha=0.7)
plt.xlabel(’Time [s]’)
plt.xlim(82, 83.5)
plt.suptitle(tr.stats.starttime)
plt.legend()
plt.show()Merging Seismograms（合并）
下面的例子显示了如何合并然后绘制三个重叠的地震图，最长的一个被认为是正确的。 请参阅merge（）方法的文档。
import numpy as np
import matplotlib.pyplot as plt
import obspy
# Read in all files starting with dis.
st = obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.1")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.2")
# sort
st.sort([’starttime’])
# start time in plot equals 0
dt = st[0].stats.starttime.timestamp
# Go through the stream object, determine time range in julian seconds
# and plot the data with a shared x axis
ax = plt.subplot(4, 1, 1) # dummy for tying axis
for i in range(3):
plt.subplot(4, 1, i + 1, sharex=ax)
t = np.linspace(st[i].stats.starttime.timestamp - dt,st[i].stats.endtime.timestamp - dt,st[i].stats.npts)
plt.plot(t, st[i].data)
# Merge the data together and show plot in a similar way
st.merge(method=1)
plt.subplot(4, 1, 4, sharex=ax)
t = np.linspace(st[0].stats.starttime.timestamp - dt,
st[0].stats.endtime.timestamp - dt,
st[0].stats.npts)
plt.plot(t, st[0].data, ’r’)
plt.show()Beamforming - FK Analysis （FK分析）
下列代码展示了如何使用obspy进行FK分析。这些数据来自····
我们应用下列设置执行array_processing()：
- 设置慢度网络拐点-3.0：3.0s/km，步长sl_s=0.03
- 窗口长1s，步长0.05
- 数据经过1~8Hz的带通滤波，不使用预白化
- semb_thres和vel_thres设置为无穷小的数字，不得更改。
- 时间戳timestamp设置为“mlabday”他可以被我们的绘图程序直接读取
- stime和etime定为UTCDate Time格式
输出被存储在out中。
下半部分展示了如何绘制输出。我们使用array_processing()输出的out，它是包含了时间戳、相对功率、绝对功率、反量和慢度的numpy ndarray。色条代表相对功率。
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import obspy
from obspy.core.util import AttribDict
from obspy.imaging.cm import obspy_sequential
from obspy.signal.invsim import corn_freq_2_paz
from obspy.signal.array_analysis import array_processing
# Load data
st = obspy.read("https://examples.obspy.org/agfa.mseed")
# Set PAZ and coordinates for all 5 channels
st[0].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[0].stats.coordinates = AttribDict({
'latitude': 48.108589,
'elevation': 0.450000,
'longitude': 11.582967})
st[1].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[1].stats.coordinates = AttribDict({
'latitude': 48.108192,
'elevation': 0.450000,
'longitude': 11.583120})
st[2].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 250000000.0,
'gain': 1.0})
st[2].stats.coordinates = AttribDict({
'latitude': 48.108692,
'elevation': 0.450000,
'longitude': 11.583414})
st[3].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j)],
'zeros': [0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[3].stats.coordinates = AttribDict({
'latitude': 48.108456,
'elevation': 0.450000,
'longitude': 11.583049})
st[4].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j), (-2.105 + 0j)],
'zeros': [0j, 0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[4].stats.coordinates = AttribDict({
'latitude': 48.108730,
'elevation': 0.450000,
'longitude': 11.583157})
# Instrument correction to 1Hz corner frequency
paz1hz = corn_freq_2_paz(1.0, damp=0.707)
st.simulate(paz_remove='self', paz_simulate=paz1hz)
# Execute array_processing
stime = obspy.UTCDateTime("20080217110515")
etime = obspy.UTCDateTime("20080217110545")
kwargs = dict(
# slowness grid: X min, X max, Y min, Y max, Slow Step
sll_x=-3.0, slm_x=3.0, sll_y=-3.0, slm_y=3.0, sl_s=0.03,
# sliding window properties
win_len=1.0, win_frac=0.05,
# frequency properties
frqlow=1.0, frqhigh=8.0, prewhiten=0,
# restrict output
semb_thres=-1e9, vel_thres=-1e9, timestamp='mlabday',
stime=stime, etime=etime
)
out = array_processing(st, **kwargs)
# Plot
labels = ['rel.power', 'abs.power', 'baz', 'slow']
xlocator = mdates.AutoDateLocator()
fig = plt.figure()
for i, lab in enumerate(labels):
ax = fig.add_subplot(4, 1, i + 1)
ax.scatter(out[:, 0], out[:, i + 1], c=out[:, 1],alpha=0.6,
edgecolors='none', cmap=obspy_sequential)
ax.set_ylabel(lab)
ax.set_xlim(out[0, 0], out[-1, 0])
ax.set_ylim(out[:, i + 1].min(), out[:, i + 1].max())
ax.xaxis.set_major_locator(xlocator)
ax.xaxis.set_major_formatter(mdates.AutoDateFormatter(xlocator))
fig.suptitle('AGFA skyscraper blasting in Munich %s' % (
stime.strftime('%Y-%m-%d'), ))
fig.autofmt_xdate()
fig.subplots_adjust(left=0.15, top=0.95, right=0.95, bottom=0.2, hspace=0)
plt.show()
也可以使用极坐标图表示，其将网格箱中的相对功率相加，每个网格箱都由背调角和要分析的部分信号的慢度定义。反向量从北向开始顺时针计数，慢度限制可以手动设定。
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colorbar import ColorbarBase
from matplotlib.colors import Normalize
import obspy
from obspy.core.util import AttribDict
from obspy.imaging.cm import obspy_sequential
from obspy.signal.invsim import corn_freq_2_paz
from obspy.signal.array_analysis import array_processing
# Load data
st = obspy.read("https://examples.obspy.org/agfa.mseed")
# Set PAZ and coordinates for all 5 channels
st[0].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[0].stats.coordinates = AttribDict({
'latitude': 48.108589,
'elevation': 0.450000,
'longitude': 11.582967})
st[1].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 205479446.68601453,
'gain': 1.0})
st[1].stats.coordinates = AttribDict({
'latitude': 48.108192,
'elevation': 0.450000,
'longitude': 11.583120})
st[2].stats.paz = AttribDict({
'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],
'zeros': [0j, 0j],
'sensitivity': 250000000.0,
'gain': 1.0})
st[2].stats.coordinates = AttribDict({
'latitude': 48.108692,
'elevation': 0.450000,
'longitude': 11.583414})
st[3].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j)],
'zeros': [0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[3].stats.coordinates = AttribDict({
'latitude': 48.108456,
'elevation': 0.450000,
'longitude': 11.583049})
st[4].stats.paz = AttribDict({
'poles': [(-4.39823 + 4.48709j), (-4.39823 - 4.48709j), (-2.105 + 0j)],
'zeros': [0j, 0j, 0j],
'sensitivity': 222222228.10910088,
'gain': 1.0})
st[4].stats.coordinates = AttribDict({
'latitude': 48.108730,
'elevation': 0.450000,
'longitude': 11.583157})
# Instrument correction to 1Hz corner frequency
paz1hz = corn_freq_2_paz(1.0, damp=0.707)
st.simulate(paz_remove='self', paz_simulate=paz1hz)
# Execute array_processing
kwargs = dict(
# slowness grid: X min, X max, Y min, Y max, Slow Step
sll_x=-3.0, slm_x=3.0, sll_y=-3.0, slm_y=3.0, sl_s=0.03,
# sliding window properties
win_len=1.0, win_frac=0.05,
# frequency properties
frqlow=1.0, frqhigh=8.0, prewhiten=0,
# restrict output
semb_thres=-1e9, vel_thres=-1e9,
stime=obspy.UTCDateTime("20080217110515"),
etime=obspy.UTCDateTime("20080217110545")
)
out = array_processing(st, **kwargs)
# Plot
cmap = obspy_sequential
# make output human readable, adjust backazimuth to values between 0 and 360
t, rel_power, abs_power, baz, slow = out.T
baz[baz < 0.0] += 360
# choose number of fractions in plot (desirably 360 degree/N is an integer!)
N = 36
N2 = 30
abins = np.arange(N + 1) * 360. / N
sbins = np.linspace(0, 3, N2 + 1)
# sum rel power in bins given by abins and sbins
hist, baz_edges, sl_edges = \
np.histogram2d(baz, slow, bins=[abins, sbins], weights=rel_power)
# transform to radian
baz_edges = np.radians(baz_edges)
# add polar and colorbar axes
fig = plt.figure(figsize=(8, 8))
cax = fig.add_axes([0.85, 0.2, 0.05, 0.5])
ax = fig.add_axes([0.10, 0.1, 0.70, 0.7], polar=True)
ax.set_theta_direction(-1)
ax.set_theta_zero_location("N")
dh = abs(sl_edges[1] - sl_edges[0])
dw = abs(baz_edges[1] - baz_edges[0])
# circle through backazimuth
for i, row in enumerate(hist):
bars = ax.bar(left=(i * dw) * np.ones(N2),
height=dh * np.ones(N2),
width=dw, bottom=dh * np.arange(N2),
color=cmap(row / hist.max()))
ax.set_xticks(np.linspace(0, 2 * np.pi, 4, endpoint=False))
ax.set_xticklabels(['N', 'E', 'S', 'W'])
# set slowness limits
ax.set_ylim(0, 3)
[i.set_color('grey') for i in ax.get_yticklabels()]
ColorbarBase(cax, cmap=cmap,
norm=Normalize(vmin=hist.min(), vmax=hist.max()))
plt.show()
Seismogram Envelopes（信号包络）
以下脚本显示如何过滤地震数据并将其与其包络信号一起绘制图像。本示例使用零相移带通滤波器（1~3Hz）进行滤波。然后我们计算出包络并将其与Trace一起绘制。
import numpy as np
import matplotlib.pyplot as plt
import obspy
import obspy.signal
st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new")
data = st[0].data
npts = st[0].stats.npts
samprate = st[0].stats.sampling_rate
# Filtering the Stream object
st_filt = st.copy()
st_filt.filter('bandpass', freqmin=1, freqmax=3, corners=2, zerophase=True)
# Envelope of filtered data
data_envelope = obspy.signal.filter.envelope(st_filt[0].data)
# The plotting, plain matplotlib
t = np.arange(0, npts / samprate, 1 / samprate)
plt.plot(t, st_filt[0].data, 'k')
plt.plot(t, data_envelope, 'k:')
plt.title(st[0].stats.starttime)
plt.ylabel('Filtered Data w/ Envelope')
plt.xlabel('Time [s]')
plt.xlim(80, 90)
plt.show()
Plotting Spectrograms（绘制频谱图）
下列代码展示了如何使用obspy Stream对象绘制频谱图。其中的很多选项都可以自定义，参考spectrogram（）获取更多细节。比如，可以很轻松的通过调用matplotlib.crm中预定义的色彩图用来调整所需，这里提供几个教程：
- http://www.physics.ox.ac.uk/users/msshin/science/code/matplotlib_cm/
- http://matplotlib.org/examples/color/colormaps_reference.html
- https://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
- import obspy st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new") st.spectrogram(log=True, title='BW.RJOB ' + str(st[0].stats.starttime))
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
