基于IASP91地球模型的天然地震人机交互数据处理系统

徐金荣， 王有学， 熊 彬， 蒋婵君， 乃振龙， 张 琪， 王心宇， 曾 成， 胡锦锋

(桂林理工大学 地球科学学院，桂林 541006)

基于IASP91地球模型的天然地震人机交互数据处理系统

徐金荣， 王有学， 熊 彬， 蒋婵君， 乃振龙， 张 琪， 王心宇， 曾 成， 胡锦锋

(桂林理工大学 地球科学学院，桂林 541006)

提出了一种基于IASP91地球模型确定各种地震波震相及其旅行时间拾取的人机交互方法，通过Intel Fortran平台设计并实现了天然地震记录可视化数据处理及多震相旅行时间拾取的交互式处理系统。该处理系统可以快速高效地处理天然地震数据、确定及拾取天然地震各震相的旅行时间。实际应用表明，该系统操作简单，准确快速，灵活方便，是处理天然地震数据的有效工具。

天然地震； 人机交互； 数据处理

0 引言

地震学的任务之一是通过对天然地震的观测及地震数据的分析，获得各个地震事件所产生的不同地震波震相的传播特征及旅行时间，进而研究地球内部的结构特征[1]。然而对于一个流动观测台阵的观测数据，其观测数据量是惊人的。在海量的天然地震观测数据中寻找出天然地震信号，并且确定该次事件地震记录中的各个震相，是一件单一、重复性的且非常耗时的工作。同时，由于不同宽频带地震仪的观测数据格式各不相同，诸如SAC、Miniseed、CSS、DMX等，使得天然地震观测数据的处理与分析变得更为复杂。此时高效的天然地震数据处理系统的开发就变得非常必要。

震相地拾取是处理天然地震数据较为重要的一个环节，准确的震相拾取是进行地震波速度成像等地震数据处理方法的前提[2-3]。目前，在处理天然地震记录的过程中，地震波震相的拾取主要有2种方式：①手动拾取；②计算机自动拾取。

手动拾取主要有：①直观检验法；②综合检验法；③达曲线检验法；④合成地震图检验法等[4]，后来还有对地震记录做一些放大、仿真、滤波之类的简单处理的半自动识别方法。这些方法都是靠识别者的经验来判断、识别，不仅费时，而且会产生较大的人为读数误差；另一方面，即使对于同一个地震事件记录，不同的处理人员处理结果也会有所差别。

计算机自动拾取震相的方法主要有：①时域分析法(STA/LTA法、能量比法、振幅比法等)；②频域分析法(Fourier变换、Walsh变换、最大熵谱等)；③模式识别法等[5-8]，其主要特点就是工作效率相对较高。但是它的各种误判也不容忽视，比如在地震记录中存在较强的干扰，采用自动拾取方法会产生误判，难以获取准确的震相波至。因此，在地震记录数据量非常大的情况下，只靠人工手动拾取不仅效率低，而且是不切合实际的；自动拾取方法虽然工作效率高，但对于存在较强噪音背景的地震事件的震相会产生误判。因此，研究和开发一个稳定、实用和高效的人机交互式天然地震数据处理系统,具有一定的现实意义和应用价值。

笔者使用Intel Fortran平台进行程序语言编译[9]，基于windows系统搭建人机交互界面[10]，以便更简洁、高效地实现天然地震数据的可视化显示，更直观通畅地实现人机交互[11-12]。同样，通过人机交互界面实现地震震相的准确识别及拾取等功能。

1 人机交互系统的总体设计

为了提高宽频带地震数据的处理速度及工作效率，笔者拟设计一套基于IASP91地球模型的天然地震人机交互数据处理系统，该系统采用人机交互方式及数据可视化技术，不仅能够对宽频带地震数据进行数据可视化、预处理、震相走时计算及拾取、数据的旋转及重采样等，而且能够将处理结果进行保存、输出，以便进一步应用于地球内部结构的研究。天然地震人机交互数据处理系统的总体设计如图1所示。

图1 天然地震人机交互数据处理系统总设计图Fig.1 Site-plan of the interactive processing system

2 数据可视化及预处理

要实现天然地震数据的震相拾取及特征分析

等，首先需要将地震数据导入计算机内存，然后根据其波形特征进行人机交互的数据处理，完成数据分析；同时，天然地震数据记录中可能存在坏道及严重的干扰等，需要对实际观测的地震数据进行滤波、坏道剔除等预处理。

2.1 地震数据的可视化

