分析地震目录#

  • 本节贡献者: 田冬冬(作者)、姚家园(审稿)

  • 最近更新日期: 2022-07-31

  • 预计花费时间: 15 分钟


在获取地震目录之后,通常还需要对地震目录做一些简单的绘图和分析。这一节演示如何 绘制最简单的地震深度直方图、地震震级直方图以及震级-频次关系图。

首先,需要把地震目录准备好。这里我们选择下载 2000–2010 年间震级大于 4.8 级, 震源深度大于 50 km 的地震。这一地震目录中共计约 7000 个地震:

from obspy.clients.fdsn import Client

client = Client("USGS")
cat = client.get_events(
    starttime="2000-01-01",
    endtime="2010-01-01",
    minmagnitude=4.8,
    mindepth=50,
)
print(cat)
---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In [1], line 4
      1 from obspy.clients.fdsn import Client
      3 client = Client("USGS")
----> 4 cat = client.get_events(
      5     starttime="2000-01-01",
      6     endtime="2010-01-01",
      7     minmagnitude=4.8,
      8     mindepth=50,
      9 )
     10 print(cat)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/clients/fdsn/client.py:554, in Client.get_events(self, starttime, endtime, minlatitude, maxlatitude, minlongitude, maxlongitude, latitude, longitude, minradius, maxradius, mindepth, maxdepth, minmagnitude, maxmagnitude, magnitudetype, eventtype, includeallorigins, includeallmagnitudes, includearrivals, eventid, limit, offset, orderby, catalog, contributor, updatedafter, filename, **kwargs)
    552     data_stream.close()
    553 else:
--> 554     cat = obspy.read_events(data_stream, format="quakeml")
    555     data_stream.close()
    556     return cat

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/decorator.py:232, in decorate.<locals>.fun(*args, **kw)
    230 if not kwsyntax:
    231     args, kw = fix(args, kw, sig)
--> 232 return caller(func, *(extras + args), **kw)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/core/util/decorator.py:291, in map_example_filename.<locals>._map_example_filename(func, *args, **kwargs)
    289                 except IOError:
    290                     pass
--> 291 return func(*args, **kwargs)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/core/event/catalog.py:809, in read_events(pathname_or_url, format, **kwargs)
    807     return _create_example_catalog()
    808 else:
--> 809     return _generic_reader(pathname_or_url, _read, format=format, **kwargs)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/core/util/base.py:623, in _generic_reader(pathname_or_url, callback_func, **kwargs)
    619 if not isinstance(pathname_or_url, str):
    620     # not a string - we assume a file-like object
    621     try:
    622         # first try reading directly
--> 623         generic = callback_func(pathname_or_url, **kwargs)
    624     except TypeError:
    625         # if this fails, create a temporary file which is read directly
    626         # from the file system
    627         pathname_or_url.seek(0)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/decorator.py:232, in decorate.<locals>.fun(*args, **kw)
    230 if not kwsyntax:
    231     args, kw = fix(args, kw, sig)
--> 232 return caller(func, *(extras + args), **kw)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/core/util/decorator.py:142, in uncompress_file(func, filename, *args, **kwargs)
    140     return func(filename, *args, **kwargs)
    141 if not isinstance(filename, str):
--> 142     return func(filename, *args, **kwargs)
    143 elif not Path(filename).exists():
    144     msg = "File not found '%s'" % (filename)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/core/event/catalog.py:817, in _read(filename, format, **kwargs)
    812 @uncompress_file
    813 def _read(filename, format=None, **kwargs):
    814     """
    815     Reads a single event file into a ObsPy Catalog object.
    816     """
--> 817     catalog, format = _read_from_plugin('event', filename, format=format,
    818                                         **kwargs)
    819     for event in catalog:
    820         event._format = format

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/core/util/base.py:422, in _read_from_plugin(plugin_type, filename, format, **kwargs)
    420     raise TypeError(msg % (format_ep.name, ', '.join(eps)))
    421 # read
--> 422 list_obj = read_format(filename, **kwargs)
    423 return list_obj, format_ep.name

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/io/quakeml/core.py:1833, in _read_quakeml(filename)
   1810 def _read_quakeml(filename):
   1811     """
   1812     Reads a QuakeML file and returns an ObsPy Catalog object.
   1813 
   (...)
   1831     2006-09-10T04:26:33.610000Z |  +9.614, +121.961 | 9.8  MS
   1832     """
-> 1833     return Unpickler().load(filename)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/io/quakeml/core.py:149, in Unpickler.load(self, file)
    140 """
    141 Reads QuakeML file into ObsPy catalog object.
    142 
   (...)
    146 :returns: ObsPy Catalog object.
    147 """
    148 self.xml_doc = _xml_doc_from_anything(file)
--> 149 return self._deserialize()

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/io/quakeml/core.py:949, in Unpickler._deserialize(self)
    947     warnings.warn(msg, UserWarning)
    948     continue
