SeisProc：基于ObsPy庫的地震數(shù)據(jù)處理軟件

李曉姝 郭 飚 周 聰

1）中國西安 710054 中國地震局第二監(jiān)測中心

2）中國北京 100029 中國地震局地質(zhì)研究所地震動力學(xué)國家重點實驗室

0 引言

近年來，隨著大規(guī)模流動地震臺陣的廣泛布設(shè)以及固定地震臺網(wǎng)密度的增加，海量地震數(shù)據(jù)能否被高效、及時的處理，成為地震解釋過程中重要的技術(shù)問題之一。而在進(jìn)行地震解釋時，數(shù)據(jù)處理是必要環(huán)節(jié)。目前，處理地震數(shù)據(jù)的主流軟件有SAC、GIANT、Seismic Handler 等，多在Unix 或Linux 操作系統(tǒng)上運(yùn)行，且多由Fortran 或C 語言編寫，開發(fā)接口有限，二次開發(fā)難度較大。此外，這些軟件較少提供靈活的圖形界面，不便于操作，且不能跨平臺運(yùn)行。

基于Numpy和Scipy庫，Krischer等（2015）開發(fā)了針對地震數(shù)據(jù)處理的Python庫——ObsPy，基于Python 的龐大生態(tài)系統(tǒng)，地震學(xué)軟件開發(fā)更為便利。ObsPy 整合了全球主要地震數(shù)據(jù)中心發(fā)布的數(shù)據(jù)格式、數(shù)據(jù)接口、地震數(shù)據(jù)處理軟件庫，并可使用簡單易用的接口統(tǒng)一調(diào)用所有功能。ObsPy 支持對大部分?jǐn)?shù)據(jù)格式的讀寫，取代了文件格式轉(zhuǎn)換工具，便于地震學(xué)研究軟件的快速開發(fā)。基于Python 高效、便捷、易懂的特點，ObsPy 通俗易懂，便捷使用，是一款良好的科研輔助軟件。

目前，有諸多基于ObsPy 庫開發(fā)的軟件和處理系統(tǒng)，如ObsPyck、SeisHub、WavePicker等（https://github.com/obspy/obspy/wiki），其中，ObsPyck作為GUI應(yīng)用程序，主要用于數(shù)據(jù)檢索及數(shù)據(jù)分析，可在Linux和Mac平臺使用；SeisHub是一種基于web的數(shù)據(jù)庫，用于存檔、處理和共享地球物理數(shù)據(jù)和元數(shù)據(jù)；WavePicker是一種針對ObsPy流對象的地震波震相到時拾取器。這些程序一般專注于地震數(shù)據(jù)處理的某個過程，功能比較單一，或需要借助命令行進(jìn)行操作，且對運(yùn)行環(huán)境和系統(tǒng)也有限制。

基于目前發(fā)展較為成熟且應(yīng)用廣泛的ObsPy 地震學(xué)庫，設(shè)計并實現(xiàn)一款包含去傾、去偏、濾波、去儀器響應(yīng)、震相拾取等功能，并且簡單易用可持續(xù)開發(fā)和拓展的地震數(shù)據(jù)處理軟件seismic data processing program（SeisProc）。該軟件實現(xiàn)了基礎(chǔ)地震數(shù)據(jù)處理功能，可基于軟件API 接口和框架進(jìn)行功能拓展。

1 ObsPy 庫

在科學(xué)計算領(lǐng)域，Python 提供了完善的基礎(chǔ)代碼庫，覆蓋了網(wǎng)絡(luò)、文件、GUI、數(shù)據(jù)庫、文本等大量內(nèi)容。使用Python 開發(fā)，許多功能可直接調(diào)用內(nèi)置函數(shù)來實現(xiàn)。除了內(nèi)置庫外，Python還有大量第三方庫，如NumPy、SciPy、Pandas、Matplotlib、TensorFlow等，提供諸如數(shù)值計算、數(shù)據(jù)庫、繪圖、深度學(xué)習(xí)等功能。ObsPy 是Python 中用于訪問和處理地震波形數(shù)據(jù)和元數(shù)據(jù)的一個第三方庫（http://www.obspy.org）（Beyreuther et al，2010），整合了獲取全球主要地震數(shù)據(jù)中心數(shù)據(jù)發(fā)布方法，且可將不同數(shù)據(jù)中心的數(shù)據(jù)和不同格式的數(shù)據(jù)統(tǒng)一處理，免除了格式轉(zhuǎn)換和數(shù)據(jù)移動等步驟。此外，ObsPy 集成了地震學(xué)界所用的專有庫，提供了包含信號處理、數(shù)據(jù)分析、數(shù)據(jù)可視化在內(nèi)的多個函數(shù)。使用ObsPy 開發(fā)地震數(shù)據(jù)處理軟件，可專注于數(shù)據(jù)分析解釋。

