查小惠 張廣偉 楊雪超

1）中國(guó)南昌 330039 江西省地震局

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

0 引言

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

1 程序運(yùn)行平臺(tái)和功能

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

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

2 程序運(yùn)行流程

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

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

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

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

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

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

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

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

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

3 關(guān)鍵子程序

3.1 arcfetch

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

3.2 ref2sac

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

3.3 taup_time

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

3.4 saclab程序包

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

3.5 datestr 和datenum

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

4 結(jié)束語(yǔ)

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

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

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

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

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