--> 949 self._set_enum('typeCertainty', event_el,
    950                event, 'event_type_certainty')
    951 event.creation_info = self._creation_info(event_el)
    952 event.event_descriptions = self._event_description(event_el)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/io/quakeml/core.py:187, in Unpickler._set_enum(self, xpath, element, obj, key)
    185 if not isinstance(obj_type, Enum):  # pragma: no cover
    186     raise ValueError
--> 187 value = self._xpath2obj(xpath, element)
    188 try:
    189     setattr(obj, key, value)

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/io/quakeml/core.py:164, in Unpickler._xpath2obj(self, xpath, element, convert_to, namespace)
    163 def _xpath2obj(self, xpath, element=None, convert_to=str, namespace=None):
--> 164     q = self._xpath(xpath, element=element, namespace=namespace)
    165     if not q:
    166         return None

File /usr/share/miniconda3/envs/seismo-learn/lib/python3.10/site-packages/obspy/io/quakeml/core.py:210, in Unpickler._xpath(self, xpath, element, namespace)
    207     xpath = "b:%s" % xpath
    208     namespaces = {"b": self.nsmap[None]}
--> 210 return element.xpath(xpath, namespaces=namespaces)

KeyboardInterrupt: 

为了进一步对数据做处理以及绘图,我们还需要导入 NumPyMatplotlib 模块:

import matplotlib.pyplot as plt
import numpy as np

地震深度分布直方图#

为了绘制地震深度直方图,我们先从地震目录中提取出地震深度信息,并保存到数组 depth 中。 需要注意的是,ObsPy 的地震目录中地震深度的单位为 m,所以需要除以 1000.0 转换为 m。

depth = np.array([event.origins[0].depth/1000 for event in cat])

我们可以使用 Matplotlib 的 hist() 函数绘制直方图: 这里我们设置的直方的最小值为 50,最大值 700,间隔为 10,同时设置了 Y 轴以对数方式显示。 从图中可以明显看到地震深度随着深度的变化:

fig, ax = plt.subplots(1, 1)
ax.hist(depth, bins=np.arange(50, 700, 10), log=True)
ax.set_xlabel("Depth (km)")
ax.set_ylabel("Counts")
plt.show()

地震震级直方图#

同理,从地震目录中提取出地震震级信息,保存到数组 mag 中并绘图。从图中可以明显看到, 震级越小地震数目越多:

mag = np.array([event.magnitudes[0].mag for event in cat])

fig, ax = plt.subplots(1, 1)
ax.hist(mag, bins=np.arange(4.5, 8, 0.1))
ax.set_xlabel("Magnitude")
ax.set_ylabel("Counts")
plt.show()

地震震级-频度关系#

在地震学中,有一个著名的定律,叫 Gutenberg–Richter 定律(简称 GR law),该定律描述了 震级与某一地区大于等于该震级的地震数量之间的关系。该定律的表达式是:\(\log_{10} N = a - b M\)

其中,\(M\) 表示震级,\(N\) 表示震级大于等于 \(M\) 的地震数量,\(a\)\(b\) 是常数。

为了绘制地震震级-频度关系图,首先需要计算得到公式里的 \(N\),即大于等于某个特定震级 \(M\) 的地震数目。 这里,我们选择 \(M\) 的取值范围为 4.0 到 8.0,间隔为 0.1。计算 \(N\) 的方法有很多,下面的 方法使用了 Python 的列表表达式以及 numpy.sum() 函数来实现:

mw = np.arange(4.0, 8.0, 0.1)
counts = np.array([(mag >= m).sum() for m in mw])

绘图脚本如下,注意图中 Y 轴是对数坐标。从图中可以明显看到 \(\log_{10}N\)\(M\) 在 4.8 - 7.6 级 之间存在线性关系。由于我们使用的地震目录里只有 4.8 级以上的地震,所以在 4.7 级以下偏离了线性关系,而大地震由于数目太少也偏离了线性关系。

fig, ax = plt.subplots(1, 1)
ax.semilogy(mw, counts, "o")
ax.set_xlabel("Magnitude")
ax.set_ylabel("Cumulative Number")
plt.show()

更进一步,我们可以对 4.8-7.6 级之间的数据进行线性拟合,得到 GR 定律中的系数 \(a\)\(b\)。 这里我们采用 numpy.logical_and 函数找到数组 mw 中所有满足条件的元素的索引, 并使用 NumPy 的布尔索引功能筛选出满足条件的震级 mw[idx] 和对应的 counts[idx],再 使用 numpy.polyfit() 函数拟合一元一次多项式,最后绘图:

idx = np.logical_and(mw >=4.8, mw<= 7.5)
# fitting y = p[0] * x + p[1]
p = np.polyfit(mw[idx], np.log10(counts[idx]), 1)
N_pre = 10**(mw * p[0] + p[1])

fig, ax = plt.subplots(1, 1)
ax.semilogy(mw, counts, "o")
ax.semilogy(mw, N_pre, color="red", label=f"$\log_{{10}} N={p[1]:.2f}{p[0]:.2f}M$")
ax.legend()
ax.set_xlabel("Magnitude")
ax.set_ylabel("Cumulative Number")
plt.show()