REFTEK連續(xù)波形截取程序

查小惠 張廣偉 楊雪超

1）中國南昌 330039 江西省地震局

2）中國北京100085 中國地震局地殼應(yīng)力研究所

0 引言

近年來，流動地震觀測越來越多，流動地震觀測數(shù)據(jù)的存儲管理和使用問題，受到越來越多的關(guān)注。2010年4月到2011年7月，中國地震局地殼應(yīng)力研究所在云南地區(qū)布設(shè)21個寬頻帶野外流動觀測臺站，所用地震計為CMG-3ESPC，數(shù)據(jù)采集器為REFTEK-130B，最終獲得233.7 G連續(xù)波形資料（Lei et al, 2012）。在資料處理階段，由于研究目的和手段的差異，有時需要采用不同的地震目錄，選用不同的地震參數(shù)，反復從連續(xù)波形數(shù)據(jù)中截取所需事件波形資料?；谘芯啃枰?，嘗試開發(fā)一套簡易的REFTEK連續(xù)波形截取程序，以便從REFTEK格式的連續(xù)波形數(shù)據(jù)中截取任意長度的SAC格式的地震事件波形。

1 程序運行平臺和功能

REFTEK連續(xù)波形截取程序可以在個人計算機上安裝運行，基于Linux系統(tǒng)下的Matlab平臺開發(fā)，因為：①linux系統(tǒng)下有現(xiàn)成的可以截取REFTEK連續(xù)波形數(shù)據(jù)的arcfetch程序和將REFTEK格式數(shù)據(jù)轉(zhuǎn)換為SAC格式數(shù)據(jù)的ref2sac程序；②Matlab平臺有很多高級函數(shù)，比如：datenum和datestr函數(shù)可以很好地解決波形截取過程中遇到的時間運算問題。在處理SAC文件上，也有現(xiàn)成的saclab工具包可以使用，而且在Matlab中可以利用Unix命令方便調(diào)用Linux系統(tǒng)的程序?；贚inux系統(tǒng)下的Matlab平臺編寫程序，可以同時利用Matlab高級函數(shù)和Linux系統(tǒng)程序兩者的優(yōu)勢，高效實現(xiàn)程序功能。

本程序的功能是從REFTEK連續(xù)波形數(shù)據(jù)中截取事件波形，原始數(shù)據(jù)為流動觀測記錄的REFTEK格式的連續(xù)波形數(shù)據(jù)，截取得到SAC格式的三分量事件地震波形數(shù)據(jù)。目前程序的測試平臺為Ubuntu12.04和Matlab2009b，只要確保Matlab相關(guān)程序包的安裝和Linux系統(tǒng)子程序的成功調(diào)用，程序即可順利運行。程序的運行時間和地震目錄的個數(shù)及截取波形的長度有關(guān)，雖然調(diào)用了較多外部程序，但是程序整體運行效率可以接受的。

2 程序運行流程

在調(diào)用Matlab程序截取波形前，需要進行相關(guān)數(shù)據(jù)準備。首先是原始連續(xù)波形數(shù)據(jù)，程序要求每個臺站記錄的REFTEK原始波形數(shù)據(jù)單獨存放在一個文件夾下，而且為了方便訪問該文件夾，使用阿拉伯數(shù)字給文件夾編號命名。然后是地震目錄的整理，地震目錄不限行數(shù)，每行8列，1到3列分別為年、月、日，第4列為地震發(fā)震時刻，格式為時分秒，5到8列分別為震源的緯度、經(jīng)度、深度和震級大小，各列之間以空格分開。最后是準備臺站信息，包括4列，分別為臺站編號、臺站經(jīng)度、臺站緯度和臺站高程。臺站信息文件直接控制截取哪些臺站的事件波形，如果不需要截取相關(guān)臺站的波形，則將該臺站的信息從臺站信息文件中刪除。臺站信息文件中的臺站編號應(yīng)該和存放該臺原始連續(xù)波形的文件夾命名編號一致，方便在波形截取中尋找對應(yīng)臺站的數(shù)據(jù)文件。

程序?qū)嶋H為一個Matlab函數(shù)，調(diào)用該函數(shù)需要輸入4個參數(shù)，分別為P波（實際為初至波）到達前時間、截取波形的總長度、相對于臺站的最小震中距和最大震中距（用來篩選地震目錄）。截取事件波形需要一個參考時間，本程序中參考時間是P波（初至波）到達臺站的理論到時，函數(shù)輸入的第1個參數(shù)是相對于該參考時間的提前時間，結(jié)合截取波形的總時間長度，計算截取波形的時間終點。程序運行流程見圖1。