ObsPy庫包含包（Packages）、腳本（Scripts）和數(shù)據(jù)庫/服務(wù)器客戶端（Database or Web Service Access Clients）3 個部分，其中客戶端提供遠(yuǎn)程訪問波形數(shù)據(jù)功能，腳本則為命令行程序提供支持。ObsPy 包含7 個重要的包，具體功能見表1。

表1 ObsPy 框架中提供的7 個包Table 1 7 packages in the ObsPy library

在7 個包中，核心包是其他ObsPy 模塊之間的粘合劑，主要功能是整理數(shù)據(jù)并便于在其他處理模塊之間進(jìn)行數(shù)據(jù)傳遞。在核心包中，ObsPy 提供3 種數(shù)據(jù)結(jié)構(gòu)（圖1）：波形數(shù)據(jù)結(jié)構(gòu)（Stream類）、事件數(shù)據(jù)結(jié)構(gòu)（Catalog類）和臺站數(shù)據(jù)結(jié)構(gòu)（Inventory類）。Stream類和Catalog 類具有相似架構(gòu)，均可看作數(shù)據(jù)組成單元與操作函數(shù)封裝的集合體，只不過Stream 類的數(shù)據(jù)組成單元為波形數(shù)據(jù)（trace），Catalog 類的則為事件數(shù)據(jù)（events）。波形數(shù)據(jù)和必要的描述儲存在Stream 類中，當(dāng)需要使用波形數(shù)據(jù)時可以調(diào)用該類的實例進(jìn)行參數(shù)傳遞。該類是本軟件實現(xiàn)信號處理例程、數(shù)據(jù)可視化、人機(jī)交互等核心功能的基礎(chǔ)。Catalog 類則用于批量管理事件，Catalog.Event 類中包含震源位置、震級大小、發(fā)震時刻、震相到時等變量，用于描述一個地震事件。與前2類不同，Inventory類是一種層次關(guān)系結(jié)構(gòu)，頂層是Inventory，其下為若干個地震臺網(wǎng)（networks），每個地震臺網(wǎng)均有其臺網(wǎng)代碼和描述。在每個臺網(wǎng)下又有若干地震臺站（stations），設(shè)有臺站代碼、經(jīng)緯度、海拔、起始時間等信息。之后是通道，在每個臺站下一般有3 個通道（channel），分別對應(yīng)N、E、Z 三個方向分量。每個通道除設(shè)有代碼、位置、采樣率等基礎(chǔ)信息之外，還有儀器響應(yīng)信息、方位角信息等，該類為實現(xiàn)本軟件去儀器響應(yīng)的功能提供了基礎(chǔ)。

圖1 數(shù)據(jù)結(jié)構(gòu)示意Fig.1 Major data structures

2 SeisProc 設(shè)計及功能實現(xiàn)

SeisProc 是以PyQt5 作為GUI 框架，基于MVC（Model View Controller）設(shè)計模式（圖2），實現(xiàn)具有3 層架構(gòu)的一套地震數(shù)據(jù)處理軟件。3層架構(gòu)指模型層（Model）、視圖層（View）和控制層（Controller）。在控制層，定義包括地震數(shù)據(jù)讀取、波形文件預(yù)處理等在內(nèi)的基本觸發(fā)指令，這些指令以PyQt 中特殊的信號與槽的形式連接到模型層中不同的功能模塊。這些功能模塊對指令進(jìn)行邏輯判斷和數(shù)據(jù)庫存取，對數(shù)據(jù)進(jìn)行實際處理，最終視圖層將當(dāng)前操作的指令結(jié)果進(jìn)行可視化顯示。該設(shè)計的優(yōu)點在于，各模塊間是獨(dú)立解耦的，可以在部分功能搭建維護(hù)時避免全局代碼的大規(guī)模改寫。

圖2 MVC 設(shè)計模式Fig.2 MVC architectural pattern

不同于現(xiàn)有地震數(shù)據(jù)處理軟件，SeisProc 的特點是具有人機(jī)交互的可視化界面。在使用PyQt 制作界面時通常有2 種方法：①通過代碼新建一個主窗口對象，繼承PyQt 的mainwindows 類，重寫setupUi 方法添加控件，控件布局可以通過Qt 自帶的Layout 布局添加器進(jìn)行設(shè)置；②通過QtDesigner 進(jìn)行自定義設(shè)計和繪制（霍亞飛，2017）。本程序設(shè)計采用第2 種方法，即對界面進(jìn)行自定義劃分和控件添加。