地震数据的高效可视化，可以提高地震数据处理工作的效率。笔者介绍的数据处理系统可以读取宽频地震仪记录的多种格式(如SAC、DMX、CSS)等，同时能够对该数据进行二维可视化。在进行地震数据二维可视化过程中，首先要对地震数据进行归一化处理(道归一)[13]，以避免各记录地相互干扰，然后将地震数据进行显示。数据可视化的形式包括：①波形图；②变面积图。

1)波形图是用曲线的形式表示一系列离散的地震数据。仪器采集的地震数据是一组离散的值，在屏幕坐标中是一个个的离散点。在进行地震波的波形图绘制前常需要根据屏幕分辨率和采样点的个数关系对地震数据进行处理。如果采样点较少就会出现折线，因此，在进行地震曲线的绘制时需要对数据进行插值使其圆滑。常用的插值方法有：①拉格朗日插值；②三次样条插值；③二维插值等[15].文中的数据处理系统利用三次样条算法进行插值，在算法中只需要提供分段样条曲线的控制点坐标值，由节点内部封闭的算法来完成样条曲线的绘制。

例如已知控制点p1、p2、p3、p4，利用三次样条函数，p1为起始点，p2 和p3为对应控制点，p4为结束点。经过三次样条插值，既能使数据曲线圆滑，又可实现较小的插值误差(图2(a)和图2(b))。

2)变面积图是将地震波形数据的正值或负值区域用蓝色笔刷进行填充。在进行地震数据显示时，为了让波形及同相轴更加清晰，易于分辨，常将地震数据的正值或者负值区域进行填充(图2(c))。

2.2 地震数据的预处理

一方面由于故障可能引起的地震仪器的工作异常，地震记录形成坏道；另一方面，较强的干扰也会影响地震信号的信噪比，影响地震信号的识别与分析。这就需要对数据进行必要的预处理，包括坏道剔除、滤波、增益控制等。

图2 数据的可视化Fig.2 Visudization of seismic data(a)连接p1、p2、p3、p4的折线图；(b)连接p1、p2、p3、p4三次样条图；(c)波形的变面积显示

1)坏道剔除。坏道的存在不仅影响处理结果的准确度，同时也使得系统的处理速度变慢。文中的数据处理系统使用人机交互的方式进行坏道剔除，在地震波形图中鼠标左键点击进行坏道选择，再执行系统中的坏道剔除动作，即可完成坏道剔除；如果用户发现错选了良好的坏道，也可交互取消选中的道集，重新进行坏道的拾取、剔除。值得注意的是，在选择坏道过程中，即使每个台站选择任意一个分量，在剔除坏道过程中，不管屏幕显示的是哪个分量，被选择的坏道所在台站的所有数据将从该次读入中移除，不再参加处理。

2) 滤波处理。在滤波处理过程中，本文数据处理系统采用带通滤波器(图3)，并且利用快速傅里叶变换(FFT)对数据进行滤波处理，以提高地震记录的信噪比[16]。带通滤波器的频率响应H(f)的数学表达式为式(1)。

H(f)=

(1)

其中：fl为低频截止频率；fh为高频截止频率；f1=fl/2，f2=fh+(fh-fl)/2。

图3 带通滤波器示意图Fig.3 Schematic of band-pass filter

3)增益处理。为了使有效的弱信号更为明显，使强弱不一的有效信号均能在屏幕上均匀地显示，常需要对地震信号进行增益控制。增益是对地震记录的一种时变比例均衡，这种比例函数是根据所期望的规则确定的。如果地震数据的时间序列为f(t)，增益函数为c(t)，那么，经增益处理后的地震记录F(t)为式(2)。

(2)

本系统所使用的增益函数c(t)为一常数，并且用户可以通过人机交互界面的地震记录的干扰情况进行分析，通过设置项选择合适的值进行设定。

3 系统功能

3.1 震相的拾取

天然地震信号中常见的地震震相有P、S、PKP、PKS、SKP、SKS、PKIKP、PKIKS、SKIKP、SKIKS以及LR等。基于IASP91模型[14-15]，在已知震源点及接收台站的情况下，利用打靶法这一传统的射线追踪方法[16]，即可获得各震相的理论旅行时间，为震相人机交互识别和拾取提供参考。基于此思想，我们的处理系统可快速计算某一地震对应震相到达各台站的时间，同时软件对该时间点进行标识并显示该标识点于屏幕上，这也是地震记录重新排列等功能的基础。IASP91模型及地球内部的地震震相如表1及图4所示。

图4 IASP91模型中的各种震相示意图Fig.4 Phases in the IASP91 model

3.2 地震数据的旋转