圖1 程序流程Fig.1 Program fl ow diagram

（1）從存放臺站信息的文件中，讀入臺站信息，包括臺站編號、臺站經(jīng)緯度和高程。然后從第1行開始循環(huán)，即波形數(shù)據(jù)的截取是逐臺進行的。

（2）讀入地震目錄，結(jié)合臺站經(jīng)緯度信息和地震目錄中地震事件的經(jīng)緯度信息，可以對地震目錄進行挑選，選擇標準是調(diào)用主程序時輸入的最小和最大震中距。如果截取的事件波形用于接收函數(shù)研究，一般將震中距限制在30°—95°。

（3）選取地震目錄后，計算到時目錄。主要計算需要截取波形的前后兩個時間點的時間。利用TAUP軟件計算P波理論走時，發(fā)震時刻加上P波理論走時就得到P波理論到時，即波形截取的參考時間。然后根據(jù)調(diào)用函數(shù)時輸入的截取波形提前時間和截取波形總時間長度，計算截取波形前后兩個端點的時間。

（4）從第1個地震目錄開始，使用arcfetch程序，從REFTEK連續(xù)波形中截取已計算兩時間點間的事件波形，使用ref2sac程序，將REFTEK格式數(shù)據(jù)轉(zhuǎn)換為SAC格式數(shù)據(jù)。

（5）使用saclab程序包處理SAC數(shù)據(jù)，主要重寫SAC數(shù)據(jù)頭文件，需要注意SAC頭文件中的時間處理。

（6）利用相關(guān)文件處理命令，對SAC文件進行存儲，并返回流程（4）截取下一個地震事件，地震目錄截取完畢，返回流程（1）截取下一個臺站。最后獲得的SAC數(shù)據(jù)仍然以臺站為單位進行存放，也就是，同一個臺站截取的不同地震事件的SAC文件存放在同一個文件夾。如果想改為一個地震事件文件夾下存放多個臺站文件的數(shù)據(jù)格式，需要編寫相應(yīng)的腳本。

3 關(guān)鍵子程序

3.1 arcfetch

arcfetch程序用來從原始的連續(xù)REFTEK波形數(shù)據(jù)中截取一段波形數(shù)據(jù)，波形截取的起始時間和波形長度可以在程序輸入?yún)?shù)中指定，且截取數(shù)據(jù)格式仍然為REFTEK數(shù)據(jù)格式。使用該程序時，需要指定原始連續(xù)波形數(shù)據(jù)的存放路徑。為保證程序成功執(zhí)行，在存放連續(xù)波形數(shù)據(jù)的文件夾中需要一個archive.sta文件，該文件可以通過arccreate程序生成。以上所述兩個程序和程序使用的說明文檔可以在REFTEK官網(wǎng)下載。arcfetch程序截取連續(xù)波形會有3種結(jié)果：①截取到完整正確的REFTEK數(shù)據(jù)文件；②沒有截取到REFTEK數(shù)據(jù)文件；③截取到REFTEK數(shù)據(jù)文件，但數(shù)據(jù)為空。對于后兩種情況需編寫相關(guān)測試程序，以避免后續(xù)的ref2sac程序出錯。

3.2 ref2sac

使用arcfetch截取得到REFTEK事件波形，需要使用ref2sac程序轉(zhuǎn)換為SAC格式。該程序為PASSCAL（大陸巖石圈地震臺陣研究計劃）提供的rpm軟件包中的一個子程序，可以將該子程序包在Linux系統(tǒng)下單獨編譯使用。該程序有兩個參數(shù)，分別是輸入文件夾和輸出文件夾。它將輸入文件夾下REFTEK格式數(shù)據(jù)轉(zhuǎn)換為SAC格式數(shù)據(jù)，并存放在輸出文件夾下。當該程序運行正確時，返回值為0，如果返回值為其他值，表明數(shù)據(jù)格式轉(zhuǎn)換不成功。

3.3 taup_time

波形截取時，以理論P波到達臺站的時間點為截取波形的參考時間，選擇截取該參考時間點前后各多少秒的數(shù)據(jù)。TAUP軟件是目前地震學中計算地震波走時的常用軟件（Crotwell et al，1999）。本程序利用taup_time計算地震目錄中地震事件傳播到臺站的理論走時，將理論走時加上發(fā)震時刻得到理論P波到時，作為波形截取的參考時間。具體實現(xiàn)過程中發(fā)現(xiàn)，并不是所有震中距都是P波為初至波，實際上參考時間應(yīng)該是理論初至波到時。在使用taup_time程序時，只給定參考模型和震中距參數(shù)，比時TAUP會計算出所有震相走時，而且初至波走時顯示在第1行，本程序提取第1行的震相走時作為實際使用的走時。所以，文中所說P波和初至波是等價的。