將SeisProc 的主圖形界面劃分為5 個功能區(qū)：基本功能欄、預(yù)處理功能區(qū)、圖像顯示區(qū)、列表數(shù)目錄區(qū)以及緩存條（圖3）。位于界面頂部的基本功能欄主要布置常用功能開關(guān)，包括對數(shù)據(jù)圖像放大、縮小、拖動等操作功能和實時數(shù)據(jù)顯示功能。界面左側(cè)為含有濾波、去傾、去儀器響應(yīng)的數(shù)據(jù)預(yù)處理功能區(qū)。界面中間大部分區(qū)域作為圖像顯示區(qū)，以使圖像顯示最大化。右側(cè)為列表樹目錄區(qū)。可通過parent 引用快速訪問數(shù)據(jù)管理框架的類。

圖3 界面布局設(shè)計Fig.3 Interface layout design

控制層和視圖層用于與用戶交互，核心且復(fù)雜的是模型層，用于數(shù)據(jù)處理。基于ObsPy 提供的數(shù)據(jù)框架對數(shù)據(jù)組織進(jìn)行分層管理，分別為Channel 類、Station 類、File 類和Files 類，每層僅需關(guān)心當(dāng)前抽象所得模型所需實現(xiàn)的功能即可（圖4），不僅兼顧開發(fā)效率，也提升了程序的可讀性，便于后續(xù)二次拓展。

圖4 模型層——數(shù)據(jù)管理框架Fig.4 Model—Data management framework

Channel 類是對ObsPy 中Trace 類的拓展和再封裝，主要承擔(dān)地震數(shù)據(jù)波形通道層次的基本操作，包括震相拾取、振幅單位轉(zhuǎn)換、列表樹節(jié)點創(chuàng)建等。Station 類是對地震臺站的抽象，是對大部分處理函數(shù)的底層實現(xiàn)，主要有去傾、濾波、去儀器響應(yīng)、自動拾取震相、計算理論到時等。File 類提供了數(shù)據(jù)輸出接口，震相拾取、理論到時結(jié)果以文件為單位在該類完成導(dǎo)出。Files 類是數(shù)據(jù)管理框架的最高一級，主要提供控制文件的添加和移除、根據(jù)指定標(biāo)準(zhǔn)對文件進(jìn)行排序、處理后的批量導(dǎo)出、文件顯示與隱藏等操作，本身不涉及函數(shù)處理的實現(xiàn)，用于保存打開的文件信息。

3 SeisProc 運(yùn)行測試

3.1 地震波形繪制

已有地震預(yù)處理程序需在命令行操作下實現(xiàn)數(shù)據(jù)的批量視圖化完成，如SAC。而SeisProc 可批量讀入波形數(shù)據(jù)，通過Matplotlib 庫的繪圖功能顯示出來，直觀簡便。其主界面中顯示了相應(yīng)的臺站名、震中距、理論震相到時（若存在），標(biāo)記在繪圖通道中。利用SeisProc 打開通道文件示例見圖5，可見讀入波形文件，選擇單屏顯示數(shù)量為6，即可見每個臺站記錄的地震數(shù)據(jù)均按照Z、N、E 通道順序繪圖顯示，選擇文件導(dǎo)出菜單，波形數(shù)據(jù)則以.mseed 和.ascii 的格式輸出。

圖5 繪圖操作界面Fig.5 Implementation interface for drawing operations

為提升波形繪制界面的交互性，對波形子圖設(shè)計縮放、拖動、顯示光標(biāo)的功能。顯示縮放功能設(shè)計2 種方式：①僅縮放縱坐標(biāo)，實現(xiàn)圖形振幅的放大縮小功能[圖6（a）]；②以鼠標(biāo)為中心的滾輪縮放，同時完成橫坐標(biāo)與縱坐標(biāo)的縮放[圖6（b）]。縱坐標(biāo)縮放用于調(diào)節(jié)視圖中波形的振幅比例，滾輪縮放可用于觀察信號的某一部分片段，如標(biāo)記震相到時。

圖6 子圖縮放功能的實現(xiàn)（a）振幅縮放；（b）以鼠標(biāo)為中心的縮放Fig.6 Subgraph scaling function implementation

拖動功能可以通過調(diào)節(jié)波形位置來達(dá)到更好的顯示效果，如截取數(shù)據(jù)，信號集中在通道圖像前端，此時即需要拖動波形調(diào)整位置。顯示光標(biāo)可以實現(xiàn)所有通道對齊某一時刻，該功能可用于不同儀器按震中距排序的相同類型通道，即在對第1 個通道進(jìn)行震相標(biāo)記后，可在光標(biāo)后其余按震中距排序的通道上查找震相。

