ObsPy：將地震學引入科學Python生態(tài)系統(tǒng)的橋梁

Lion Krischer　Tobias Megies　Robert Barsch　Moritz Beyreuther　Thomas Lecocq　Corentin Caudron　Joachim Wassermann

?

ObsPy：將地震學引入科學Python生態(tài)系統(tǒng)的橋梁

Lion KrischerTobias MegiesRobert BarschMoritz BeyreutherThomas LecocqCorentin CaudronJoachim Wassermann

摘要NumPy和SciPy這兩個Python庫是十分強大的數(shù)值處理和分析工具，適用于多種應用程序。我們開發(fā)了一個Python庫ObsPy(http://obspy.org)，目的是使地震學軟件包和工作流程的發(fā)展更為便利，也利用這些功能為地震學進入更大的科學Python生態(tài)系統(tǒng)建橋鋪路。許多領(lǐng)域的科學家希望轉(zhuǎn)化他們現(xiàn)有的工具和程序，以便能夠利用Python所提供的這類平臺環(huán)境，但是經(jīng)常遇到下述困擾，例如特殊的文件格式、未知的專業(yè)術(shù)語，以及找不到合適的辦法來替代軟件中的某一個重要功能。我們提出一種方案，即在科學的NumPy包上層實現(xiàn)特定領(lǐng)域的時間序列庫。據(jù)此，我們顯示了一個時間序列數(shù)據(jù)的內(nèi)部抽象表現(xiàn)的具體化實現(xiàn)，它能支持各種不同文件格式的讀寫。隨后我們仔細描述了已充分發(fā)揮作用的舊代碼的集成與改造，使它們能夠在Python編寫的現(xiàn)代工作流程之中繼續(xù)發(fā)揮作用。最后我們舉例研究如何將科研代碼整合到ObsPy中，使其受眾更為廣泛。雖然本文給出的例子針對的是地震學，但是其中許多概念和抽象方法都可以直接應用于其他學科，特別是那些重點放在時間序列分析上的學科。

關(guān)鍵詞地震學時間序列分析Python地震信號處理NumPySciPy

0引言

圍繞NumPy[23]和SciPy建立的科學Python生態(tài)系統(tǒng)為所有科學、數(shù)學和工程領(lǐng)域帶來很大的發(fā)展可能。它使得創(chuàng)建多效而強大的工作流程和應用程序成為可能。多年來，很多科學領(lǐng)域中都發(fā)展了它們自身的一套文件格式、工具和分析軟件，適于完成自身特定的任務，而試圖將其擴展至更靈活的用途時，常常會太困難甚至不可能。因此，一種像SciPy庫(http://scipy.org/stackspec.html)提供的更為普適的方式是受歡迎的。

ObsPy[2]對地震學界內(nèi)通常使用的幾乎所有文件格式提供讀寫支持，它取代了大量的文件格式轉(zhuǎn)換工具；在這樣廣泛的輸入/輸出支持的基礎(chǔ)上，它用地震學家之間交流的專業(yè)術(shù)語提供信號處理程序；第三個里程碑即整合了獲取由世界范圍的地震數(shù)據(jù)中心發(fā)布數(shù)據(jù)的方法。最后，它集成了大量地震學界所用的專有庫，并用易用的接口統(tǒng)一了所有功能的調(diào)用。

1特定專業(yè)領(lǐng)域的時間序列分析

ObsPy庫包括一個特定專業(yè)領(lǐng)域的時間序列分析工具包，它能幫助地震學家以他們熟悉的注記符號來建立處理流程。它也為那些可能還在猶豫是否要投入學習Python及其科學生態(tài)系統(tǒng)的業(yè)內(nèi)專家提供了一個由NumPy和SciPy支持的多項功能的接口。第一部分包括一個技術(shù)描述，說明了大量不同的時間序列文件格式是如何由統(tǒng)一的方式進行處理的，以及怎樣使用插件系統(tǒng)來將由NumPy和SciPy提供的通用信號處理程序映射到上述方便的接口。這種辦法的可行性已經(jīng)得到證明，我們相信它很快會被應用到其他領(lǐng)域。

1.1地震波形

地震學主要研究在固體介質(zhì)中傳播的彈性波。由天然源(構(gòu)造地震、火山、海浪)或人工源(核爆、礦山爆破、誘發(fā)地震)產(chǎn)生的地震波動結(jié)果被地震儀器測量并記錄。它們是對多達3個互相垂直方向的彈性波場進行單點觀測。每個方向或分量都以一維時變信號的形式測量地面運動的位移、速度或加速度。產(chǎn)生的等間隔采樣時間序列被稱作地震圖或地震波形。

過去的幾十年中，隨著數(shù)字地震學的出現(xiàn)，大量地震波形文件格式涌現(xiàn)出來。其中一些，比如用作數(shù)據(jù)存檔和數(shù)據(jù)流的MiniSEED格式[18]，是考慮一些特殊目的而設(shè)計的，其他大多數(shù)格式是作為采用特定輸入/輸出格式的地震信號處理或分析軟件包的副產(chǎn)品而出現(xiàn)的。對其中一些軟件包的應用導致它們的輸入/輸出格式變成廣泛使用的數(shù)據(jù)交換格式，盡管這并不是當初的設(shè)計目標。例如有個反復出現(xiàn)的問題就是字節(jié)順序：多數(shù)格式并不指明它們是以大字節(jié)存儲順序還是小字節(jié)存儲順序?qū)懭胛募?，兩者都采用的情況也有所見。各式各樣的文件格式以及軟件只能接受一種特殊文件格式的情況，使得許多格式轉(zhuǎn)換程序被編寫出來并發(fā)布。

每個可用的波形格式本質(zhì)上都是儲存了一個或多個時間序列，以及與這個序列相關(guān)的元信息。這些格式的必要組成部分可以提煉為一個公共的信息子集，從而能夠使用統(tǒng)一的內(nèi)部表示。ObsPy實現(xiàn)了上述想法，我們將在下文對其進行具體描述。

1.2波形數(shù)據(jù)格式的特征

假設(shè)一個地震數(shù)據(jù)記錄從A時間開始，B時間結(jié)束，并且進行等間隔采樣，沒有斷點和重復點。那么它可以由以下信息唯一地表示：第一個采樣點的時間、采樣率、記錄儀器以及包含這個信號的一個數(shù)組。這是公共核心的信息，所有波形文件必須以某種形式包含這些信息。地震學中地震接收儀器可用包括4個短字符串的所謂SEED標識符[18]來唯一識別。全球范圍內(nèi)都協(xié)調(diào)努力來維持這個系統(tǒng)不變。這樣產(chǎn)生的直接結(jié)果就是大多數(shù)波形格式都符合這個系統(tǒng)要求，從而解決了接收器的位置問題。

地震數(shù)據(jù)經(jīng)常會出現(xiàn)斷點和重復點，這意味著數(shù)據(jù)并不是均勻采樣，前面的假設(shè)也就不再成立。數(shù)據(jù)斷點有可能由記錄臺站短時間的斷電或者數(shù)據(jù)通信問題而導致，而時鐘漂移和校正則可能是造成數(shù)據(jù)重復的原因之一。一些數(shù)據(jù)格式和ObsPy解決這個問題的方法是允許波形數(shù)據(jù)的多段式存儲，每段進行等間隔采樣，這樣在單個段內(nèi)通?？梢员３至己脿顟B(tài)。

在ObsPy中，波形文件格式可能支持的其他附加信息是分開存儲的，這將在下一部分進行詳細介紹。

1.3內(nèi)部數(shù)據(jù)表示

在ObsPy內(nèi)部，波形數(shù)據(jù)是用一個Stream對象來表示的，它的行為就像一個可以容納任意數(shù)量Trace對象的容器。ObsPy定義一條Trace為包括一個單一的連續(xù)的在一個時窗內(nèi)等間隔采樣的波形數(shù)據(jù)，以及與之相伴的必要的元信息。每個Trace對象具有一個data屬性，即一個一維的NumPy數(shù)組。其他更多的信息都放到一個字典式的stats屬性中。

stats對象存儲了4個SEED標識符，包括network，station，location和channel，指明了記錄器的物理位置和設(shè)備。同時還包括第一個和最后一個采樣點、數(shù)據(jù)采樣率、采樣間隔以及采樣點數(shù)。這一信息部分是多余的，ObsPy負責使其保持一致。例如，endtime屬性是只讀的，并且會隨著起始時間、采樣率或數(shù)組中采樣點數(shù)的變化而自動調(diào)整。其他特殊文件格式可能包含未按此抽象方式來轉(zhuǎn)述的額外信息，它們將儲存到stats對象內(nèi)部的一個容器中。見圖1中描述的內(nèi)部表示的示意圖。

1.4通過插件系統(tǒng)拓寬的數(shù)據(jù)格式支持

為了支持盡可能多的波形文件格式，我們選擇了一種模塊化插件的辦法，每個支持的數(shù)據(jù)格式都是其子模塊，并且在pkg-util的插件系統(tǒng)幫助下向ObsPy進行注冊。列表1顯示了一個波形插件示例的注冊。文件格式子模塊必須實現(xiàn)有指定接口的兩個函數(shù)，第三個函數(shù)為可選：

?is_format(filename)：格式識別功能，若傳來的文件是子模塊的格式，返回True，否則為False。

圖1　顯示ObsPy內(nèi)部波形數(shù)據(jù)表示的圖解：由多個Trace對象匯入一個Stream對象組成(原圖為彩色圖——譯注)。每個Trace對象表示一個不間斷的、等間隔采樣的時間序列，其采樣點都儲存到一個NumPy數(shù)組中；附加的元數(shù)據(jù)都放到stats屬性下。stats對象通過利用定制的條目設(shè)置程序使元信息保持自身一致。一個特定專業(yè)領(lǐng)域信號處理方法的集合使得科學家能夠用他們所熟悉的專業(yè)用語創(chuàng)建工作流程

列表1摘自O(shè)bsPy的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對象，包括文件數(shù)據(jù)的表示。

?write_format(stream,filename，**kwargs)：將Stream對象寫入給定文件。該函數(shù)為可選項，如果不提供，那么該格式的寫出支持就不存在。

這種策略上的強制模塊化產(chǎn)生一個干凈的關(guān)注點分離，即每個子模塊都能獨立地實現(xiàn)和測試。此外，它還允許用戶擴展ObsPy來對那些因不常使用而沒有整合到ObsPy主庫中的格式給予輸入/輸出支持。一個常見的例子就是針對數(shù)值波形求解程序的輸出格式。pkg-utils的插件系統(tǒng)會將這些額外的格式與ObsPy的其他部分無縫地集成到一起。

1.5格式自動檢測與用法

ObsPy有一個頂層函數(shù)read()，是一個讀入波形數(shù)據(jù)時的單一入口點。見列表2中的用法舉例。

read(filename，**kwargs) 程序調(diào)用所有注冊格式的is_fileformat()函數(shù)，直到有一個返回True。根據(jù)正在探尋的格式，is_fileformat()程序解析開始的幾個字節(jié)或執(zhí)行更復雜的經(jīng)驗方法。出于執(zhí)行效率的原因，格式檢測程序的速度必須要快。于是，ObsPy按照人工挑選的清單進行格式檢測，這樣最常用的格式將首先檢測，從而提高了平均執(zhí)行效率。在格式得到確認后，就會調(diào)用相應格式的read_fileformat()函數(shù)，對文件進行解析并返回Stream對象。也可以直接為read()程序提供格式參數(shù)，這樣可以省略格式檢測程序。

通過將某一功能實現(xiàn)為read()程序的一部分，這種結(jié)構(gòu)允許將此功能對所有格式進行共享。例如當檢測到有效的HTTP通用資源地址時，文件會自動被下載，以及一系列不同的檔案格式的解壓。

世界各地數(shù)據(jù)中心提供的Web服務是波形數(shù)據(jù)的重要來源。ObsPy實現(xiàn)客戶端能夠與諸多數(shù)據(jù)中心[20]進行交互操作。一條波形數(shù)據(jù)的請求將以其被存入Stream對象結(jié)束，所以無論數(shù)據(jù)來源于何處，得到數(shù)據(jù)后的工作流程都是相同的。

1.6特定專業(yè)領(lǐng)域的便捷方法

統(tǒng)一的內(nèi)部數(shù)據(jù)表示開啟了定義數(shù)據(jù)轉(zhuǎn)換方法的可能性。對此，ObsPy提供了全面的信號處理程序集合，這些地震學中經(jīng)常使用的程序，非常依賴于NumPy和SciPy提供的功能。它們的實現(xiàn)既滿足了地震學家對某個處理操作的需要，同時也保持數(shù)據(jù)自身的一致性。例如Trace.decimate()方法會應用濾波對數(shù)據(jù)進行抽采，然后修改時間序列的元數(shù)據(jù)的采樣率。由于采樣點數(shù)可能會很多，所以大多數(shù)操作是按照對Stream和Trace對象進行原地修改而實現(xiàn)的。

列表2一段交互的Python會話的截屏，用以表明ObsPy中read()函數(shù)的用法。它將會檢測文件格式并調(diào)用相應的讀入程序。如果read()流程檢測到有效的HTTP通用資源地址(URL)，它就會在讀入之前先下載該資源。如果檢測到存檔格式，它就先進行解壓

>>>importobspy>>>st=obspy.read("filename")>>>st>>>printst3Trace(s)inStream:BW.RJOB..EHZ∣2009-08-24T00:20:03.000000Z-...∣100.0Hz,3000samplesBW.RJOB..EHN∣2009-08-24T00:20:03.000000Z-...∣100.0Hz,3000samplesBW.RJOB..EHE∣2009-08-24T00:20:03.000000Z-...∣100.0Hz,3000samples表1 Stream和Trace對象的處理方法的摘選。大多數(shù)Trace的方法也可用于Stream對象,就是將其應用Stream對象包含的所有子對象Stream方法merge() 合并具有相同ID的Trace對象rotate()旋轉(zhuǎn)兩分量或三分量的Stream對象select()選出滿足條件的Trace對象組成的新Stream對象Trace方法decimate()根據(jù)整數(shù)因子對數(shù)據(jù)抽采detrend()去除數(shù)據(jù)中的線性趨勢differentiate()數(shù)據(jù)相對于時間的微分filter()對數(shù)據(jù)進行濾波integrate()數(shù)據(jù)相對于時間的積分normalize()將數(shù)據(jù)歸一化到它的絕對峰值remove_response()反卷積儀器響應resample()在頻率域?qū)?shù)據(jù)進行重采樣taper()數(shù)據(jù)兩端尖滅截斷trigger()對數(shù)據(jù)運行觸發(fā)算法trim()由指定的起始時間和結(jié)束時間來截取數(shù)據(jù)