地震记录有东西(EW)、南北(NS)、垂直(VZ)三个分量，在地震数据的研究工作(如接收函数、横波分裂等)中，需要将地震波的两个水平分量分解为径向(SV波)及横向(SH波)两个分量。在地震事件已经确定的前提下，利用台站相对于震中的方位角β，可以将实际观测的南北分量N(t)、东西分量E(t)进行旋转，得到径向分量R(t)及横向分量T(t)为式(3)。

表1 IASP91地球速度模型参数表

注：表中x=r/R，其中r为距地心的半径，R为地球的半径，且R=6 371 km

(3)

3.3 地震数据的合并及存储

天然地震数据以台站为单位存储观测的特点，使得数据的处理及分析，尤其是利用台站间的互相关进行地震波的震相确定造成很大的不便。故在设计该地震数据处理系统时，使得经过系统的预处理、地震事件确定、震相识别等，将系统窗口内显示的全部数据合并存储于DAT格式的文件中，并以该事件发生的年月日+时分秒+设定震相+扩展名(.dat)进行命名。该功能使得后期对该地震记录的处理大大提高了效率。

4 系统的人机交互设计

系统的人机交互贯彻了数据输入到数据处理完成的全过程(图5)。

图5 震相拾取过程的人机交互示意图Fig.5 The interaction of phases picking

经过了该系统对地震数据的导入、预处理，接着进行人机交互的震相识别及拾取。在进行各地震震相的识别和拾取时，需要响应动作①：系统首先从准备好的地震事件列表文件中，自动筛选出地震数据所在日期的所有地震事件，人工选择地震记录所对应的地震事件；然后手动输入需要识别的震相。这时，系统执行对应震相的识别计算，通过IASP91模型的速度参数，结合射线追踪理论，便可确认对应震相在每个独立台站的到时。接着响应动作②：机器拾取的震相需要人工判断其震相识别正误。如果识别正确，则进行下一步的震相人机交互精确拾取，保存拾取结果，结束交互动作；如果识别错误，则返回动作①，从今选择地震事件，直至系统识别到正确的震相。

5 应用实例

选用南岭地区天然地震流动台站的观测数据的一部分(图6)，以日本福岛外海2013年10月26日(北京东八区时间)7.1级地震为例，检验该处理系统的震相识别及拾取效果。

5.1 参数设定

在系统的参数设定界面Paramerters Setting 中，对所需要读取的数据段进行设置(图7)。其中包括数据段开始时刻、时段长度、滤波系数、分量显示、波形放大、变面积图显示、记录排序等设置项。由于福岛地震距离台网最近台站约为3 000 km，到达台站的时刻大约在地震发生后的1 100 s附近，故在示例中参数设定数据段长度设为30 min,起始时间当延后与发震时刻，保证初至波及后续波都被读取进来。完成了上述的参数设置后，通过菜单栏的file->open，并且选择相应的数据格式文件，即可导入与之对应时段的各台站地震数据，并在屏幕上进行数据的可视化。图8为进过预处理后的地震数据显示。

图6 南岭地区天然地震流动台阵Fig.6 Portable seismic array in Nanling

5.2 处理效果

完成了数据的预处理后，通过参数设置窗口，并且完成响应动作①及响应动作②后(图4)，即可进一步进去的震相精确拾取界面(图9)。

图9可以看出，经过地震记录的震相识别、波形放大、变面积图显示，可以清楚地确定该地震的P震相的同相轴，红色十字线处的同相轴提供了清晰的拾取标准。为了保证拾取的准确性，用户可进行多次的拾取选定，且以对同一台站的最后一次拾取为准。这样既可快速处理海量的地震数据，又保证了震相到时人工拾取的精确度。

图7 参数设置及事件选择Fig.7 Parameters setting and seismic event selection

图8 预处理后的地震记录Fig.8 Seismogram after bad traces elimination and filtering

图9 人机交互拾取P震相示意图Fig.9 Interactive picking of Seismic phases

6 结论

根据天然地震数据处理的需求，笔者基于Intel Fortran 开发平台，在Windows环境下设计并实现了一套天然地震人机交互式的数据处理系统，该系统能够实现多种格式的地震数据导入、预处理、震相识别及到时拾取等功能。实际测试结果表明，该人机交互系统操作方便，效果良好，为天然地震数据的处理及分析提供了一个高效、灵活、方便的工具。

[1] 傅淑芳,刘宝成,李文艺.地震学教程[M].北京：地震出版社,1984. FU S F, LIU B C, LI W Y, et al. Seismology course book [M]. Beijing: Seismological Press, 1984.(In Chinese)

