REFTEK连续波形截取程序

查小惠 张广伟 杨雪超

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

2）中国北京100085 中国地震局地壳应力研究所

0 引言

近年来，流动地震观测越来越多，流动地震观测数据的存储管理和使用问题，受到越来越多的关注。2010年4月到2011年7月，中国地震局地壳应力研究所在云南地区布设21个宽频带野外流动观测台站，所用地震计为CMG-3ESPC，数据采集器为REFTEK-130B，最终获得233.7 G连续波形资料（Lei et al, 2012）。在资料处理阶段，由于研究目的和手段的差异，有时需要采用不同的地震目录，选用不同的地震参数，反复从连续波形数据中截取所需事件波形资料。基于研究需要，尝试开发一套简易的REFTEK连续波形截取程序，以便从REFTEK格式的连续波形数据中截取任意长度的SAC格式的地震事件波形。

1 程序运行平台和功能

REFTEK连续波形截取程序可以在个人计算机上安装运行，基于Linux系统下的Matlab平台开发，因为：①linux系统下有现成的可以截取REFTEK连续波形数据的arcfetch程序和将REFTEK格式数据转换为SAC格式数据的ref2sac程序；②Matlab平台有很多高级函数，比如：datenum和datestr函数可以很好地解决波形截取过程中遇到的时间运算问题。在处理SAC文件上，也有现成的saclab工具包可以使用，而且在Matlab中可以利用Unix命令方便调用Linux系统的程序。基于Linux系统下的Matlab平台编写程序，可以同时利用Matlab高级函数和Linux系统程序两者的优势，高效实现程序功能。

本程序的功能是从REFTEK连续波形数据中截取事件波形，原始数据为流动观测记录的REFTEK格式的连续波形数据，截取得到SAC格式的三分量事件地震波形数据。目前程序的测试平台为Ubuntu12.04和Matlab2009b，只要确保Matlab相关程序包的安装和Linux系统子程序的成功调用，程序即可顺利运行。程序的运行时间和地震目录的个数及截取波形的长度有关，虽然调用了较多外部程序，但是程序整体运行效率可以接受的。

2 程序运行流程

在调用Matlab程序截取波形前，需要进行相关数据准备。首先是原始连续波形数据，程序要求每个台站记录的REFTEK原始波形数据单独存放在一个文件夹下，而且为了方便访问该文件夹，使用阿拉伯数字给文件夹编号命名。然后是地震目录的整理，地震目录不限行数，每行8列，1到3列分别为年、月、日，第4列为地震发震时刻，格式为时分秒，5到8列分别为震源的纬度、经度、深度和震级大小，各列之间以空格分开。最后是准备台站信息，包括4列，分别为台站编号、台站经度、台站纬度和台站高程。台站信息文件直接控制截取哪些台站的事件波形，如果不需要截取相关台站的波形，则将该台站的信息从台站信息文件中删除。台站信息文件中的台站编号应该和存放该台原始连续波形的文件夹命名编号一致，方便在波形截取中寻找对应台站的数据文件。

程序实际为一个Matlab函数，调用该函数需要输入4个参数，分别为P波（实际为初至波）到达前时间、截取波形的总长度、相对于台站的最小震中距和最大震中距（用来筛选地震目录）。截取事件波形需要一个参考时间，本程序中参考时间是P波（初至波）到达台站的理论到时，函数输入的第1个参数是相对于该参考时间的提前时间，结合截取波形的总时间长度，计算截取波形的时间终点。程序运行流程见图1。

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

（1）从存放台站信息的文件中，读入台站信息，包括台站编号、台站经纬度和高程。然后从第1行开始循环，即波形数据的截取是逐台进行的。

（2）读入地震目录，结合台站经纬度信息和地震目录中地震事件的经纬度信息，可以对地震目录进行挑选，选择标准是调用主程序时输入的最小和最大震中距。如果截取的事件波形用于接收函数研究，一般将震中距限制在30°—95°。

（3）选取地震目录后，计算到时目录。主要计算需要截取波形的前后两个时间点的时间。利用TAUP软件计算P波理论走时，发震时刻加上P波理论走时就得到P波理论到时，即波形截取的参考时间。然后根据调用函数时输入的截取波形提前时间和截取波形总时间长度，计算截取波形前后两个端点的时间。

（4）从第1个地震目录开始，使用arcfetch程序，从REFTEK连续波形中截取已计算两时间点间的事件波形，使用ref2sac程序，将REFTEK格式数据转换为SAC格式数据。

（5）使用saclab程序包处理SAC数据，主要重写SAC数据头文件，需要注意SAC头文件中的时间处理。

（6）利用相关文件处理命令，对SAC文件进行存储，并返回流程（4）截取下一个地震事件，地震目录截取完毕，返回流程（1）截取下一个台站。最后获得的SAC数据仍然以台站为单位进行存放，也就是，同一个台站截取的不同地震事件的SAC文件存放在同一个文件夹。如果想改为一个地震事件文件夹下存放多个台站文件的数据格式，需要编写相应的脚本。

3 关键子程序

3.1 arcfetch