多數(shù)處理方法是對單個Trace定義的，對Stream對象施加同樣的方法將會對其每個子Trace對象調(diào)用相應的方法。這就允許對三分量甚至更多分量的數(shù)據(jù)進行處理。一些方法是只針對Stream對象的。參考表1選擇可行的信號處理程序。

為了提供多種不同函數(shù)，ObsPy又一次使用了插件系統(tǒng)。列表3～5以Trace.taper()方法為例表明了這一點。所用的插件系統(tǒng)通過將不同模塊定義的功能集成進簡單的科學家用的API特定領(lǐng)域的應用程序接口之中，進而，它還鼓勵了重復代碼剔除和代碼再利用，從而使得ObsPy能夠提供更多種不同函數(shù)，而不用實現(xiàn)和測試大量實例中的細節(jié)。

2舊代碼的集成

對于某些特定的任務，地震學界經(jīng)常依靠已經(jīng)使用十幾年甚至更久的代碼。這種大量公開和應用于各類問題，保證了這些代碼在大多數(shù)情況下能按照期望來工作，因為大部分問題早已被發(fā)現(xiàn)和修復。雖然這些代碼常常不便使用，但仍然希望使其保持活力并從過去開發(fā)的投入中獲益。ObsPy集成了一些舊代碼，使其能夠使用現(xiàn)代工作流程和新近發(fā)展的數(shù)據(jù)處理方法，同時依靠穩(wěn)定的、經(jīng)過充分測試的代碼完成特定功能。這部分內(nèi)容將通過兩個例子給予展示。

