python3 操作obs obspy教程

转载

IT剑客之家 2023-08-02 13:39:17

文章标签 python3 操作obs obspy 地震 python git 文章分类 Python 后端开发

由于篇幅过长，一共分为八个文档，此为第二部分，内容如下目录：

Downsampling Seismograms（下采样）

Merging Seismograms（合并）

Beamforming - FK Analysis （FK分析）

Plotting Spectrograms（绘制频谱图）

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

python3 操作obs obspy教程_git

也可以使用极坐标图表示，其将网格箱中的相对功率相加，每个网格箱都由背调角和要分析的部分信号的慢度定义。反向量从北向开始顺时针计数，慢度限制可以手动设定。

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

python3 操作obs obspy教程_python_02

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

python3 操作obs obspy教程_git_03

Plotting Spectrograms（绘制频谱图）

下列代码展示了如何使用obspy Stream对象绘制频谱图。其中的很多选项都可以自定义，参考spectrogram（）获取更多细节。比如，可以很轻松的通过调用matplotlib.crm中预定义的色彩图用来调整所需，

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

python3 操作obs obspy教程_python_04

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。

上一篇：android 虚线圆虚线工具

下一篇：docker持续交付持续部署实现 docker 持续集成

提问和评论都可以，用心的回复会被更多人看到评论

发布评论

相关文章

官方博客	全部文章	热门标签	班级博客
了解我们	网站地图	意见反馈

鸿蒙开发者社区	51CTO学堂
51CTO	软考资讯