arcfetch程序用来从原始的连续REFTEK波形数据中截取一段波形数据，波形截取的起始时间和波形长度可以在程序输入参数中指定，且截取数据格式仍然为REFTEK数据格式。使用该程序时，需要指定原始连续波形数据的存放路径。为保证程序成功执行，在存放连续波形数据的文件夹中需要一个archive.sta文件，该文件可以通过arccreate程序生成。以上所述两个程序和程序使用的说明文档可以在REFTEK官网下载。arcfetch程序截取连续波形会有3种结果：①截取到完整正确的REFTEK数据文件；②没有截取到REFTEK数据文件；③截取到REFTEK数据文件，但数据为空。对于后两种情况需编写相关测试程序，以避免后续的ref2sac程序出错。

3.2 ref2sac

使用arcfetch截取得到REFTEK事件波形，需要使用ref2sac程序转换为SAC格式。该程序为PASSCAL（大陆岩石圈地震台阵研究计划）提供的rpm软件包中的一个子程序，可以将该子程序包在Linux系统下单独编译使用。该程序有两个参数，分别是输入文件夹和输出文件夹。它将输入文件夹下REFTEK格式数据转换为SAC格式数据，并存放在输出文件夹下。当该程序运行正确时，返回值为0，如果返回值为其他值，表明数据格式转换不成功。

3.3 taup_time

波形截取时，以理论P波到达台站的时间点为截取波形的参考时间，选择截取该参考时间点前后各多少秒的数据。TAUP软件是目前地震学中计算地震波走时的常用软件（Crotwell et al，1999）。本程序利用taup_time计算地震目录中地震事件传播到台站的理论走时，将理论走时加上发震时刻得到理论P波到时，作为波形截取的参考时间。具体实现过程中发现，并不是所有震中距都是P波为初至波，实际上参考时间应该是理论初至波到时。在使用taup_time程序时，只给定参考模型和震中距参数，比时TAUP会计算出所有震相走时，而且初至波走时显示在第1行，本程序提取第1行的震相走时作为实际使用的走时。所以，文中所说P波和初至波是等价的。

3.4 saclab程序包

ref2sac程序只是将REFTEK格式的地震数据转换为SAC格式，但是SAC数据的头文件并不正确，所以在转换得到SAC格式地震数据后，需要进一步修改SAC数据的头文件。saclab是由Andreas Wustefeld使用Matlab语言编写的程序包，便于SAC文件读写和头文件修改。需要注意的是，O、B、A、E这4个SAC头文件，因为涉及到波形数据的时间信息，一旦错误，后续利用该SAC波形数据进行的相关走时研究都将会出现系统性错误。

3.5 datestr 和datenum

Matlab中表示日期的时间有3种格式：日期字符串、连续的日期数值和日期向量（王正林等，2007）。日期字符串格式使用较多，也较为直观，如2013-10-1 10∶10∶10，但是不方便进行时间的加减运算。连续的日期数值格式以公元元年1月1日为起点的，用一个数值表示当前时间到该起点以天为单位的时间距离，该数值为浮点数，所以方便进行时间的加减运算。Matlab中的datenum函数可以实现将日期字符串和日期向量格式转化为连续的日期数字格式，而datestr函数可以实现把日期数字和日期向量格式转换成日期字符串格式，波形截取程序中，需要进行儒略日的计算，利用datenum函数使用一条语句就可以实现， julianday=datenum（year，month，day）-datenum（year，1，1）+1。波形截取程序需要进行时间的加减运算，以确定截取波形的时间起始点，但arcfetch在截取波形时需要输入的时间格式是日期字符串格式，因为日期字符串表示的时间进行加减运算比较困难。本程序采取的方法是，在Matlab中，先将日期字符串格式转化为连续的日期数值格式，然后进行相关的时间运算，再转换回日期字符串的格式，便于时间运算。

4 结束语

REFTEK连续波形截取程序，是一个Matlab函数，功能就是从REFTEK连续波形数据中截取事件波形。该程序连接多个重要的子程序，并使用相应的Matlab工具包，分别完成特定功能，利用Matlab的高端平台特性，以简单编程高效完成波形截取任务。在Linux系统下，利用Matlab平台连接各种处理程序的编程方式，减少编程时间，快速完成特定功能，但程序的移植性和易用性较差，不过作为个人的科研使用仍然具有一定使用价值。利用该程序截取得到的事件波形，可以进行震相到时读取，进而开展走时相关研究工作（黎源等，2012），也可以利用三分量波形数据开展包括接收函数（查小惠等，2013）在内的一系列波形研究工作。目前使用该程序截取的波形所开展的相关研究工作，得到了较好的结果（Lei et al，2012；黎源等，2012； 查小惠等，2013），进一步证明了该程序的可靠性。

程序编写得到中国科学院青藏高原研究所刘启民博士和中国地震局地壳应力研究所孙长青博士的帮助，在此表示感谢。

黎源， 雷建设. 青藏高原东缘上地幔顶部Pn波速度结构及各向异性研究[J]. 地球物理学报， 2012， 55（11）：3 615-3 624.

王正林， 刘明. 精通MATLAB7[M].北京：电子工业出版社， 2007.

查小惠， 雷建设. 云南地区地壳厚度和泊松比研究[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.