3.2 地震數(shù)據(jù)預(yù)處理

地震原始波形數(shù)據(jù)在應(yīng)用前需要去除噪聲干擾、儀器零點漂移等因素的影響，以使得地震事件識別和進(jìn)一步數(shù)據(jù)處理可靠、高效。SeisProc 將去偏、去傾、濾波、去儀器響應(yīng)整合為預(yù)處理功能。在理想情況下，地震記錄應(yīng)以零線為基準(zhǔn)線上下波動。但在實際觀測中，地震波形基準(zhǔn)線并不總是零線，常見零線誤差有：常量型、傾斜型、波動型（焦煜媛等，2017）。SeisProc 主要考慮了常量型和傾斜型，利用ObsPy 提供的Trace.detrend（type=’simple’,**options）接口來實現(xiàn)。在濾波器選擇上，SeisProc 提供了低通濾波和帶通濾波2 種方式，均通過Filter（type=‘type’，**args）函數(shù)來實現(xiàn)。為了更好地顯示地震信號，選擇合適的地震儀對其他類型干擾進(jìn)行壓制，但有時也存在去除儀器響應(yīng)的需求，如獲取某個臺站絕對振幅值、在不同儀器響應(yīng)的臺站之間進(jìn)行波形對比等，ObsPy 提供了Remove_response（inventory）接口來實現(xiàn)。

3.3 震相拾取

除了基本預(yù)處理功能，SeisProc 也可以在菜單欄觸發(fā)震相自動拾取功能，調(diào)用ObsPy.Trigger 模塊下的Ar_pick（）函數(shù)，即可完成自動震相拾取，見圖7，CA.WT10 臺站圖中標(biāo)出三通道P 波、S 波自動拾取結(jié)果。

圖7 震相自動拾取結(jié)果Fig.7 Phase auto-pickup results

采用默認(rèn)參數(shù)，P 波震相拾取結(jié)果較為準(zhǔn)確，接近真實情況。由于S 波的特殊性，其震相自動拾取存在較大誤差。測試過程中發(fā)現(xiàn)，部分誤差可通過修改lta_s 和sta_s 參數(shù)小幅度減少，但不能完全消除，具體改進(jìn)方法需進(jìn)一步研究。在波形數(shù)據(jù)信噪比較低、射線路徑復(fù)雜等條件下，自動識別的可靠性大幅降低。SeisProc 同時支持人工標(biāo)記和拾取震相，并將此信息保存至地震數(shù)據(jù)頭文件。如圖8（a）所示，將要標(biāo)記到時的通道數(shù)據(jù)進(jìn)行合理的局部放大，與傳統(tǒng)的標(biāo)記準(zhǔn)則一樣，把第1 個起跳時刻作為要標(biāo)記的震相到時，右鍵可見事先添加的震相名稱，確定到時后點擊要標(biāo)注的震相，即可見到時標(biāo)記信息，點擊保存則被記錄在頭文件中見圖8（b）—（c）對比結(jié)果，可見新生成的數(shù)據(jù)頭文件中已包含標(biāo)記的震相走時（震相到時與起始時間的差）。

圖8 人工拾取震相標(biāo)記到時示意（a）人工標(biāo)記到時；（b）標(biāo)記前數(shù)據(jù)頭文件；（c）標(biāo)記后數(shù)據(jù)頭文件Fig.8 Schematic diagram of manually picking up the mark

4 拓展性測試

基于模型層—視圖層—控制層分離的軟件架構(gòu)，結(jié)合ObsPy 清晰的數(shù)據(jù)結(jié)構(gòu)框架，SeisProc 實現(xiàn)了波形繪制、濾波等基本預(yù)處理功能。為驗證該軟件當(dāng)前接口的拓展性能，進(jìn)行拓展測試。

作為地震學(xué)研究中處于基礎(chǔ)而關(guān)鍵的環(huán)節(jié)，初至震相拾取速度和精度直接影響著地震精確定位、震相識別、震源機(jī)制及破裂過程研究、地震勘探以及地震層析成像中的應(yīng)用效率和精度（王彩霞等，2013）。雖然隨著計算機(jī)技術(shù)的發(fā)展，震相拾取逐漸向半自動化和自動化過渡，但有時仍需人工判定，易產(chǎn)生誤判，若在人工拾取過程中提供一些輔助信息，則會減少誤讀、誤判的發(fā)生。因此，在現(xiàn)有功能基礎(chǔ)上添加STA/ LTA 類型觸發(fā)器和頻譜圖來輔助人工拾取震相。前者可通過長短時窗平均值比值的變化反映信號能量的變化，當(dāng)比值高于閾值時，則認(rèn)為此時可能發(fā)生地震（孫印等，2018）。而頻譜圖可較好地描述信號和噪聲隨時間變化的規(guī)律，即可以細(xì)致地刻畫各類波的時頻結(jié)構(gòu)和性質(zhì)（蔡希玲等，2005）。