列表3摘自O(shè)bsPy的setup.py文件。安裝完成后，這些入口點會進行注冊并使pkg-utils可見。下面的代碼片段通過注冊不同的兩端尖滅窗來展示這一點。在文本寫作的時候，ObsPy顯示了18個不同的尖滅窗，它們主要來自scipy.signal模塊。用法說明見列表5

..."obspy.plugin.taper":[ "cosine=obspy.signal.invsim:cosTa-per", "barthann=scipy.signal:barthann", ...]...

列表4來自Trace()對象內(nèi)部方法taper()的代碼片段。它顯示了表3中展示的注冊過的功能如何調(diào)用，從而獲得不同的兩端尖滅窗，以及可選的關(guān)鍵變量是如何進行傳遞的。用法展示見列表5

...retrievefunctioncallfromentrypointsfunc=_getFunctionFromEntryPoint( "taper",type)...taper_sides=func(2*wlen,**kwargs)...

列表5列表3和列表4中注冊過和調(diào)用過的函數(shù)的用法。注意實際的兩端尖滅窗是如何形成的，在這里它們是在Python不同模塊中定義的。在此基礎(chǔ)上，此表還顯示了一套鏈式方法調(diào)用，這之所以可能是因為這些方法會返回相應的對象。copy()的調(diào)用是必要的。由于大多數(shù)方法是原處處理方法，因此可能會改變對象。大多數(shù)情況下，這個節(jié)省內(nèi)存的方法是需要的，它是為適應通常使用情況而有意設(shè)計的

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計算走時

