查小惠 張廣偉 楊雪超

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

程序?qū)嶋H為一個Matlab函數(shù)，調(diào)用該函數(shù)需要輸入4個參數(shù)，分別為P波（實際為初至波）到達(dá)前時間、截取波形的總長度、相對于臺站的最小震中距和最大震中距（用來篩選地震目錄）。截取事件波形需要一個參考時間，本程序中參考時間是P波（初至波）到達(dá)臺站的理論到時，函數(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)緯度信息，可以對地震目錄進行挑選，選擇標(biāo)準(zhǔn)是調(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ù)，并存放在輸出文件夾下。當(dāng)該程序運行正確時，返回值為0，如果返回值為其他值，表明數(shù)據(jù)格式轉(zhuǎn)換不成功。

3.3 taup_time

波形截取時，以理論P波到達(dá)臺站的時間點為截取波形的參考時間，選擇截取該參考時間點前后各多少秒的數(shù)據(jù)。TAUP軟件是目前地震學(xué)中計算地震波走時的常用軟件（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ù)值表示當(dāng)前時間到該起點以天為單位的時間距離，該數(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），進一步證明了該程序的可靠性。

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

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

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

查小惠， 雷建設(shè). 云南地區(qū)地殼厚度和泊松比研究[J]. 中國科學(xué)：地球科學(xué)， 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.