由ObsPy 中Classic_sta_lta（）的API（表2）可知，Classic_sta_lta（）函數(shù)需要輸入待處理通道數(shù)據(jù)、短時窗平均值和長時窗平均值3 個參數(shù)，最終仍返回一個數(shù)組。而Spectrogram（data,samp_rate,per_lap,wlen,log,outfile,fmt,axes,dbscale,mult,cmap,zorder,title,show,sphinx,clip）的接口參數(shù)較多，可根據(jù)查詢相應(yīng)接口參數(shù)意義按照需求進(jìn)行調(diào)試。

表2 典型STA/LTA 函數(shù)接口Table 2 Classic STA/LTA function interface

考慮到以上操作并非針對所有通道，結(jié)合數(shù)據(jù)結(jié)構(gòu)框架可知，Channel 類主要承擔(dān)地震數(shù)據(jù)波形通道層次的基本操作，只需在Channel類中建立新接口，并在合適位置調(diào)用即可。為了減少對軟件運(yùn)行速度的影響，將2 個操作暫時添加到右鍵菜單功能表，只需更改控制層中Popqmlmenu（）函數(shù)，增加Spectrogram（）和Classic_sta_lta（）的觸發(fā)指令，以及_Onspectrogram（）和_Onclassic_sta_lta（）的響應(yīng)函數(shù)（圖9）。上述2 個功能的試運(yùn)行結(jié)果見圖10。

圖9 拓展功能響應(yīng)函數(shù)添加代碼示例Fig.9 The addition of a feature response function

比值圖[圖10（a）]所示突躍點被認(rèn)為可能是一個地震出現(xiàn)的標(biāo)識，可在突躍點附近尋找震相，為人工震相拾取提供了有效參考。圖10（b）即為針對QH.DAT.BHZ 數(shù)據(jù)進(jìn)行的時頻分析結(jié)果。當(dāng)波形數(shù)據(jù)、STA/LTA、時頻分析圖三者聯(lián)合分析時，可以更易識別地震到時，且震相到時拾取范圍更明確。

圖10 STA/LTA 和時頻圖添加后界面顯示（a）添加STA/LTA；（b）添加時頻圖Fig.10 Graphic drawing after adding STA/LTA and time-frequency diagram

5 結(jié)束語

針對現(xiàn)有地震數(shù)據(jù)處理軟件易用性不強(qiáng)、功能較單一的問題，利用ObsPy 庫設(shè)計實現(xiàn)一款簡單、易用、功能可拓展的地震數(shù)據(jù)處理軟件SeisProc。該軟件基于Python 語言和PyQt 庫，包含交互式圖形用戶界面的地震數(shù)據(jù)預(yù)處理軟件包，可完成數(shù)據(jù)的下載、讀寫、預(yù)處理，還可根據(jù)需求添加新功能。

相對于其它地震處理軟件，SeisProc 有以下幾個主要特點：①簡單，所有功能統(tǒng)一到一個圖形用戶界面，且調(diào)用實現(xiàn)只需使用鼠標(biāo)點擊操作；②高效，可同時加載多個臺站信息；③可擴(kuò)展性，通過軟件接口可持續(xù)開發(fā)；④跨平臺，軟件基于Python 語言編寫，可在Unix、Linux 及Windows 等平臺運(yùn)行，無需改。

此外，通過定義事件、分發(fā)事件和響應(yīng)事件3 個步驟完成控制層設(shè)計，利用Qt 特有的信號與槽的機(jī)制連接相應(yīng)處理函數(shù)，最終在視圖層完成事件響應(yīng)，實現(xiàn)地震數(shù)據(jù)處理的基礎(chǔ)功能，例如濾波、去傾、去儀器響應(yīng)等。另外，考慮到程序的美觀及后續(xù)的可持續(xù)性開發(fā)，設(shè)計主窗口—子窗口的2 層架構(gòu)。除可實現(xiàn)基礎(chǔ)功能，該軟件具有較強(qiáng)的實用性和較大的發(fā)展空間，值得推廣。