2022, Vol. 43
2) 中国北京 100029 中国地震局地质研究所地震动力学国家重点实验室
2) State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China
近年来，随着大规模流动地震台阵的广泛布设以及固定地震台网密度的增加，海量地震数据能否被高效、及时的处理，成为地震解释过程中重要的技术问题之一。而在进行地震解释时，数据处理是必要环节。目前，处理地震数据的主流软件有SAC、GIANT、Seismic Handler等，多在Unix或Linux操作系统上运行，且多由Fortran或C语言编写，开发接口有限，二次开发难度较大。此外，这些软件较少提供灵活的图形界面，不便于操作，且不能跨平台运行。
基于Numpy和Scipy库，Krischer等（2015）开发了针对地震数据处理的Python库——ObsPy，基于Python的庞大生态系统，地震学软件开发更为便利。ObsPy整合了全球主要地震数据中心发布的数据格式、数据接口、地震数据处理软件库，并可使用简单易用的接口统一调用所有功能。ObsPy支持对大部分数据格式的读写，取代了文件格式转换工具，便于地震学研究软件的快速开发。基于Python高效、便捷、易懂的特点，ObsPy通俗易懂，便捷使用，是一款良好的科研辅助软件。
目前，有诸多基于ObsPy库开发的软件和处理系统，如ObsPyck、SeisHub、WavePicker等（https://github.com/obspy/obspy/wiki），其中，ObsPyck作为GUI应用程序，主要用于数据检索及数据分析，可在Linux和Mac平台使用；SeisHub是一种基于web的数据库，用于存档、处理和共享地球物理数据和元数据；WavePicker是一种针对ObsPy流对象的地震波震相到时拾取器。这些程序一般专注于地震数据处理的某个过程，功能比较单一，或需要借助命令行进行操作，且对运行环境和系统也有限制。
基于目前发展较为成熟且应用广泛的ObsPy地震学库，设计并实现一款包含去倾、去偏、滤波、去仪器响应、震相拾取等功能，并且简单易用可持续开发和拓展的地震数据处理软件seismic data processing program（SeisProc）。该软件实现了基础地震数据处理功能，可基于软件API接口和框架进行功能拓展。
1 ObsPy库在科学计算领域，Python提供了完善的基础代码库，覆盖了网络、文件、GUI、数据库、文本等大量内容。使用Python开发，许多功能可直接调用内置函数来实现。除了内置库外，Python还有大量第三方库，如NumPy、SciPy、Pandas、Matplotlib、TensorFlow等，提供诸如数值计算、数据库、绘图、深度学习等功能。ObsPy是Python中用于访问和处理地震波形数据和元数据的一个第三方库（http://www.obspy.org）（Beyreuther et al，2010），整合了获取全球主要地震数据中心数据发布方法，且可将不同数据中心的数据和不同格式的数据统一处理，免除了格式转换和数据移动等步骤。此外，ObsPy集成了地震学界所用的专有库，提供了包含信号处理、数据分析、数据可视化在内的多个函数。使用ObsPy开发地震数据处理软件，可专注于数据分析解释。
ObsPy库包含包（Packages）、脚本（Scripts）和数据库/服务器客户端（Database or Web Service Access Clients）3个部分，其中客户端提供远程访问波形数据功能，脚本则为命令行程序提供支持。ObsPy包含7个重要的包，具体功能见表 1。
|
表 1 ObsPy框架中提供的7个包 Table 1 7 packages in the ObsPy library |
在7个包中，核心包是其他ObsPy模块之间的粘合剂，主要功能是整理数据并便于在其他处理模块之间进行数据传递。在核心包中，ObsPy提供3种数据结构（图 1）：波形数据结构（Stream类）、事件数据结构（Catalog类）和台站数据结构（Inventory类）。Stream类和Catalog类具有相似架构，均可看作数据组成单元与操作函数封装的集合体，只不过Stream类的数据组成单元为波形数据（trace），Catalog类的则为事件数据（events）。波形数据和必要的描述储存在Stream类中，当需要使用波形数据时可以调用该类的实例进行参数传递。该类是本软件实现信号处理例程、数据可视化、人机交互等核心功能的基础。Catalog类则用于批量管理事件，Catalog.Event类中包含震源位置、震级大小、发震时刻、震相到时等变量，用于描述一个地震事件。与前2类不同，Inventory类是一种层次关系结构，顶层是Inventory，其下为若干个地震台网（networks），每个地震台网均有其台网代码和描述。在每个台网下又有若干地震台站（stations），设有台站代码、经纬度、海拔、起始时间等信息。之后是通道，在每个台站下一般有3个通道（channel），分别对应N、E、Z三个方向分量。每个通道除设有代码、位置、采样率等基础信息之外，还有仪器响应信息、方位角信息等，该类为实现本软件去仪器响应的功能提供了基础。
|
图 1 数据结构示意 Fig.1 Major data structures |
SeisProc是以PyQt5作为GUI框架，基于MVC（Model View Controller）设计模式（图 2），实现具有3层架构的一套地震数据处理软件。3层架构指模型层（Model）、视图层（View）和控制层（Controller）。在控制层，定义包括地震数据读取、波形文件预处理等在内的基本触发指令，这些指令以PyQt中特殊的信号与槽的形式连接到模型层中不同的功能模块。这些功能模块对指令进行逻辑判断和数据库存取，对数据进行实际处理，最终视图层将当前操作的指令结果进行可视化显示。该设计的优点在于，各模块间是独立解耦的，可以在部分功能搭建维护时避免全局代码的大规模改写。
|
图 2 MVC设计模式 Fig.2 MVC architectural pattern |
不同于现有地震数据处理软件，SeisProc的特点是具有人机交互的可视化界面。在使用PyQt制作界面时通常有2种方法：①通过代码新建一个主窗口对象，继承PyQt的mainwindows类，重写setupUi方法添加控件，控件布局可以通过Qt自带的Layout布局添加器进行设置；②通过QtDesigner进行自定义设计和绘制（霍亚飞，2017）。本程序设计采用第2种方法，即对界面进行自定义划分和控件添加。
将SeisProc的主图形界面划分为5个功能区：基本功能栏、预处理功能区、图像显示区、列表数目录区以及缓存条（图 3）。位于界面顶部的基本功能栏主要布置常用功能开关，包括对数据图像放大、缩小、拖动等操作功能和实时数据显示功能。界面左侧为含有滤波、去倾、去仪器响应的数据预处理功能区。界面中间大部分区域作为图像显示区，以使图像显示最大化。右侧为列表树目录区。可通过parent引用快速访问数据管理框架的类。
|
图 3 界面布局设计 Fig.3 Interface layout design |
控制层和视图层用于与用户交互，核心且复杂的是模型层，用于数据处理。基于ObsPy提供的数据框架对数据组织进行分层管理，分别为Channel类、Station类、File类和Files类，每层仅需关心当前抽象所得模型所需实现的功能即可（图 4），不仅兼顾开发效率，也提升了程序的可读性，便于后续二次拓展。
|
图 4 模型层——数据管理框架 Fig.4 Model—Data management framework |
Channel类是对ObsPy中Trace类的拓展和再封装，主要承担地震数据波形通道层次的基本操作，包括震相拾取、振幅单位转换、列表树节点创建等。Station类是对地震台站的抽象，是对大部分处理函数的底层实现，主要有去倾、滤波、去仪器响应、自动拾取震相、计算理论到时等。File类提供了数据输出接口，震相拾取、理论到时结果以文件为单位在该类完成导出。Files类是数据管理框架的最高一级，主要提供控制文件的添加和移除、根据指定标准对文件进行排序、处理后的批量导出、文件显示与隐藏等操作，本身不涉及函数处理的实现，用于保存打开的文件信息。
3 SeisProc运行测试 3.1 地震波形绘制已有地震预处理程序需在命令行操作下实现数据的批量视图化完成，如SAC。而SeisProc可批量读入波形数据，通过Matplotlib库的绘图功能显示出来，直观简便。其主界面中显示了相应的台站名、震中距、理论震相到时（若存在），标记在绘图通道中。利用SeisProc打开通道文件示例见图 5，可见读入波形文件，选择单屏显示数量为6，即可见每个台站记录的地震数据均按照Z、N、E通道顺序绘图显示，选择文件导出菜单，波形数据则以.mseed和.ascii的格式输出。
|
图 5 绘图操作界面 Fig.5 Implementation interface for drawing operations |
为提升波形绘制界面的交互性，对波形子图设计缩放、拖动、显示光标的功能。显示缩放功能设计2种方式：①仅缩放纵坐标，实现图形振幅的放大缩小功能[图 6（a）]；②以鼠标为中心的滚轮缩放，同时完成横坐标与纵坐标的缩放[图 6（b）]。纵坐标缩放用于调节视图中波形的振幅比例，滚轮缩放可用于观察信号的某一部分片段，如标记震相到时。
|
图 6 子图缩放功能的实现 （a）振幅缩放；（b）以鼠标为中心的缩放 Fig.6 Subgraph scaling function implementation |
拖动功能可以通过调节波形位置来达到更好的显示效果，如截取数据，信号集中在通道图像前端，此时即需要拖动波形调整位置。显示光标可以实现所有通道对齐某一时刻，该功能可用于不同仪器按震中距排序的相同类型通道，即在对第1个通道进行震相标记后，可在光标后其余按震中距排序的通道上查找震相。
3.2 地震数据预处理地震原始波形数据在应用前需要去除噪声干扰、仪器零点漂移等因素的影响，以使得地震事件识别和进一步数据处理可靠、高效。SeisProc将去偏、去倾、滤波、去仪器响应整合为预处理功能。在理想情况下，地震记录应以零线为基准线上下波动。但在实际观测中，地震波形基准线并不总是零线，常见零线误差有：常量型、倾斜型、波动型（焦煜媛等，2017）。SeisProc主要考虑了常量型和倾斜型，利用ObsPy提供的Trace.detrend（type=’simple’, **options）接口来实现。在滤波器选择上，SeisProc提供了低通滤波和带通滤波2种方式，均通过Filter（type=‘type’，**args）函数来实现。为了更好地显示地震信号，选择合适的地震仪对其他类型干扰进行压制，但有时也存在去除仪器响应的需求，如获取某个台站绝对振幅值、在不同仪器响应的台站之间进行波形对比等，ObsPy提供了Remove_response（inventory）接口来实现。
3.3 震相拾取除了基本预处理功能，SeisProc也可以在菜单栏触发震相自动拾取功能，调用ObsPy.Trigger模块下的Ar_pick（）函数，即可完成自动震相拾取，见图 7，CA.WT10台站图中标出三通道P波、S波自动拾取结果。
|
图 7 震相自动拾取结果 Fig.7 Phase auto-pickup results |
采用默认参数，P波震相拾取结果较为准确，接近真实情况。由于S波的特殊性，其震相自动拾取存在较大误差。测试过程中发现，部分误差可通过修改lta_s和sta_s参数小幅度减少，但不能完全消除，具体改进方法需进一步研究。在波形数据信噪比较低、射线路径复杂等条件下，自动识别的可靠性大幅降低。SeisProc同时支持人工标记和拾取震相，并将此信息保存至地震数据头文件。如图 8（a）所示，将要标记到时的通道数据进行合理的局部放大，与传统的标记准则一样，把第1个起跳时刻作为要标记的震相到时，右键可见事先添加的震相名称，确定到时后点击要标注的震相，即可见到时标记信息，点击保存则被记录在头文件中见图 8（b）—（c）对比结果，可见新生成的数据头文件中已包含标记的震相走时（震相到时与起始时间的差）。
|
图 8 人工拾取震相标记到时示意 （a）人工标记到时；（b）标记前数据头文件；（c）标记后数据头文件 Fig.8 Schematic diagram of manually picking up the mark |
基于模型层—视图层—控制层分离的软件架构，结合ObsPy清晰的数据结构框架，SeisProc实现了波形绘制、滤波等基本预处理功能。为验证该软件当前接口的拓展性能，进行拓展测试。
作为地震学研究中处于基础而关键的环节，初至震相拾取速度和精度直接影响着地震精确定位、震相识别、震源机制及破裂过程研究、地震勘探以及地震层析成像中的应用效率和精度（王彩霞等，2013）。虽然随着计算机技术的发展，震相拾取逐渐向半自动化和自动化过渡，但有时仍需人工判定，易产生误判，若在人工拾取过程中提供一些辅助信息，则会减少误读、误判的发生。因此，在现有功能基础上添加STA / LTA类型触发器和频谱图来辅助人工拾取震相。前者可通过长短时窗平均值比值的变化反映信号能量的变化，当比值高于阈值时，则认为此时可能发生地震（孙印等，2018）。而频谱图可较好地描述信号和噪声随时间变化的规律，即可以细致地刻画各类波的时频结构和性质（蔡希玲等，2005）。
由ObsPy中Classic_sta_lta（）的API（表 2）可知，Classic_sta_lta（）函数需要输入待处理通道数据、短时窗平均值和长时窗平均值3个参数，最终仍返回一个数组。而Spectrogram（data, samp_rate, per_lap, wlen, log, outfile, fmt, axes, dbscale, mult, cmap, zorder, title, show, sphinx, clip）的接口参数较多，可根据查询相应接口参数意义按照需求进行调试。
|
|
表 2 典型STA/LTA函数接口 Table 2 Classic STA/LTA function interface |
考虑到以上操作并非针对所有通道，结合数据结构框架可知，Channel类主要承担地震数据波形通道层次的基本操作，只需在Channel类中建立新接口，并在合适位置调用即可。为了减少对软件运行速度的影响，将2个操作暂时添加到右键菜单功能表，只需更改控制层中Popqmlmenu（）函数，增加Spectrogram（）和Classic_sta_lta（）的触发指令，以及_Onspectrogram（）和_Onclassic_sta_lta（）的响应函数（图 9）。上述2个功能的试运行结果见图 10。
|
图 9 拓展功能响应函数添加代码示例 Fig.9 The addition of a feature response function |
|
图 10 STA/LTA和时频图添加后界面显示 （a）添加STA/LTA；（b）添加时频图 Fig.10 Graphic drawing after adding STA/LTA and time-frequency diagram |
比值图[图 10（a）]所示突跃点被认为可能是一个地震出现的标识，可在突跃点附近寻找震相，为人工震相拾取提供了有效参考。图 10（b）即为针对QH.DAT.BHZ数据进行的时频分析结果。当波形数据、STA/LTA、时频分析图三者联合分析时，可以更易识别地震到时，且震相到时拾取范围更明确。
5 结束语针对现有地震数据处理软件易用性不强、功能较单一的问题，利用ObsPy库设计实现一款简单、易用、功能可拓展的地震数据处理软件SeisProc。该软件基于Python语言和PyQt库，包含交互式图形用户界面的地震数据预处理软件包，可完成数据的下载、读写、预处理，还可根据需求添加新功能。
相对于其它地震处理软件，SeisProc有以下几个主要特点：①简单，所有功能统一到一个图形用户界面，且调用实现只需使用鼠标点击操作；②高效，可同时加载多个台站信息；③可扩展性，通过软件接口可持续开发；④跨平台，软件基于Python语言编写，可在Unix、Linux及Windows等平台运行，无需改。
此外，通过定义事件、分发事件和响应事件3个步骤完成控制层设计，利用Qt特有的信号与槽的机制连接相应处理函数，最终在视图层完成事件响应，实现地震数据处理的基础功能，例如滤波、去倾、去仪器响应等。另外，考虑到程序的美观及后续的可持续性开发，设计主窗口—子窗口的2层架构。除可实现基础功能，该软件具有较强的实用性和较大的发展空间，值得推广。
蔡希玲, 吕英海. 地震数据时频分析与分频处理[J]. 勘探地球物理进展, 2005, 28(4): 265-270. |
霍亚飞. QT Creator快速入门[M]. 北京: 北京航空航天大学出版社, 2017.
|
焦煜媛, 马克博. 地震波数据格式及数据预处理技术介绍[J]. 甘肃科技, 2017, 33(20): 45-47. DOI:10.3969/j.issn.1000-0952.2017.20.017 |
孙印, 潘素珍, 刘明军. 天然地震识别与震相自动拾取技术进展[J]. 中国地震, 2018, 34(4): 606-620. DOI:10.3969/j.issn.1001-4683.2018.04.002 |
王彩霞, 白超英, 王馨. 地震震相初至自动检测技术综述[J]. 地球物理学进展, 2013, 28(5): 2363-2375. |
Beyreuther M, Barsch R, Krischer L, et al. ObsPy: a python toolbox for seismology[J]. Seismological Research Letters, 2010, 81(3): 530-533. DOI:10.1785/gssrl.81.3.530 |
Krischer L, Megies T, Barsch R, et al. ObsPy: a bridge for seismology into the scientific Python ecosystem[J]. Computational Science & Discovery, 2015, 8(1): 014003. |