[2] AKI K, LEE W H. Determination of the three-dimensional seismic structure of the lithosphere [J]. Geophys. Res.,1977,82(1):277-296.

[3] RAWLISON N, SAMBRIDGE M S. Seismic traveltime tomography of the crust and lithosphere [J]. Advances in Geophysics, 2003, 46: 81-198.

[4] 王继,陈九辉.流动地震台阵观测初至震相的自动检测[J].地震学报,2006,28(1)：42-51. WANG J, CHEN J H. Automatic onset phase picking for portable seismic array observation [J]. Acta Seismological sinica, 2006, 28(1): 42-51. (In Chinese)

[5] STEVENSON P. R.. Micro-earthquakes at Flathead Lake, Montana: A study using automatic earthquake processing [J]. Bull. Seism. Soc. Am., 1976, 66(1): 61-80.

[6] 高淑芳,李山有.一种改进的STA/LTA震相自动识别方法[J].世界地震工程,2008,24(2)：37-41. GAO S F, LI S Y. A modified STA/LTA method to identified earthquake phase automatically [J]. World Earthquake Engineering, 2008, 24(2):37-41. (In Chinese)

[7] COPPENS F. First arrival picking on common-offset trace collections for automatic estimation of static corrections [J]. Geophysical Prospecting, 1985, 33(8): 1212-1231.

[8] TOM GOFORTH, EUGENE HERRIN. An automatic seismic signal detection algorithm based on the Walsh transform [J]. Bull. seism.Soc.Am., 1981, 71(4):1351-1360.

[9] 刘卫国,蔡旭辉.Fortran 90程序设计教程[M].北京：北京邮电大学出版社,2007. LIU W G, CAI X H. Programming tutorial of Fortran 90 [M]. Beijing: Beijing university of posts and telecommunications publishing house, 2007. (In Chinese)

[10]张宏林,孔艳,王哲．按实例学 Visual Basic 6.0[M]．北京：人民邮电出版社，2000. ZHANG H L, KONG Y, WANG Z. Learning Visual Basic 6.0 based on examples [M]. Beijing: Posts and Telecom Press, 2000. (In Chinese)

[11]李乐山.人机界面设计[M].北京：科学出版社.2004. LI L S. Interface Design [M]. Beijing: Science Press. 2004. (In Chinese)

[12]RASMUSSEN J. Information Processing and Human-Machine Interaction: An Approach to Cognitive Engineering [M].New York: Elsevier, 1986.

[13]葛文庚.数据可视化及其在地学中的应用[J].工程地球物理学报,2006(6):465-469. GE WG. Data visualization and its application in the Seismology [J]. Journal of engineering geophysics,2006(6):465-469. (In Chinese)

[14]CERVENY V. Seismic Ray Theory [M]. Cambridge: Cambridge University Press, 2001.

[15]王有学.地震波理论基础[M].北京：地质出版社,2011. WANG Y X. Introduction of seismic wave theory [M]. Beijing: Geology Publishing House, 2011. (In Chinese)

[16]JULIAN B R, GUBBINS D. Three dimensional seismic ray tracing[J].J Geophys, 1977, 43: 95-119.

An interactive system based on the IASP91 model for earthquake data processing

XU Jinrong, WANG Youxue, XIONG Bin, JIANG Chanjun, NAI Zhenlong, ZHANG Qi, WANG Xinyu, ZENG Cheng, HU Jinfeng

(School of Earth Sciences, Guilin University of Technology, Guilin 541006, China)

This paper puts forward a new interactive processing system which based on the IASP91 model and can determine earthquake's phases and picks visually the travel times. This system is designed on the Window platform by using the Intel Fortran compiler to realize the visual seismic data processing, phase identifications and traveltime picking, so it can process seismic data, identify phases and pick seismic arrivals efficiently. The manipulation shows that the friendly system is an effective tools for seismic data process quickly, accurately and easily.

earthquake; interactive; data processing

2016-11-19 改回日期：2017-03-07

国家自然科学基金项目(41574030,41174077，41674075)；广西自然科学基金重点基金项目(2016GXNSFDA380014);大陆构造与动力学国家重点实验室开放课题(k201505)

徐金荣(1989-)，男，硕士，研究方向为固体地球物理，E-mail：10360928@qq.com。

王有学(1961-)，男，博士，教授，主要研究方向为地球物理、地球内部结构及动力学，E-mail：uxue.wang@glut.edu.cn。

1001-1749(2017)02-0275-07

P631.4

A

10.3969/j.issn.1001-1749.2017.02.19