地震走時的計算是地震學中的常見問題。當?shù)卣鸢l(fā)生在X點，我們的問題是，什么震相在何時到達Y點。使用整個地球的全三維模擬[29]的成本對于許多應用來說都過高，用精度較低但更快的近似方法常常就足夠了。如果一維球狀均勻分層地球模型是適用的，那么常用的Buland和Chapman[3]方法即可非?？焖俚赜嬎阕邥r。一個通常被稱作iaspei-tau[14，26]的Fortran程序包已經(jīng)實現(xiàn)了這個方法。功能更多的Java版本[4]也有。obspy.taup子模塊使用了Python標準庫里的ctypes來提供iaspei-tau的Python接口，見列表6。圖2繪制了一個地震波穿過地球的走時，橫坐標以度為單位的地理距離。

圖2　使用一維地球模型ak135的地震波走時(原圖為彩色圖——譯注)。上圖為由obspy.taup模塊計算的一些地震波震相走時。下圖顯示了由TauP工具包[4]和obspy.taup模塊分別計算的P波走時差。該差值主要是由內(nèi)部坐標系不同帶來的，地震學界早已了解[15]。偏差結(jié)果強烈地依賴于所計算的地震震相和震源-接收點的幾何位置

列表6由obspy.taup得到震相走時的使用方法，震相由震中距25°，深度10km的一個震源所產(chǎn)生。震中距(delta)單位為球狀地球模型假設(shè)下的度數(shù)，深度(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

遍布全球的地震接收儀器是為了盡可能精確地測量并記錄地震動。影響最終波形特性的因素有很多，其中包括物理儀器的頻率響應、任何放大器的效應、模擬或數(shù)字濾波器以及數(shù)字化帶來的影響。對許多研究地球的應用來說，為了盡可能地估計真實地震動，去除這些因素帶來的影響是相當重要的。因為地震接收儀器和處理過程的影響而進行數(shù)據(jù)校正的第一步就是計算記錄系統(tǒng)的頻率響應。之后用地震圖對該頻率響應進行反褶積，從而能得到具有物理單位、且在研究頻帶內(nèi)不受儀器效應影響的地震圖。這個過程在地震學中被稱為儀器校正。

地震的記錄系統(tǒng)是用由不同階段或元素所組成的線性鏈來描述的。SEED數(shù)據(jù)格式[18]可以精確地表示這個流程，這已經(jīng)被地震學界廣泛地接受并從1990年代初開始使用。最近，SEED的指定接替者StationXML[19] 已經(jīng)被開發(fā)出來，并且已經(jīng)開始得到地震學界的接受。作為XML格式，就支持工具、方便的數(shù)據(jù)發(fā)布和人工可讀性方面來說，它比二進制的SEED格式更有優(yōu)勢。除了一些小的變化外，這兩種格式儲存的信息是相同的。

列表7C結(jié)構(gòu)體表示由零點和極點組成的儀器物理響應。用ctypes再次創(chuàng)建此內(nèi)容需要定義如列表8中所示的Python類

structcomplex{doublereal;doubleimag;};structpole_zeroType{intnzeros;intnpoles;doublea0;doublea0_freq;structcomplex*zeros;structcomplex*poles;};

列表8用ctypes庫創(chuàng)建的Python類，用來表示列表7中所示的零、pole_zeroTypeC結(jié)構(gòu)體。本例中的復數(shù)(complex_number)已經(jīng)定義過

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文件計算地震接收系統(tǒng)頻率響應的標準流程包括用rdseed[12]生成RESP文件，該文件是一個文本的SEED的真子集。然后將這些RESP文件輸給evalresp[11]，得到頻率響應的結(jié)果。之所以要轉(zhuǎn)而通過RESP文件得到結(jié)果，是因為它是evalresp所能接受的唯一一種輸入格式。用存儲于StationXML文件中的元數(shù)據(jù)進行儀器校正還包括另外一個步驟，即用另外一個工具將StationXML文件轉(zhuǎn)換成SEED文件。

在Python中實施這一工作包括大量系統(tǒng)調(diào)用和不必要的輸入/輸出運算，這使它變得不適用于現(xiàn)代化大型數(shù)據(jù)工作流程。因為在計算儀器響應的時候存在很多可能會對最終結(jié)果產(chǎn)生很大影響的隱患，所以取代evalresp將成為主要的努力方向。而ObsPy則提供了從StationXML(通過lxml解析)到evalresp的直接通道，借助于ctypes調(diào)用evlaresp內(nèi)部函數(shù)。這種方法對于那些仍在積極開發(fā)中的代碼并不可行，因為內(nèi)部應用程序接口的穩(wěn)定性得不到保證。不過evalresp除了需要偶爾的維護，并不再進行主動的開發(fā)。

ObsPy對StationXML數(shù)據(jù)的內(nèi)部表示用于導出嵌套的C結(jié)構(gòu)體，符合evalresp內(nèi)部函數(shù)的期望。它們是用ctypes按照列表7和8中表示的那樣進行定義和初始化的。Ctypes還需要如列表9和10中所列函數(shù)頭的定義。圖3是由StationXML得到的振幅和相位響應，同時展示了記錄系統(tǒng)不同階段的效應。

地震接收儀的記錄系統(tǒng)通常包括許多步驟，最終在evalresp內(nèi)部生成一個相當復雜的表示。evalresp的文件解析步驟負責將SEED結(jié)構(gòu)翻譯成相應的內(nèi)部結(jié)構(gòu)，有時候做出一些不明顯的決定。為確保evalresp集成無誤，我們進行了廣泛的測試[16]。我們定義了作用在StationXML文件上的集成evalresp的正確性，即如果最終的響應和由StationXML轉(zhuǎn)到SEED，再由SEED轉(zhuǎn)換得到RESP文件，用evalresp對RESP文件進行計算的結(jié)果相同，就證明集成是正確的。我們從最大的地震學數(shù)據(jù)發(fā)布機構(gòu)美國地震學聯(lián)合研究協(xié)會下載了幾乎所有可用的StationXML庫存數(shù)據(jù)。這些數(shù)據(jù)來自27 000多個臺站，其中大多數(shù)具有被不同時間周期所定義的多個記錄通道，這樣就有超過100 000個儀器響應。我們反復執(zhí)行程序，直到通過測試，以此證明我們的解決方案可以可靠地應用于各地地震社團中遇到的任何數(shù)據(jù)。

圖3　XM.05..HHZ通道的振幅和相位響應，用來將速度(m·s-1)記錄的地震圖轉(zhuǎn)換為數(shù)字計數(shù)(原圖為彩色圖——譯注)。記錄儀為GuralpCMG-6TD地震計。綠線表示只考慮物理儀器和模數(shù)轉(zhuǎn)換增益的響應；黑線除此以外還合并了3個相繼的數(shù)字抽樣步驟；紅線是這一通道的柰奎斯特頻率；綠線和黑線結(jié)果是用evalresp的ObsPy整合工具計算的。為了進行對比，我們用橘黃色線表示由JEvalResp[13]計算的相應結(jié)果，evalresp的Java版本是很少的能夠進行這個計算的軟件包之一。結(jié)果一致性非常好，奈奎斯特頻率后面的差別是由繪圖所用離散頻點的間距引起的

列表9evalrespC代碼中的函數(shù)說明。這是用于計算響應的主函數(shù)。注意輸入如何使用了各種不同的參數(shù)，從簡單的整數(shù)到自定義結(jié)構(gòu)的數(shù)組。complex的結(jié)構(gòu)體在Python里并不需要定義，因為numpy.complex128的dtype具有相同的內(nèi)存布局

...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模塊提供的數(shù)組類型會生成自動的類型、維數(shù)和函數(shù)調(diào)用的檢測標志。這產(chǎn)生出在調(diào)用共享庫之前Python一側(cè)的便捷的調(diào)用語法和錯誤處理

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目標

與上述那些已經(jīng)久負盛譽、幾乎無漏洞的舊代碼不同，許多科學家將Python用于他們的分析工作流程，傾向于開發(fā)零散代碼用于他們的輸入數(shù)據(jù)，因此第一眼看上去對其他人沒什么用。當輸入數(shù)據(jù)有不同的格式，而輸出結(jié)果需要與不同下游處理軟件兼容等情況存在時，這個一般性斷言在地震學中也是真的。本節(jié)中，我們首先給一個案例，研究如何翻譯一個能讀取奇異格式Win的C代碼，然后呈現(xiàn)一個博士論文期間開發(fā)的科研代碼，該代碼現(xiàn)在已作為一個新模塊被整合進ObsPy中。由于ObsPy的插件功能，只有必要的處理步驟需要編寫。

3.2讀取奇異的格式——obspy.win

WIN是由Hakusan日本公司構(gòu)想的一種格式，是他們Datamark數(shù)據(jù)采集器系列的缺省存儲格式。一分鐘的WIN文件數(shù)據(jù)是高度壓縮的，每秒的數(shù)據(jù)和前一秒相比都有不同程度的壓縮。截至目前，有3種方法可以將WIN格式轉(zhuǎn)換為更標準的SAC格式，包括Linux下的Win2Sac和japan2sac程序，以及Windows下的有圖形用戶界面的Win2Sac程序。Win2Sac的代碼來自日本，而japan2sac是Chambery(法國)一個科研項目的一部分，尚未全部完成。為了用最新的程序和軟件來處理這些數(shù)據(jù)，例如用基于ObsPy的腳本，使用者還得將全部波形存檔轉(zhuǎn)換為SAC格式，從而可以用Obs-Py讀取。這個過程的每一步都使硬盤上存檔的體積加倍，同時也可能在最終產(chǎn)出結(jié)果中加入了錯誤。

幸好Win2Sac的源代碼[22]可在線獲得，同時還提供了描述WIN格式的數(shù)據(jù)表單[30]，我們已經(jīng)可以寫ObsPy的插件來對其進行讀取。最終的目標是能夠直接不必重復地從波形存檔中讀取數(shù)據(jù)，并且與處理其他格式的數(shù)據(jù)一樣來處理這些數(shù)據(jù)。即使有的目標是將WIN轉(zhuǎn)換成MiniSEED格式，但是現(xiàn)在最好的解決辦法還是用Obs-Py寫Python腳本，而不再依靠一系列格式轉(zhuǎn)換步驟。

3.3將科研代碼改寫成ObsPy的擴展模塊

一座火山的狀態(tài)及其活躍性可以被其側(cè)翼布設(shè)的地震計記錄到的振幅大小來表明。實時地震振幅測量[RSAM(1)]最早是由Endo和Murray[8]引入到火山噴發(fā)預測和火山活動性評估的工作中。我們也可以計算一定時間段內(nèi)信號的振幅平方標準偏差，或者其實時地震能量測量[RSEM(2)][5]，再除以觀測個數(shù)：

(1)

(2)

圖4　在一個火山活動比較安靜的時期，PSG臺在5～20Hz頻帶內(nèi)的地震振幅(RSAM)每日和每周的漲落。星期五(做禮拜的日子)的變化顯而易見。振幅單位為counts(原圖為彩色圖——譯注)

圖5　從火山活動性比較活躍的時期到比較安靜的時期火山口附近(POS臺和PUN臺站)以及火山登山露營基地附近臺站(PSG臺站)記錄到的過渡變化。盡管計算的是記錄第十個百分點的振幅，但從PSG臺站仍然能看到人類活動的影響(原圖為彩色圖——譯注)

式中Ai是地震信號的振幅，Aavg是N個采樣點的平均振幅。

當多個震源重疊時，火山地震學家通常是在事先定義好的頻帶內(nèi)分別計算RSEM或RSAM值，即地震能量譜測量(SSEM)[28]，類似于Stephens[27]推薦的SSAM。首先對數(shù)據(jù)進行去平均和二階巴特沃思帶通濾波。然后在30s的時窗(100×30個采樣點)內(nèi)計算每個值。最后計算每天第10個和第25個百分點處的值，以及中位值，用以部分去除不需要的臨時信號和地震事件，它們被看成是火山脈動的干擾[7]。

obspy.signal.ssxm子模塊中實現(xiàn)了RSAM，RSEM，SSAM和SSEM算法。最終的代碼利用了Pandas(Python數(shù)據(jù)分析庫)的重采樣功能，完成后SSXM()方法代碼的總行數(shù)不到50行。

3.3.1事例結(jié)果

用不同頻帶對信號進行濾波的好處可以用卡瓦伊真火山(印度尼西亞東爪哇)地震數(shù)據(jù)進行描述。地震時間序列通常會被人類日間活動所污染，特別是在高頻段(>10Hz)以及那些靠近“人文”噪聲源的臺站，比如停車場和火山登山的露營基地。我們甚至能觀測到每日的工作活動在做禮拜的日子變?nèi)?即圖4中的星期五)。為了最小化瞬時信號的影響，人們可以在特定的頻帶內(nèi)(圖5)使用第10個百分點處(p10)的振幅。除了火山活躍期(圖5)，與POS和PUN臺站相比，PSG臺站的p10值仍然被人類活動噪聲嚴重地影響著，這兩個臺站位于火山口附近，遠離人類活動。

