基于Matlab的三峡水库地震数据处理与分析1
2014-05-05马文涛
蔺 永 马文涛
(中国地震局地质研究所,北京 100029)
基于Matlab的三峡水库地震数据处理与分析1
蔺 永 马文涛
(中国地震局地质研究所,北京 100029)
本文以Matlab为平台,编制了相应的程序代码,实现了SAC二进制数据读取、地震观测报告数据检验、频谱分析、地震数据的互相关计算和各种地震数据成图等项功能,并将它们运用于长江三峡水库地震加密台网数据的处理。结果表明,三峡库区地震都属于水库诱发地震,在仙女山断裂过江段、九湾溪断裂附近的地震属于断层因蓄水诱发的错动事件;巴东库区神龙溪两岸地震明显呈现出3条线性分布,这是由于水库蓄水后库水从神龙溪两岸等地下暗河渗入而诱发的地震;在长江三峡库区泄滩乡以西的两岸则存在着一些与蓄水相关的塌陷地震。
Matlab 水库地震 SAC数据读取与分析 长江三峡水库
引言
Matlab计算软件是美国MathWorks公司出品的一种商业数学软件,它具有矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言等强大的计算功能,可以用于算法开发、数据可视化、数据分析、数值计算等的高级技术计算语言和交互式环境,用户也可以直接调用Matlab函数库中的各种函数,简单地编制自己的实用程序,广泛地应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。在地震数据分析方面也得到一些运用,牟垒育等(2006)编制了Matlab版的Geiger地震定位软件;万永革(2007)在Matlab上给出了数字信号处理具体应用实例;谭雨文等(2011)利用Matlab解决了台站数据的可视化问题;万永革等(2012)在Geiger地震定位中考虑了数据的不确定性问题。2009年3—12月,由中国地震局地质研究所布设的三峡库区地震台网共记录到了地震事件2995个,其中震相数据45999条。面对如此多的地震数据,单靠人工按照地震事件为单位逐一处理显然是不现实的。因此,本文以三峡水库地震加密台站数据为例,编制了一些基于Matlab平台的处理小程序,开展了Matlab的水库地震数据处理与分析和图形可视化工作,其目的是推进水库地震的研究水平。
1 长江三峡水库地震加密台站数据分析
长江三峡水库是现今世界上发电量最大的发电站,坝高185m,最大库容393亿m3。2003年5月开始蓄水,库区于6月10日开始出现诱发地震,历经145m、156m和175m三期蓄水过程,总共发生水库诱发地震数千次(廖武林等,2009;马文涛等,2010),地震地质灾害明显增多。本文依据国家科技重点课题“典型水库诱发地震危险性评定技术及预警技术研究”(2008BAC38B04),在三峡水库地震台网的基础之上,新建了26个数字地震台站组成的长江三峡水库地震加密台网,其中包括26个美国Reftak130型宽频记录仪、21个英国L-22E/3D型短周期检波器、5个美国GuralpCMG-3ESPC型宽频检波器,能有效控制在三峡库区湖北段库区的水库地震。在2009年3—12月,该台网共记录到发生在三峡库区湖北段的2995次ML0.8—2.9级地震,采样率设计为200Hz,选择minseed格式存储地震数据。使用批处理程序,可以方便地将地震数据转变成SAC二进制数据,便于编制基于Matlab软件平台的后续处理软件,开展地震数据读取、处理与分析。
2 Matlab的数据处理分析与数据成图
在将地震数据调入内存之后,可以利用Matlab的各种模块对数据进行各种分析与成图处理。
2.1 Matlab的SAC二进制数据读取
SAC二进制数据是一种asc码的数据格式,它通常主要包括两个部分:用来说明数据存储方式的固定长度的头段区以及一个或两个数据区。前者一般说明仪器选取的参数、台站位置及时间服务参数等,后者是多通道数据流。具体头段变量说明详见SAC使用手册(Seismic Analysis Code Users Manual)。
要读取SAC数据,在matlab中可以直接调用文件I/O函数:fopen、fclose、fread、fwrite(详细用法参见matlab帮助文档)。
下面给出读取SAC几个常用变量的读取方式:
其中,hd.delta是等间隔数据的采样间隔;hd.b是自变量的起始值;hd.stla和hd.stlo分别是台站的经纬度;hd.evla和hd.evlo分别是事件的经纬度;hd.dist是震中距;hd.az和hd.baz分别是事件到台站的方位角和台站到事件的方位角;hd.npts是数据点数(其他变量信息参见SAC用户手册)。
继续以上的代码读取SAC文件的数据段:
图1 2009年7月1日9点的地震波形图与该地震频谱分析图Fig. 1 Waveform and spectrum of earthquake at 9:00 July 1st, 2009
2.2 利用Matlab对地震观测报告数据进行检验
地震观测报告是基本的地震数据,它包括了每个地震相对于各个台站的P、S波到时、震幅值。但由于波型震相辨别上差异等原因,有一部分数据存在着问题。如何评价这些数据?可以从直达P波、S波的不同走时曲线来解读。利用Matlab文本读取函数textread、textscan、fgetl等对地震观测报告数据中的2995个地震事件的45999条震相记录进行读取,并画出走时图对数据进行挑选(图2中a和b),去掉非P波、非S波等不好的地震数据,及不好的地震数据。挑选后的数据利用Matlab生成hypoDD地震目录数据输入文件(程序代码不再列出)。
图2 走时图对比(a:初始走时图;b:处理数据后走时图)Fig. 2 Comparison of travel time curves (a: initial travel time curve; b: processed travel time curve)
2.3 Matlab对地震波形数据进行频谱分析
频谱分析就是通过快速傅里叶变换实现地震信号从时间域到频率域的转化,从而获得地震事件的优势频率及最大震幅值、拐角频率等参数,进而对其波形进行滤波或者对事件类型进行分析。因此,地震数据的频谱分析在地震类型识别与爆破分析中有着非常广泛的应用。基于Matlab的频谱分析可以对data数组中的数据进行处理:
图1通过频谱分析,可以根据数字滤波器等原理,利用Matlab实现数字滤波,去除干扰信号或干扰源,从波形上初步区分地震与其它振动事件。
2.4 地震数据的Matlab互相关计算
对滤波后的波形数据利用Matlab中的互相关函数xcorr( )进行互相关计算(程序代码不再列出),根据设定相关系数的阈值提取事件对的到时差。最后生成hypoDD所需要的dt.cc文件。结合上面挑选后的地震目录数据就可以利用hypoDD软件进行双差定位了,定位后的成果图如图4—图7所示。
2.5 Matlab的地震数据成图
利用Matlab可以绘制震中分布图、深度分布图、深度频次图、M-T图等。图3是三峡库区地震加密台网的地震数据处理结果的成果图。
图3 双差定位的三峡库区北段2009年3—12月地震分布(为便于结合构造地质条件分析此图由其他绘图软件绘制)Fig. 3 Earthquake distribution of the northern section of Three Gorges Reservoir (data of 2009.3-12)
图3中灰色区为砂板页岩为主的碎屑岩;澄色区是以花岗岩类为主的结晶岩;白色区为以灰岩为主的岩溶;黄色圆圈为震源深度<2km的地震;红色圆圈为震源深度≥2km的地震;斧头叉为煤矿。
图4 双差定位的深度分布图Fig. 4 Depth distribution of hypoDD
图5 深度频次图Fig. 5 The histogram of depth
图6 M-T图Fig. 6 M-T plot
图7 长江三峡水库水位图Fig. 7 Water levelof Three Gorges reservoir with time
对长江三峡水库地区的2995次地震双差定位的结果表明,水平方向平均误差为0.083km,深度方向平均误差为0.107km。地震震中分布明显呈线性集中现象,震源深度主要在6km以上。高精度的定位出小震震群分布图像呈线性分布或团块状丛集分布,长江三峡水库湖北段地震主要集中分布在香溪河附近的仙女山断裂北端及九湾溪断裂、泄滩乡以西的长江两岸和巴东北岸神龙溪及附近地区,震源深度<10km,平均在4km左右。库区地震活动频次与库水位升降过程正相关,说明属于水库诱发地震。
在仙女山断裂过江段、九湾溪断裂和泄滩乡、沙镇溪镇西部地区等的地震可能与仙女山断裂带、牛口断裂或顺层解理等不连续结构面软化,导致岩体失稳而产生水库诱发地震,巴东库区神龙溪两岸地震明显呈现出3条线性分布,通过对比该地区碳酸盐岩的分布特征,发现是由于水库蓄水后,库水从神龙溪两岸等地下暗河渗入而诱发地震,在长江三峡库区泄滩乡以西的两岸存在着一些与蓄水相关的塌陷地震。
3 结论
(1)本文通过在Matlab平台,开展了相关的程序编制,实现了SAC二进制数据读取、地震观测报告数据检验、频谱分析、地震数据的互相关计算和各种地震数据成图等功能,有助于充分利用Matlab平台开展水库地震的研究工作。水库地震研究的基础工作就是要首先确定地震时、空、强的基本要素,本文利用Matlab对三峡库区2995个地震事件的45999条震相记录数据进行了根据走时曲线的数据挑选工作,并将挑选后的数据利用Matlab生成定位程序的输入文件,进而利用Matlab对定位结果进行初步成图。本文在Matlab的数据批量化处理和数据图形可视化上进行了相关的工作,使得水库地震研究的处理效率得到了极大的提高。
(2)作为一种面向科学与工程计算的高级语言,Matlab允许使用数学形式的语言编写程序,这让本文在许多数学计算的代码编写上变得更加简单直观。Matlab还含有丰富的库函数,可直接调用,开展复杂的数学运算,在编制代码过程中,直接调用了如:fft函数、xcorr函数以及文本处理函数等一些库函数,提高了编程效率,充分展现了Matlab其优势所在。
(3)本文利用Matlab对地震观测报告数据利用走时曲线进行了数据挑选和校验工作,并生成了HYPODD程序的输入文件进行定位,定位后又利用Matlab进行成图工作。开展了地质构造条件分析。结果表明,库区地震活动频次与库水位升降过程正相关,都属于水库诱发地震。在仙女山断裂过江段、九湾溪断裂附近的地震属于该断层因蓄水引发的水库诱发地震;巴东库区神龙溪两岸地震明显呈现出3条线性分布,通过对比该地区碳酸盐岩的分布特征,发现是由于水库蓄水后,库水从神龙溪两岸等地下暗河渗入而诱发地震;在长江三峡库区泄滩乡以西的两岸存在着一些与蓄水相关的塌陷地震。
廖武林,张丽芬,姚运生,2009.三峡水库地震活动特征研究.地震地质,31(4):707—714.
马文涛,徐长朋,李海鸥等,2010.长江三峡水库诱发地震加密观测及地震成因初步分析.地震地质,32(4):552—562.
牟垒育,赵仲和,张伟等,2006.用INGLADA和GEIGER方法实现近震精定位.中国地震,22(3):294—302.
谭雨文,刘国明,2011.Matlab在地震信号处理中的应用实例.防灾减灾学报,27(3):61—66.
万永革,2007.数字信号处理的Matlab实现.北京:科学出版社.
万永革,盛书中,程万正等,2012.考虑到时误差的地震定位算法及其在四川地区2001—2008年地震定位的应用.地震地质,34(1):1—10.
Matlab-based Seismic Data Processing and Analysis of Three Gorges Reservoir
Lin Yong and Ma Wentao
(Institute of Geology, China Earthquake Administration, Beijing 100029, China)
Taking Matlab as a platform to read SAC binary data, we tested the seismic report data, performed spectral analysis, calculated cross-correlation and seismic data mapping. We also used our Matlab code to process the seismic data of Three Gorges Reservoir. Our results show that earthquakes around Three Gorges Reservoir are induced by the reservoir. Along Xiannvshan fault and Jiuwanxi fault earthquakes are induced by fault storage. Along Shenlong river in the reservoir area, earthquakes displayed three linear distributions in northern Badong county, which was related to water seeping into the ground through underground rivers. There also exist some collapse induced earthquakes in west Three Gorges Reservoir near Xietan village.
Matlab; Reservoir earthquake; SAC data reading and analysis; Three Gorges Reservoir
蔺永,马文涛,2014.基于Matlab的三峡水库地震数据处理与分析.震灾防御技术,9(3):447—453.
10.11899/zzfy20140311
中国地震局课题(2200408)和国家科技支撑课题(2008BAC38B04)资助
2014-02-12
蔺永,男,生于1984年。中国地震局地质研究所硕士研究生。专业方向:水库诱发地震。E-mail:yong-89 @qq.com。
马文涛,男,生于1958年。中国地震局地质研究所副研究员。主要从事水库诱发地震研究。E-mail:wentaoma @sina.com。