3.4 saclab程序包

ref2sac程序只是將REFTEK格式的地震數(shù)據(jù)轉(zhuǎn)換為SAC格式，但是SAC數(shù)據(jù)的頭文件并不正確，所以在轉(zhuǎn)換得到SAC格式地震數(shù)據(jù)后，需要進一步修改SAC數(shù)據(jù)的頭文件。saclab是由Andreas Wustefeld使用Matlab語言編寫的程序包，便于SAC文件讀寫和頭文件修改。需要注意的是，O、B、A、E這4個SAC頭文件，因為涉及到波形數(shù)據(jù)的時間信息，一旦錯誤，后續(xù)利用該SAC波形數(shù)據(jù)進行的相關(guān)走時研究都將會出現(xiàn)系統(tǒng)性錯誤。

3.5 datestr 和datenum

Matlab中表示日期的時間有3種格式：日期字符串、連續(xù)的日期數(shù)值和日期向量（王正林等，2007）。日期字符串格式使用較多，也較為直觀，如2013-10-1 10∶10∶10，但是不方便進行時間的加減運算。連續(xù)的日期數(shù)值格式以公元元年1月1日為起點的，用一個數(shù)值表示當前時間到該起點以天為單位的時間距離，該數(shù)值為浮點數(shù)，所以方便進行時間的加減運算。Matlab中的datenum函數(shù)可以實現(xiàn)將日期字符串和日期向量格式轉(zhuǎn)化為連續(xù)的日期數(shù)字格式，而datestr函數(shù)可以實現(xiàn)把日期數(shù)字和日期向量格式轉(zhuǎn)換成日期字符串格式，波形截取程序中，需要進行儒略日的計算，利用datenum函數(shù)使用一條語句就可以實現(xiàn)， julianday=datenum（year，month，day）-datenum（year，1，1）+1。波形截取程序需要進行時間的加減運算，以確定截取波形的時間起始點，但arcfetch在截取波形時需要輸入的時間格式是日期字符串格式，因為日期字符串表示的時間進行加減運算比較困難。本程序采取的方法是，在Matlab中，先將日期字符串格式轉(zhuǎn)化為連續(xù)的日期數(shù)值格式，然后進行相關(guān)的時間運算，再轉(zhuǎn)換回日期字符串的格式，便于時間運算。

4 結(jié)束語

REFTEK連續(xù)波形截取程序，是一個Matlab函數(shù)，功能就是從REFTEK連續(xù)波形數(shù)據(jù)中截取事件波形。該程序連接多個重要的子程序，并使用相應(yīng)的Matlab工具包，分別完成特定功能，利用Matlab的高端平臺特性，以簡單編程高效完成波形截取任務(wù)。在Linux系統(tǒng)下，利用Matlab平臺連接各種處理程序的編程方式，減少編程時間，快速完成特定功能，但程序的移植性和易用性較差，不過作為個人的科研使用仍然具有一定使用價值。利用該程序截取得到的事件波形，可以進行震相到時讀取，進而開展走時相關(guān)研究工作（黎源等，2012），也可以利用三分量波形數(shù)據(jù)開展包括接收函數(shù)（查小惠等，2013）在內(nèi)的一系列波形研究工作。目前使用該程序截取的波形所開展的相關(guān)研究工作，得到了較好的結(jié)果（Lei et al，2012；黎源等，2012； 查小惠等，2013），進一步證明了該程序的可靠性。

程序編寫得到中國科學院青藏高原研究所劉啟民博士和中國地震局地殼應(yīng)力研究所孫長青博士的幫助，在此表示感謝。

黎源， 雷建設(shè). 青藏高原東緣上地幔頂部Pn波速度結(jié)構(gòu)及各向異性研究[J]. 地球物理學報， 2012， 55（11）：3 615-3 624.

王正林， 劉明. 精通MATLAB7[M].北京：電子工業(yè)出版社， 2007.

查小惠， 雷建設(shè). 云南地區(qū)地殼厚度和泊松比研究[J]. 中國科學：地球科學， 2013， 43（3）： 446-456.

Crotwell H P， Owens T J and Ritsema J. The Taup ToolKit： Flexible Seismic Travel-Time and Raypath Utilities[J]. Seismological Research Letters， 1999，70（2）：154-160.

Lei J S， Zhang G W， Xie F R， et al. .Relocation of the 10 March 2011 Yingjiang， China， earthquake sequence and its tectonic implications [J]. Earthquake Science， 2012， 25： 103-110.