其他研究還通過對地震信號進行濾波而調(diào)查到一些特殊的過程。比如如果任意尺度上的持續(xù)破裂都會產(chǎn)生持續(xù)的地震信號[6]，那么SSEM的計算可能會提供應變釋放率的測量。人們還可以對微震噪聲源(0.1～1Hz)的平穩(wěn)性進行評定來估計通過由地震環(huán)境噪聲互相關(guān)技術(shù)[17]得到的速度變化的結(jié)果。

4結(jié)論

舊代碼和數(shù)據(jù)格式廣泛應用于各個學科，在可見的未來，它們將繼續(xù)發(fā)揮重要的作用。本文所提出的方法有助于將它們集成到現(xiàn)代化的計算機環(huán)境中。文中展示的大多數(shù)設(shè)計和實現(xiàn)時的選擇都可以轉(zhuǎn)移到計算科學的其他領(lǐng)域。從復雜的文件格式集合中抽取出它們的共性，這使得構(gòu)建數(shù)據(jù)源獨立的工作流程成為可能，尤其有利于數(shù)據(jù)繁重的領(lǐng)域。同時結(jié)合使用特定專業(yè)領(lǐng)域詞匯的功能庫，包括一些緊要的舊代碼，為那些甚至不熟悉Python的人提供了誘人的軟件包。

