ObsPy：将地震学引入科学Python生态系统的桥梁
2016-08-10 10:49LionKrischerTobiasMegiesRobertBarschMoritzBeyreutherThomasLecocqCorentinCaudronJoachimWassermann
地球与行星物理论评(中英文) 2016年4期
Lion Krischer Tobias Megies Robert Barsch Moritz Beyreuther Thomas Lecocq Corentin Caudron Joachim Wassermann
ObsPy：将地震学引入科学Python生态系统的桥梁
Lion KrischerTobias MegiesRobert BarschMoritz BeyreutherThomas LecocqCorentin CaudronJoachim Wassermann
摘要NumPy和SciPy这两个Python库是十分强大的数值处理和分析工具，适用于多种应用程序。我们开发了一个Python库ObsPy(http://obspy.org)，目的是使地震学软件包和工作流程的发展更为便利，也利用这些功能为地震学进入更大的科学Python生态系统建桥铺路。许多领域的科学家希望转化他们现有的工具和程序，以便能够利用Python所提供的这类平台环境，但是经常遇到下述困扰，例如特殊的文件格式、未知的专业术语，以及找不到合适的办法来替代软件中的某一个重要功能。我们提出一种方案，即在科学的NumPy包上层实现特定领域的时间序列库。据此，我们显示了一个时间序列数据的内部抽象表现的具体化实现，它能支持各种不同文件格式的读写。随后我们仔细描述了已充分发挥作用的旧代码的集成与改造，使它们能够在Python编写的现代工作流程之中继续发挥作用。最后我们举例研究如何将科研代码整合到ObsPy中，使其受众更为广泛。虽然本文给出的例子针对的是地震学，但是其中许多概念和抽象方法都可以直接应用于其他学科，特别是那些重点放在时间序列分析上的学科。
关键词地震学时间序列分析Python地震信号处理NumPySciPy
0引言
围绕NumPy[23]和SciPy建立的科学Python生态系统为所有科学、数学和工程领域带来很大的发展可能。它使得创建多效而强大的工作流程和应用程序成为可能。多年来，很多科学领域中都发展了它们自身的一套文件格式、工具和分析软件，适于完成自身特定的任务，而试图将其扩展至更灵活的用途时，常常会太困难甚至不可能。因此，一种像SciPy库(http://scipy.org/stackspec.html)提供的更为普适的方式是受欢迎的。
ObsPy[2]对地震学界内通常使用的几乎所有文件格式提供读写支持，它取代了大量的文件格式转换工具；在这样广泛的输入/输出支持的基础上，它用地震学家之间交流的专业术语提供信号处理程序；第三个里程碑即整合了获取由世界范围的地震数据中心发布数据的方法。最后，它集成了大量地震学界所用的专有库，并用易用的接口统一了所有功能的调用。
1特定专业领域的时间序列分析
ObsPy库包括一个特定专业领域的时间序列分析工具包，它能帮助地震学家以他们熟悉的注记符号来建立处理流程。它也为那些可能还在犹豫是否要投入学习Python及其科学生态系统的业内专家提供了一个由NumPy和SciPy支持的多项功能的接口。第一部分包括一个技术描述，说明了大量不同的时间序列文件格式是如何由统一的方式进行处理的，以及怎样使用插件系统来将由NumPy和SciPy提供的通用信号处理程序映射到上述方便的接口。这种办法的可行性已经得到证明，我们相信它很快会被应用到其他领域。
1.1地震波形
地震学主要研究在固体介质中传播的弹性波。由天然源(构造地震、火山、海浪)或人工源(核爆、矿山爆破、诱发地震)产生的地震波动结果被地震仪器测量并记录。它们是对多达3个互相垂直方向的弹性波场进行单点观测。每个方向或分量都以一维时变信号的形式测量地面运动的位移、速度或加速度。产生的等间隔采样时间序列被称作地震图或地震波形。
过去的几十年中，随着数字地震学的出现，大量地震波形文件格式涌现出来。其中一些，比如用作数据存档和数据流的MiniSEED格式[18]，是考虑一些特殊目的而设计的，其他大多数格式是作为采用特定输入/输出格式的地震信号处理或分析软件包的副产品而出现的。对其中一些软件包的应用导致它们的输入/输出格式变成广泛使用的数据交换格式，尽管这并不是当初的设计目标。例如有个反复出现的问题就是字节顺序：多数格式并不指明它们是以大字节存储顺序还是小字节存储顺序写入文件，两者都采用的情况也有所见。各式各样的文件格式以及软件只能接受一种特殊文件格式的情况，使得许多格式转换程序被编写出来并发布。
每个可用的波形格式本质上都是储存了一个或多个时间序列，以及与这个序列相关的元信息。这些格式的必要组成部分可以提炼为一个公共的信息子集，从而能够使用统一的内部表示。ObsPy实现了上述想法，我们将在下文对其进行具体描述。
1.2波形数据格式的特征
假设一个地震数据记录从A时间开始，B时间结束，并且进行等间隔采样，没有断点和重复点。那么它可以由以下信息唯一地表示：第一个采样点的时间、采样率、记录仪器以及包含这个信号的一个数组。这是公共核心的信息，所有波形文件必须以某种形式包含这些信息。地震学中地震接收仪器可用包括4个短字符串的所谓SEED标识符[18]来唯一识别。全球范围内都协调努力来维持这个系统不变。这样产生的直接结果就是大多数波形格式都符合这个系统要求，从而解决了接收器的位置问题。
地震数据经常会出现断点和重复点，这意味着数据并不是均匀采样，前面的假设也就不再成立。数据断点有可能由记录台站短时间的断电或者数据通信问题而导致，而时钟漂移和校正则可能是造成数据重复的原因之一。一些数据格式和ObsPy解决这个问题的方法是允许波形数据的多段式存储，每段进行等间隔采样，这样在单个段内通常可以保持良好状态。
在ObsPy中，波形文件格式可能支持的其他附加信息是分开存储的，这将在下一部分进行详细介绍。
1.3内部数据表示
在ObsPy内部，波形数据是用一个Stream对象来表示的，它的行为就像一个可以容纳任意数量Trace对象的容器。ObsPy定义一条Trace为包括一个单一的连续的在一个时窗内等间隔采样的波形数据，以及与之相伴的必要的元信息。每个Trace对象具有一个data属性，即一个一维的NumPy数组。其他更多的信息都放到一个字典式的stats属性中。
stats对象存储了4个SEED标识符，包括network，station，location和channel，指明了记录器的物理位置和设备。同时还包括第一个和最后一个采样点、数据采样率、采样间隔以及采样点数。这一信息部分是多余的，ObsPy负责使其保持一致。例如，endtime属性是只读的，并且会随着起始时间、采样率或数组中采样点数的变化而自动调整。其他特殊文件格式可能包含未按此抽象方式来转述的额外信息，它们将储存到stats对象内部的一个容器中。见图1中描述的内部表示的示意图。
1.4通过插件系统拓宽的数据格式支持
为了支持尽可能多的波形文件格式，我们选择了一种模块化插件的办法，每个支持的数据格式都是其子模块，并且在pkg-util的插件系统帮助下向ObsPy进行注册。列表1显示了一个波形插件示例的注册。文件格式子模块必须实现有指定接口的两个函数，第三个函数为可选：
•is_format(filename)：格式识别功能，若传来的文件是子模块的格式，返回True，否则为False。
图1 显示ObsPy内部波形数据表示的图解：由多个Trace对象汇入一个Stream对象组成(原图为彩色图——译注)。每个Trace对象表示一个不间断的、等间隔采样的时间序列，其采样点都储存到一个NumPy数组中；附加的元数据都放到stats属性下。stats对象通过利用定制的条目设置程序使元信息保持自身一致。一个特定专业领域信号处理方法的集合使得科学家能够用他们所熟悉的专业用语创建工作流程
列表1摘自ObsPy的setup.py脚本文件，以MiniSEED插件为例说明波形插件入口点是如何定义的。distutils将负责在安装后对插件进行注册。其他Python模块可以对极少使用的格式注册它们自己波形格式的插件
ENTRY_POINTS={ ... "obspy.plugin.waveform":[ ... "MSEED=obspy.mseed.core", ...] ... "obspy.plugin.waveform.MSEED":[ "isFormat=obspy.mseed.core:isMSEED", "readFormat=obspy.mseed.core:readM-SEED", "writeFormat=obspy.mseed.core:writeM-SEED", ], ...}
•read_format(filename，**Kwargs)：读取文件并返回Stream对象，包括文件数据的表示。
•write_format(stream,filename，**kwargs)：将Stream对象写入给定文件。该函数为可选项，如果不提供，那么该格式的写出支持就不存在。
这种策略上的强制模块化产生一个干净的关注点分离，即每个子模块都能独立地实现和测试。此外，它还允许用户扩展ObsPy来对那些因不常使用而没有整合到ObsPy主库中的格式给予输入/输出支持。一个常见的例子就是针对数值波形求解程序的输出格式。pkg-utils的插件系统会将这些额外的格式与ObsPy的其他部分无缝地集成到一起。
1.5格式自动检测与用法
ObsPy有一个顶层函数read()，是一个读入波形数据时的单一入口点。见列表2中的用法举例。
read(filename，**kwargs) 程序调用所有注册格式的is_fileformat()函数，直到有一个返回True。根据正在探寻的格式，is_fileformat()程序解析开始的几个字节或执行更复杂的经验方法。出于执行效率的原因，格式检测程序的速度必须要快。于是，ObsPy按照人工挑选的清单进行格式检测，这样最常用的格式将首先检测，从而提高了平均执行效率。在格式得到确认后，就会调用相应格式的read_fileformat()函数，对文件进行解析并返回Stream对象。也可以直接为read()程序提供格式参数，这样可以省略格式检测程序。
通过将某一功能实现为read()程序的一部分，这种结构允许将此功能对所有格式进行共享。例如当检测到有效的HTTP通用资源地址时，文件会自动被下载，以及一系列不同的档案格式的解压。
世界各地数据中心提供的Web服务是波形数据的重要来源。ObsPy实现客户端能够与诸多数据中心[20]进行交互操作。一条波形数据的请求将以其被存入Stream对象结束，所以无论数据来源于何处，得到数据后的工作流程都是相同的。
1.6特定专业领域的便捷方法
统一的内部数据表示开启了定义数据转换方法的可能性。对此，ObsPy提供了全面的信号处理程序集合，这些地震学中经常使用的程序，非常依赖于NumPy和SciPy提供的功能。它们的实现既满足了地震学家对某个处理操作的需要，同时也保持数据自身的一致性。例如Trace.decimate()方法会应用滤波对数据进行抽采，然后修改时间序列的元数据的采样率。由于采样点数可能会很多，所以大多数操作是按照对Stream和Trace对象进行原地修改而实现的。
列表2一段交互的Python会话的截屏，用以表明ObsPy中read()函数的用法。它将会检测文件格式并调用相应的读入程序。如果read()流程检测到有效的HTTP通用资源地址(URL)，它就会在读入之前先下载该资源。如果检测到存档格式，它就先进行解压
>>>importobspy>>>st=obspy.read("filename")>>>st
多数处理方法是对单个Trace定义的，对Stream对象施加同样的方法将会对其每个子Trace对象调用相应的方法。这就允许对三分量甚至更多分量的数据进行处理。一些方法是只针对Stream对象的。参考表1选择可行的信号处理程序。
为了提供多种不同函数，ObsPy又一次使用了插件系统。列表3～5以Trace.taper()方法为例表明了这一点。所用的插件系统通过将不同模块定义的功能集成进简单的科学家用的API特定领域的应用程序接口之中，进而，它还鼓励了重复代码剔除和代码再利用，从而使得ObsPy能够提供更多种不同函数，而不用实现和测试大量实例中的细节。
2旧代码的集成
对于某些特定的任务，地震学界经常依靠已经使用十几年甚至更久的代码。这种大量公开和应用于各类问题，保证了这些代码在大多数情况下能按照期望来工作，因为大部分问题早已被发现和修复。虽然这些代码常常不便使用，但仍然希望使其保持活力并从过去开发的投入中获益。ObsPy集成了一些旧代码，使其能够使用现代工作流程和新近发展的数据处理方法，同时依靠稳定的、经过充分测试的代码完成特定功能。这部分内容将通过两个例子给予展示。
列表3摘自ObsPy的setup.py文件。安装完成后，这些入口点会进行注册并使pkg-utils可见。下面的代码片段通过注册不同的两端尖灭窗来展示这一点。在文本写作的时候，ObsPy显示了18个不同的尖灭窗，它们主要来自scipy.signal模块。用法说明见列表5
..."obspy.plugin.taper":[ "cosine=obspy.signal.invsim:cosTa-per", "barthann=scipy.signal:barthann", ...]...
列表4来自Trace()对象内部方法taper()的代码片段。它显示了表3中展示的注册过的功能如何调用，从而获得不同的两端尖灭窗，以及可选的关键变量是如何进行传递的。用法展示见列表5
...retrievefunctioncallfromentrypointsfunc=_getFunctionFromEntryPoint( "taper",type)...taper_sides=func(2*wlen,**kwargs)...
列表5列表3和列表4中注册过和调用过的函数的用法。注意实际的两端尖灭窗是如何形成的，在这里它们是在Python不同模块中定义的。在此基础上，此表还显示了一套链式方法调用，这之所以可能是因为这些方法会返回相应的对象。copy()的调用是必要的。由于大多数方法是原处处理方法，因此可能会改变对象。大多数情况下，这个节省内存的方法是需要的，它是为适应通常使用情况而有意设计的
importobspyst=obspy.read("filename")Taperwithacosinewindow.st2=st.copy().taper("cosine")TaperawithmodifiedBartlett-Hannwin-dow.st3=st.copy().taper("barthann")Themethodscanbechainedforacompactnota-tion.st4=st.copy().detrend("linear").taper("cosine")
2.1用iaspei-tau计算走时
地震走时的计算是地震学中的常见问题。当地震发生在X点，我们的问题是，什么震相在何时到达Y点。使用整个地球的全三维模拟[29]的成本对于许多应用来说都过高，用精度较低但更快的近似方法常常就足够了。如果一维球状均匀分层地球模型是适用的，那么常用的Buland和Chapman[3]方法即可非常快速地计算走时。一个通常被称作iaspei-tau[14，26]的Fortran程序包已经实现了这个方法。功能更多的Java版本[4]也有。obspy.taup子模块使用了Python标准库里的ctypes来提供iaspei-tau的Python接口，见列表6。图2绘制了一个地震波穿过地球的走时，横坐标以度为单位的地理距离。
图2 使用一维地球模型ak135的地震波走时(原图为彩色图——译注)。上图为由obspy.taup模块计算的一些地震波震相走时。下图显示了由TauP工具包[4]和obspy.taup模块分别计算的P波走时差。该差值主要是由内部坐标系不同带来的，地震学界早已了解[15]。偏差结果强烈地依赖于所计算的地震震相和震源-接收点的几何位置
列表6由obspy.taup得到震相走时的使用方法，震相由震中距25°，深度10km的一个震源所产生。震中距(delta)单位为球状地球模型假设下的度数，深度(depth)单位为km，地球模型为ak135。返回值为一个列表型变量，而非字典型变量，因为对于指定的距离和震相组合，震相名字并非唯一
>>>fromobspy.taupimportgetTravelTimes>>>getTravelTimes(delta=25.0,depth=10.0,model="ak135")[{"phase_name":"P", "take-offangle":28.375334, "time":323.91006,...},{"phase_name":"pP", "take-offangle":151.61096, "time":326.94455,...},...
2.2仪器响应：Stationxml控制evalresp
遍布全球的地震接收仪器是为了尽可能精确地测量并记录地震动。影响最终波形特性的因素有很多，其中包括物理仪器的频率响应、任何放大器的效应、模拟或数字滤波器以及数字化带来的影响。对许多研究地球的应用来说，为了尽可能地估计真实地震动，去除这些因素带来的影响是相当重要的。因为地震接收仪器和处理过程的影响而进行数据校正的第一步就是计算记录系统的频率响应。之后用地震图对该频率响应进行反褶积，从而能得到具有物理单位、且在研究频带内不受仪器效应影响的地震图。这个过程在地震学中被称为仪器校正。
地震的记录系统是用由不同阶段或元素所组成的线性链来描述的。SEED数据格式[18]可以精确地表示这个流程，这已经被地震学界广泛地接受并从1990年代初开始使用。最近，SEED的指定接替者StationXML[19] 已经被开发出来，并且已经开始得到地震学界的接受。作为XML格式，就支持工具、方便的数据发布和人工可读性方面来说，它比二进制的SEED格式更有优势。除了一些小的变化外，这两种格式储存的信息是相同的。
列表7C结构体表示由零点和极点组成的仪器物理响应。用ctypes再次创建此内容需要定义如列表8中所示的Python类
structcomplex{doublereal;doubleimag;};structpole_zeroType{intnzeros;intnpoles;doublea0;doublea0_freq;structcomplex*zeros;structcomplex*poles;};
列表8用ctypes库创建的Python类，用来表示列表7中所示的零、pole_zeroTypeC结构体。本例中的复数(complex_number)已经定义过
ImportctypesasC...classpole_zeroType(C.Structure): _fields_=[ ("nzeros",C.c_int), ("npoles",C.c_int), ("a0",C.c_double), ("a0_freq",C.c_double), ("zeros",C.POINTER(complex_num-ber)), ("poles",C.POINTER(complex_num-ber)),]
由SEED文件计算地震接收系统频率响应的标准流程包括用rdseed[12]生成RESP文件，该文件是一个文本的SEED的真子集。然后将这些RESP文件输给evalresp[11]，得到频率响应的结果。之所以要转而通过RESP文件得到结果，是因为它是evalresp所能接受的唯一一种输入格式。用存储于StationXML文件中的元数据进行仪器校正还包括另外一个步骤，即用另外一个工具将StationXML文件转换成SEED文件。
在Python中实施这一工作包括大量系统调用和不必要的输入/输出运算，这使它变得不适用于现代化大型数据工作流程。因为在计算仪器响应的时候存在很多可能会对最终结果产生很大影响的隐患，所以取代evalresp将成为主要的努力方向。而ObsPy则提供了从StationXML(通过lxml解析)到evalresp的直接通道，借助于ctypes调用evlaresp内部函数。这种方法对于那些仍在积极开发中的代码并不可行，因为内部应用程序接口的稳定性得不到保证。不过evalresp除了需要偶尔的维护，并不再进行主动的开发。
ObsPy对StationXML数据的内部表示用于导出嵌套的C结构体，符合evalresp内部函数的期望。它们是用ctypes按照列表7和8中表示的那样进行定义和初始化的。Ctypes还需要如列表9和10中所列函数头的定义。图3是由StationXML得到的振幅和相位响应，同时展示了记录系统不同阶段的效应。
地震接收仪的记录系统通常包括许多步骤，最终在evalresp内部生成一个相当复杂的表示。evalresp的文件解析步骤负责将SEED结构翻译成相应的内部结构，有时候做出一些不明显的决定。为确保evalresp集成无误，我们进行了广泛的测试[16]。我们定义了作用在StationXML文件上的集成evalresp的正确性，即如果最终的响应和由StationXML转到SEED，再由SEED转换得到RESP文件，用evalresp对RESP文件进行计算的结果相同，就证明集成是正确的。我们从最大的地震学数据发布机构美国地震学联合研究协会下载了几乎所有可用的StationXML库存数据。这些数据来自27 000多个台站，其中大多数具有被不同时间周期所定义的多个记录通道，这样就有超过100 000个仪器响应。我们反复执行程序，直到通过测试，以此证明我们的解决方案可以可靠地应用于各地地震社团中遇到的任何数据。
图3 XM.05..HHZ通道的振幅和相位响应，用来将速度(m·s-1)记录的地震图转换为数字计数(原图为彩色图——译注)。记录仪为GuralpCMG-6TD地震计。绿线表示只考虑物理仪器和模数转换增益的响应；黑线除此以外还合并了3个相继的数字抽样步骤；红线是这一通道的柰奎斯特频率；绿线和黑线结果是用evalresp的ObsPy整合工具计算的。为了进行对比，我们用橘黄色线表示由JEvalResp[13]计算的相应结果，evalresp的Java版本是很少的能够进行这个计算的软件包之一。结果一致性非常好，奈奎斯特频率后面的差别是由绘图所用离散频点的间距引起的
列表9evalrespC代码中的函数说明。这是用于计算响应的主函数。注意输入如何使用了各种不同的参数，从简单的整数到自定义结构的数组。complex的结构体在Python里并不需要定义，因为numpy.complex128的dtype具有相同的内存布局
...structcomplex{ doublereal; doubleimag;};...voidcalc_resp(structchannel*chan, double*freq, intnfreqs, structcomplex*output, char*out_units, intstart_stage, intstop_stage, intuseTotalSensitivityFlag);
列表10列表9中各项在ctypes中对应的说明。numpy.ctypeslib模块提供的数组类型会生成自动的类型、维数和函数调用的检测标志。这产生出在调用共享库之前Python一侧的便捷的调用语法和错误处理
fromobspy.signal.headersimportclibevrespimportctypesasC...clibevresp.calc_resp.argtypes=[ C.POINTER(channel), np.ctypeslib.ndpointer( dtype='float64',ndim=1, flags='C_CONTIGUOUS'), C.c_int,np.ctypeslib.ndpointer( dtype='complex128',ndim=1, flags='C_CONTIGUOUS'), C.c_char_p,C.c_int,C.c_int, C.c_int]clibevresp.calc_resp.restype=C.c_void_p
3将研究代码集成到ObsPy
3.1目标
与上述那些已经久负盛誉、几乎无漏洞的旧代码不同，许多科学家将Python用于他们的分析工作流程，倾向于开发零散代码用于他们的输入数据，因此第一眼看上去对其他人没什么用。当输入数据有不同的格式，而输出结果需要与不同下游处理软件兼容等情况存在时，这个一般性断言在地震学中也是真的。本节中，我们首先给一个案例，研究如何翻译一个能读取奇异格式Win的C代码，然后呈现一个博士论文期间开发的科研代码，该代码现在已作为一个新模块被整合进ObsPy中。由于ObsPy的插件功能，只有必要的处理步骤需要编写。
3.2读取奇异的格式——obspy.win
WIN是由Hakusan日本公司构想的一种格式，是他们Datamark数据采集器系列的缺省存储格式。一分钟的WIN文件数据是高度压缩的，每秒的数据和前一秒相比都有不同程度的压缩。截至目前，有3种方法可以将WIN格式转换为更标准的SAC格式，包括Linux下的Win2Sac和japan2sac程序，以及Windows下的有图形用户界面的Win2Sac程序。Win2Sac的代码来自日本，而japan2sac是Chambery(法国)一个科研项目的一部分，尚未全部完成。为了用最新的程序和软件来处理这些数据，例如用基于ObsPy的脚本，使用者还得将全部波形存档转换为SAC格式，从而可以用Obs-Py读取。这个过程的每一步都使硬盘上存档的体积加倍，同时也可能在最终产出结果中加入了错误。
幸好Win2Sac的源代码[22]可在线获得，同时还提供了描述WIN格式的数据表单[30]，我们已经可以写ObsPy的插件来对其进行读取。最终的目标是能够直接不必重复地从波形存档中读取数据，并且与处理其他格式的数据一样来处理这些数据。即使有的目标是将WIN转换成MiniSEED格式，但是现在最好的解决办法还是用Obs-Py写Python脚本，而不再依靠一系列格式转换步骤。
3.3将科研代码改写成ObsPy的扩展模块
一座火山的状态及其活跃性可以被其侧翼布设的地震计记录到的振幅大小来表明。实时地震振幅测量[RSAM(1)]最早是由Endo和Murray[8]引入到火山喷发预测和火山活动性评估的工作中。我们也可以计算一定时间段内信号的振幅平方标准偏差，或者其实时地震能量测量[RSEM(2)][5]，再除以观测个数：
(1)
(2)
图4 在一个火山活动比较安静的时期，PSG台在5～20Hz频带内的地震振幅(RSAM)每日和每周的涨落。星期五(做礼拜的日子)的变化显而易见。振幅单位为counts(原图为彩色图——译注)
图5 从火山活动性比较活跃的时期到比较安静的时期火山口附近(POS台和PUN台站)以及火山登山露营基地附近台站(PSG台站)记录到的过渡变化。尽管计算的是记录第十个百分点的振幅，但从PSG台站仍然能看到人类活动的影响(原图为彩色图——译注)
式中Ai是地震信号的振幅，Aavg是N个采样点的平均振幅。
当多个震源重叠时，火山地震学家通常是在事先定义好的频带内分别计算RSEM或RSAM值，即地震能量谱测量(SSEM)[28]，类似于Stephens[27]推荐的SSAM。首先对数据进行去平均和二阶巴特沃思带通滤波。然后在30s的时窗(100×30个采样点)内计算每个值。最后计算每天第10个和第25个百分点处的值，以及中位值，用以部分去除不需要的临时信号和地震事件，它们被看成是火山脉动的干扰[7]。
obspy.signal.ssxm子模块中实现了RSAM，RSEM，SSAM和SSEM算法。最终的代码利用了Pandas(Python数据分析库)的重采样功能，完成后SSXM()方法代码的总行数不到50行。
3.3.1事例结果
用不同频带对信号进行滤波的好处可以用卡瓦伊真火山(印度尼西亚东爪哇)地震数据进行描述。地震时间序列通常会被人类日间活动所污染，特别是在高频段(>10Hz)以及那些靠近“人文”噪声源的台站，比如停车场和火山登山的露营基地。我们甚至能观测到每日的工作活动在做礼拜的日子变弱(即图4中的星期五)。为了最小化瞬时信号的影响，人们可以在特定的频带内(图5)使用第10个百分点处(p10)的振幅。除了火山活跃期(图5)，与POS和PUN台站相比，PSG台站的p10值仍然被人类活动噪声严重地影响着，这两个台站位于火山口附近，远离人类活动。
其他研究还通过对地震信号进行滤波而调查到一些特殊的过程。比如如果任意尺度上的持续破裂都会产生持续的地震信号[6]，那么SSEM的计算可能会提供应变释放率的测量。人们还可以对微震噪声源(0.1～1Hz)的平稳性进行评定来估计通过由地震环境噪声互相关技术[17]得到的速度变化的结果。
4结论
旧代码和数据格式广泛应用于各个学科，在可见的未来，它们将继续发挥重要的作用。本文所提出的方法有助于将它们集成到现代化的计算机环境中。文中展示的大多数设计和实现时的选择都可以转移到计算科学的其他领域。从复杂的文件格式集合中抽取出它们的共性，这使得构建数据源独立的工作流程成为可能，尤其有利于数据繁重的领域。同时结合使用特定专业领域词汇的功能库，包括一些紧要的旧代码，为那些甚至不熟悉Python的人提供了诱人的软件包。
ObsPy在地震学界享有广泛的认同度。我们相信这主要归功于它解决了一个实际的问题，那就是大量不同文件格式之前需要过多的文件格式转换工具或不同的信号处理工具。一旦人们开始使用ObsPy，他们将会很快发现Python的灵活和强大。相较MATLAB而言，使用ObsPy还能获得全功能编程语言的优势。另一个有竞争力的特性是能够使用许多Python中第三方模块，虽然它们并不与信号处理直接相关，但能够使用数据库、Web服务、机器学习库，当然还有最近在处理大数据集方面的新进展，它们在地震学领域和其他学科内将变得越来越重要。另外，Obspy的完整包是免费开源的，并且能够在几乎所有相关平台上运行。
成功使用ObsPy的研究包括事件重定位[21]、旋转[10]和时变[24]的地震学、大数据处理[1]，以及综合研究全波形反演[25]和衰减核[9]等。ObsPy的网站可通过链接http://www.obspy.org访问，其中包括详细的文档、详尽的教程和一个邮件列表，该列表的目的是为了使转换到ObsPy的过程更为容易，以及建立关于ObsPy的社交圈。本文中所描述的模块化和测试驱动的研发有利于添加新功能。借此，ObsPy不断地增加外部贡献的个数，目标是变为由整个业界的人士共同维护的代码。
参考文献
[1]Atkinson M，Baxter R，Brezany P，Corcho O，Galea M，Parsons M，Snelling D and van Hemert J 2013TheDataBonanza：ImprovingKnowledgeDiscoveryinScience，Engineering，andBusiness(Hoboken，NJ：Wiley)
[2]Beyreuther M，Barsch R，Krischer L，Megies T，Behr Y and Wassermann J 2010 ObsPy：a python toolbox for seismologySeismol.Res.Lett.81 530
[3]Buland R and Chapman C 1983 The computation of seismic travel timesBull.Seismol.Soc.Am.73 1271-302
[4]Crotwell H P，Owens T J and Ritsema J 1999 The taup toolkit：flexible seismic travel-time and ray-path utilitiesSeismol.Res.Lett.70 154-60
[5]de la Cruz-Reyna S and Reyes-Dvila G 2001 A model to describe precursory material-failure phenomena：applications to short-term forecasting at colima volcano，MexicoBull.Volcanol.63 297-308
[6]de la Cruz-Reyna S，Trraga M，Ortiz R and Martínez-Bringas A 2010 Tectonic earthquakes triggering volcanic seismicity and eruptions.Case studies at tungurahua and popocatpet volcanoesJ.Volcanol.Geotherm.Res.193 3748
[7]di Grazia G，Falsaperla S and Langer H 2006 Volcanic tremor location during the 2004 mount etna lava effusionGeophys.Res.Lett.33 4
[8]Endo E and Murray T 1991 Real-time seismic amplitude measurement：a volcano monitoring and prediction toolBull.Volcanol.53 533-45
[9]Fichtner A and van Driel M 2014 Models and Frechet kernels for frequency-(in)dependent QGeophys.J.Int.198 1878-89
[10]Hadziioannou C，Gaebler P，Schreiber U，Wassermann J and Igel H 2012 Examining ambient noise using colocated measurements of rotatio-nal and translational motionJ.Seismol.16 787-96
[11]IRIS 2014 Iris：Dms：Nodes：Dmc：Software downloads：evalresp http：//iris.edu/dms/nodes/dmc/software/downloads/evapresp/
[12]IRIS 2014 IRIS：DMS：Nodes：DMC：Software Downloads：rdseed http：//iris.edu/dms/nodes/dmc/software/downloads/rdseed/
[13]isti 2014 JEvalResp http：//isti.com/JEvalResp/
[14]Kennett B L N and Engdahl E R 1991 Traveltimes for global earthquake location and phase identificationGeophys.J.Int.105 429-65
[15]Knapmeyer M 2005 Numerical accuracy of tra-vel-time software in comparison with analytic resultsSeismol.Res.Lett.76 74-81
[16]Krischer L 2014 Stationxml test case git repository https：//github.com/obspy/sandbox/tree/master/stationxml_test
[17]Lecocq T，Caudron C and Brenguier F 2014 MSNoise：a python package for monitoring seismic velocity changes using ambient seismic noiseSeismol.Res.Lett.85 715-26
[18]Incorporated Research Institutions for Seismo-logy(IRIS)2012SEEDReferenceManual-StandardfortheExchangeofEarthquakeDatawww.fdsn.org/seed_manual/SEEDManual_V2.4.pdf
[19]The International Federation of Digital Seismograph Networks(FDSN)2014FDSNStationXMLSchemahttp：//fdsn.org/xml/station/
[20]Megies T，Beyreuther M，Barsch R，Krischer L and Wassermann J 2011 ObsPy-What can it do for data centers and observatories?Ann.Geophys.54 47-58
[21]Megies T and Wassermann J 2014 Microseismicity observed at a non-pressure-stimulated geothermal power plantGeothermicswww.sciencedirect.com/science/article/pii/S037565051 4000030
[22]Ohmi S 2014win2sac.chttp：//1.rcep.dpri.kyoto-u.ac.jp/～ohmi/utils/src/win2sac.c
[23]Oliphant T E 2007 Python for scientific computingComput.Sci.Eng.9 10-20
[24]Richter T，Sens-Schönfelder C，Kind R and Asch G 2014 Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometryJ.Geophys.Res.119 4747-65
[25]Schiemenz A and Igel H 2013 Accelerated 3D full-waveform inversion using simultaneously encoded sources in the time domain：application to Valhall ocean-bottom cable dataGeophys.J.Int.195 1970-88
[26]Snoke J A 2009 Traveltime tables for iasp91 and ak135Seismol.Res.Lett.80 260-2
[27]Stephens C D，Chouet B A，Page R A，Lahr J C and Power J A 1994 Seismological aspects of the 1989-1990 eruptions at redoubt volcano，Alaska：the SSAM perspectiveJ.Volcanol.Geotherm.Res.62 153-82
[29]Tromp J，Komattisch D and Liu Q 2008 Spectral-element and adjoint methods in seismologyCommun.Comput.Phys.3 1-32 http：//resolver.caltech.edu/CaltechAUTHORS：TROccp08
[30]winformat 2014 Manpage of winformat http：//eoc.eri.u-tokyo.ac.jp/cgi-bin/show_man_en?winformat
译 者 简 介
韩雪君(1984-)，女，中国地震台网中心工程师，主要从事地震观测数据处理工作。E-mail：hxj@seis.ac.cn。
Lion Krischer, Tobias Megies, Robert Barsch, Moritz Beyreuther,Thomas Lecocq,Corentin Caudron, Joachim Wassermann.2015.ObsPy: a bridge for seismology into the scientific Python ecosystem.ComputationalScience&Discovery.8: 014003.doi:10.1088/1749-4699/8/1/014003
韩雪君 译.2016.ObsPy:将地震学引入科学Python生态系统的桥梁.世界地震译丛.47(4):344-357.doi:10.16738/j.cnki.issn. 1003-3238.201604006
中国地震台网中心韩雪君译；马延路校
中国地震局地球物理研究所鲁来玉复校
猜你喜欢
小学生学习指导(中年级)(2021年4期)2021-04-27
课堂内外(初中版)(2020年5期)2020-06-19
电子制作(2019年22期)2020-01-14
动漫星空(2018年11期)2018-10-26
动漫星空(2018年2期)2018-10-26
动漫星空(2018年9期)2018-10-26
动漫星空(2018年5期)2018-10-26
软件(2016年6期)2017-02-06
中学生数理化·中考版(2015年10期)2015-09-10
智能建筑电气技术(2015年1期)2015-03-01