ObsPy在地震學界享有廣泛的認同度。我們相信這主要歸功于它解決了一個實際的問題，那就是大量不同文件格式之前需要過多的文件格式轉(zhuǎn)換工具或不同的信號處理工具。一旦人們開始使用ObsPy，他們將會很快發(fā)現(xiàn)Python的靈活和強大。相較MATLAB而言，使用ObsPy還能獲得全功能編程語言的優(yōu)勢。另一個有競爭力的特性是能夠使用許多Python中第三方模塊，雖然它們并不與信號處理直接相關(guān)，但能夠使用數(shù)據(jù)庫、Web服務、機器學習庫，當然還有最近在處理大數(shù)據(jù)集方面的新進展，它們在地震學領(lǐng)域和其他學科內(nèi)將變得越來越重要。另外，Obspy的完整包是免費開源的，并且能夠在幾乎所有相關(guān)平臺上運行。

成功使用ObsPy的研究包括事件重定位[21]、旋轉(zhuǎn)[10]和時變[24]的地震學、大數(shù)據(jù)處理[1]，以及綜合研究全波形反演[25]和衰減核[9]等。ObsPy的網(wǎng)站可通過鏈接http://www.obspy.org訪問，其中包括詳細的文檔、詳盡的教程和一個郵件列表，該列表的目的是為了使轉(zhuǎn)換到ObsPy的過程更為容易，以及建立關(guān)于ObsPy的社交圈。本文中所描述的模塊化和測試驅(qū)動的研發(fā)有利于添加新功能。借此，ObsPy不斷地增加外部貢獻的個數(shù)，目標是變?yōu)橛烧麄€業(yè)界的人士共同維護的代碼。

參考文獻

[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，F(xiàn)alsaperla 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-)，女，中國地震臺網(wǎng)中心工程師，主要從事地震觀測數(shù)據(jù)處理工作。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生態(tài)系統(tǒng)的橋梁.世界地震譯叢.47(4):344-357.doi:10.16738/j.cnki.issn. 1003-3238.201604006

中國地震臺網(wǎng)中心韓雪君譯；馬延路校

中國地震局地球物理研究所魯來玉復校