藺 永 馬文濤

（中國(guó)地震局地質(zhì)研究所，北京 100029）

基于Matlab的三峽水庫(kù)地震數(shù)據(jù)處理與分析1

藺 永 馬文濤

（中國(guó)地震局地質(zhì)研究所，北京 100029）

本文以Matlab為平臺(tái)，編制了相應(yīng)的程序代碼，實(shí)現(xiàn)了SAC二進(jìn)制數(shù)據(jù)讀取、地震觀測(cè)報(bào)告數(shù)據(jù)檢驗(yàn)、頻譜分析、地震數(shù)據(jù)的互相關(guān)計(jì)算和各種地震數(shù)據(jù)成圖等項(xiàng)功能，并將它們運(yùn)用于長(zhǎng)江三峽水庫(kù)地震加密臺(tái)網(wǎng)數(shù)據(jù)的處理。結(jié)果表明，三峽庫(kù)區(qū)地震都屬于水庫(kù)誘發(fā)地震，在仙女山斷裂過(guò)江段、九灣溪斷裂附近的地震屬于斷層因蓄水誘發(fā)的錯(cuò)動(dòng)事件；巴東庫(kù)區(qū)神龍溪兩岸地震明顯呈現(xiàn)出3條線性分布，這是由于水庫(kù)蓄水后庫(kù)水從神龍溪兩岸等地下暗河滲入而誘發(fā)的地震；在長(zhǎng)江三峽庫(kù)區(qū)泄灘鄉(xiāng)以西的兩岸則存在著一些與蓄水相關(guān)的塌陷地震。

Matlab 水庫(kù)地震 SAC數(shù)據(jù)讀取與分析 長(zhǎng)江三峽水庫(kù)

引言

Matlab計(jì)算軟件是美國(guó)MathWorks公司出品的一種商業(yè)數(shù)學(xué)軟件，它具有矩陣運(yùn)算、繪制函數(shù)和數(shù)據(jù)、實(shí)現(xiàn)算法、創(chuàng)建用戶界面、連接其他編程語(yǔ)言等強(qiáng)大的計(jì)算功能，可以用于算法開(kāi)發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析、數(shù)值計(jì)算等的高級(jí)技術(shù)計(jì)算語(yǔ)言和交互式環(huán)境，用戶也可以直接調(diào)用Matlab函數(shù)庫(kù)中的各種函數(shù)，簡(jiǎn)單地編制自己的實(shí)用程序，廣泛地應(yīng)用于工程計(jì)算、控制設(shè)計(jì)、信號(hào)處理與通訊、圖像處理、信號(hào)檢測(cè)、金融建模設(shè)計(jì)與分析等領(lǐng)域。在地震數(shù)據(jù)分析方面也得到一些運(yùn)用，牟壘育等（2006）編制了Matlab版的Geiger地震定位軟件；萬(wàn)永革（2007）在Matlab上給出了數(shù)字信號(hào)處理具體應(yīng)用實(shí)例；譚雨文等（2011）利用Matlab解決了臺(tái)站數(shù)據(jù)的可視化問(wèn)題；萬(wàn)永革等（2012）在Geiger地震定位中考慮了數(shù)據(jù)的不確定性問(wèn)題。2009年3—12月，由中國(guó)地震局地質(zhì)研究所布設(shè)的三峽庫(kù)區(qū)地震臺(tái)網(wǎng)共記錄到了地震事件2995個(gè)，其中震相數(shù)據(jù)45999條。面對(duì)如此多的地震數(shù)據(jù)，單靠人工按照地震事件為單位逐一處理顯然是不現(xiàn)實(shí)的。因此，本文以三峽水庫(kù)地震加密臺(tái)站數(shù)據(jù)為例，編制了一些基于Matlab平臺(tái)的處理小程序，開(kāi)展了Matlab的水庫(kù)地震數(shù)據(jù)處理與分析和圖形可視化工作，其目的是推進(jìn)水庫(kù)地震的研究水平。

1 長(zhǎng)江三峽水庫(kù)地震加密臺(tái)站數(shù)據(jù)分析

長(zhǎng)江三峽水庫(kù)是現(xiàn)今世界上發(fā)電量最大的發(fā)電站，壩高185m，最大庫(kù)容393億m3。2003年5月開(kāi)始蓄水，庫(kù)區(qū)于6月10日開(kāi)始出現(xiàn)誘發(fā)地震，歷經(jīng)145m、156m和175m三期蓄水過(guò)程，總共發(fā)生水庫(kù)誘發(fā)地震數(shù)千次（廖武林等，2009；馬文濤等，2010），地震地質(zhì)災(zāi)害明顯增多。本文依據(jù)國(guó)家科技重點(diǎn)課題“典型水庫(kù)誘發(fā)地震危險(xiǎn)性評(píng)定技術(shù)及預(yù)警技術(shù)研究”（2008BAC38B04），在三峽水庫(kù)地震臺(tái)網(wǎng)的基礎(chǔ)之上，新建了26個(gè)數(shù)字地震臺(tái)站組成的長(zhǎng)江三峽水庫(kù)地震加密臺(tái)網(wǎng)，其中包括26個(gè)美國(guó)Reftak130型寬頻記錄儀、21個(gè)英國(guó)L-22E/3D型短周期檢波器、5個(gè)美國(guó)GuralpCMG-3ESPC型寬頻檢波器，能有效控制在三峽庫(kù)區(qū)湖北段庫(kù)區(qū)的水庫(kù)地震。在2009年3—12月，該臺(tái)網(wǎng)共記錄到發(fā)生在三峽庫(kù)區(qū)湖北段的2995次ML0.8—2.9級(jí)地震，采樣率設(shè)計(jì)為200Hz，選擇minseed格式存儲(chǔ)地震數(shù)據(jù)。使用批處理程序，可以方便地將地震數(shù)據(jù)轉(zhuǎn)變成SAC二進(jìn)制數(shù)據(jù)，便于編制基于Matlab軟件平臺(tái)的后續(xù)處理軟件，開(kāi)展地震數(shù)據(jù)讀取、處理與分析。

2 Matlab的數(shù)據(jù)處理分析與數(shù)據(jù)成圖

在將地震數(shù)據(jù)調(diào)入內(nèi)存之后，可以利用Matlab的各種模塊對(duì)數(shù)據(jù)進(jìn)行各種分析與成圖處理。

2.1 Matlab的SAC二進(jìn)制數(shù)據(jù)讀取

SAC二進(jìn)制數(shù)據(jù)是一種asc碼的數(shù)據(jù)格式，它通常主要包括兩個(gè)部分：用來(lái)說(shuō)明數(shù)據(jù)存儲(chǔ)方式的固定長(zhǎng)度的頭段區(qū)以及一個(gè)或兩個(gè)數(shù)據(jù)區(qū)。前者一般說(shuō)明儀器選取的參數(shù)、臺(tái)站位置及時(shí)間服務(wù)參數(shù)等，后者是多通道數(shù)據(jù)流。具體頭段變量說(shuō)明詳見(jiàn)SAC使用手冊(cè)（Seismic Analysis Code Users Manual）。

要讀取SAC數(shù)據(jù)，在matlab中可以直接調(diào)用文件I/O函數(shù)：fopen、fclose、fread、fwrite（詳細(xì)用法參見(jiàn)matlab幫助文檔）。

下面給出讀取SAC幾個(gè)常用變量的讀取方式：

其中，hd.delta是等間隔數(shù)據(jù)的采樣間隔；hd.b是自變量的起始值；hd.stla和hd.stlo分別是臺(tái)站的經(jīng)緯度；hd.evla和hd.evlo分別是事件的經(jīng)緯度；hd.dist是震中距；hd.az和hd.baz分別是事件到臺(tái)站的方位角和臺(tái)站到事件的方位角；hd.npts是數(shù)據(jù)點(diǎn)數(shù)（其他變量信息參見(jiàn)SAC用戶手冊(cè)）。

繼續(xù)以上的代碼讀取SAC文件的數(shù)據(jù)段：

圖1 2009年7月1日9點(diǎn)的地震波形圖與該地震頻譜分析圖Fig. 1 Waveform and spectrum of earthquake at 9:00 July 1st, 2009

2.2 利用Matlab對(duì)地震觀測(cè)報(bào)告數(shù)據(jù)進(jìn)行檢驗(yàn)

地震觀測(cè)報(bào)告是基本的地震數(shù)據(jù)，它包括了每個(gè)地震相對(duì)于各個(gè)臺(tái)站的P、S波到時(shí)、震幅值。但由于波型震相辨別上差異等原因，有一部分?jǐn)?shù)據(jù)存在著問(wèn)題。如何評(píng)價(jià)這些數(shù)據(jù)？可以從直達(dá)P波、S波的不同走時(shí)曲線來(lái)解讀。利用Matlab文本讀取函數(shù)textread、textscan、fgetl等對(duì)地震觀測(cè)報(bào)告數(shù)據(jù)中的2995個(gè)地震事件的45999條震相記錄進(jìn)行讀取，并畫(huà)出走時(shí)圖對(duì)數(shù)據(jù)進(jìn)行挑選（圖2中a和b），去掉非P波、非S波等不好的地震數(shù)據(jù)，及不好的地震數(shù)據(jù)。挑選后的數(shù)據(jù)利用Matlab生成hypoDD地震目錄數(shù)據(jù)輸入文件（程序代碼不再列出）。

圖2 走時(shí)圖對(duì)比（a：初始走時(shí)圖；b：處理數(shù)據(jù)后走時(shí)圖）Fig. 2 Comparison of travel time curves (a: initial travel time curve; b: processed travel time curve)

2.3 Matlab對(duì)地震波形數(shù)據(jù)進(jìn)行頻譜分析

頻譜分析就是通過(guò)快速傅里葉變換實(shí)現(xiàn)地震信號(hào)從時(shí)間域到頻率域的轉(zhuǎn)化，從而獲得地震事件的優(yōu)勢(shì)頻率及最大震幅值、拐角頻率等參數(shù)，進(jìn)而對(duì)其波形進(jìn)行濾波或者對(duì)事件類型進(jìn)行分析。因此，地震數(shù)據(jù)的頻譜分析在地震類型識(shí)別與爆破分析中有著非常廣泛的應(yīng)用。基于Matlab的頻譜分析可以對(duì)data數(shù)組中的數(shù)據(jù)進(jìn)行處理：

圖1通過(guò)頻譜分析，可以根據(jù)數(shù)字濾波器等原理，利用Matlab實(shí)現(xiàn)數(shù)字濾波，去除干擾信號(hào)或干擾源，從波形上初步區(qū)分地震與其它振動(dòng)事件。

2.4 地震數(shù)據(jù)的Matlab互相關(guān)計(jì)算

對(duì)濾波后的波形數(shù)據(jù)利用Matlab中的互相關(guān)函數(shù)xcorr( )進(jìn)行互相關(guān)計(jì)算（程序代碼不再列出），根據(jù)設(shè)定相關(guān)系數(shù)的閾值提取事件對(duì)的到時(shí)差。最后生成hypoDD所需要的dt.cc文件。結(jié)合上面挑選后的地震目錄數(shù)據(jù)就可以利用hypoDD軟件進(jìn)行雙差定位了，定位后的成果圖如圖4—圖7所示。

2.5 Matlab的地震數(shù)據(jù)成圖

利用Matlab可以繪制震中分布圖、深度分布圖、深度頻次圖、M-T圖等。圖3是三峽庫(kù)區(qū)地震加密臺(tái)網(wǎng)的地震數(shù)據(jù)處理結(jié)果的成果圖。

圖3 雙差定位的三峽庫(kù)區(qū)北段2009年3—12月地震分布（為便于結(jié)合構(gòu)造地質(zhì)條件分析此圖由其他繪圖軟件繪制）Fig. 3 Earthquake distribution of the northern section of Three Gorges Reservoir (data of 2009.3-12)

圖3中灰色區(qū)為砂板頁(yè)巖為主的碎屑巖；澄色區(qū)是以花崗巖類為主的結(jié)晶巖；白色區(qū)為以灰?guī)r為主的巖溶；黃色圓圈為震源深度＜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 長(zhǎng)江三峽水庫(kù)水位圖Fig. 7 Water levelof Three Gorges reservoir with time

對(duì)長(zhǎng)江三峽水庫(kù)地區(qū)的2995次地震雙差定位的結(jié)果表明，水平方向平均誤差為0.083km，深度方向平均誤差為0.107km。地震震中分布明顯呈線性集中現(xiàn)象，震源深度主要在6km以上。高精度的定位出小震震群分布圖像呈線性分布或團(tuán)塊狀叢集分布，長(zhǎng)江三峽水庫(kù)湖北段地震主要集中分布在香溪河附近的仙女山斷裂北端及九灣溪斷裂、泄灘鄉(xiāng)以西的長(zhǎng)江兩岸和巴東北岸神龍溪及附近地區(qū)，震源深度＜10km，平均在4km左右。庫(kù)區(qū)地震活動(dòng)頻次與庫(kù)水位升降過(guò)程正相關(guān)，說(shuō)明屬于水庫(kù)誘發(fā)地震。

在仙女山斷裂過(guò)江段、九灣溪斷裂和泄灘鄉(xiāng)、沙鎮(zhèn)溪鎮(zhèn)西部地區(qū)等的地震可能與仙女山斷裂帶、牛口斷裂或順層解理等不連續(xù)結(jié)構(gòu)面軟化，導(dǎo)致巖體失穩(wěn)而產(chǎn)生水庫(kù)誘發(fā)地震，巴東庫(kù)區(qū)神龍溪兩岸地震明顯呈現(xiàn)出3條線性分布，通過(guò)對(duì)比該地區(qū)碳酸鹽巖的分布特征，發(fā)現(xiàn)是由于水庫(kù)蓄水后，庫(kù)水從神龍溪兩岸等地下暗河滲入而誘發(fā)地震，在長(zhǎng)江三峽庫(kù)區(qū)泄灘鄉(xiāng)以西的兩岸存在著一些與蓄水相關(guān)的塌陷地震。

3 結(jié)論

（1）本文通過(guò)在Matlab平臺(tái)，開(kāi)展了相關(guān)的程序編制，實(shí)現(xiàn)了SAC二進(jìn)制數(shù)據(jù)讀取、地震觀測(cè)報(bào)告數(shù)據(jù)檢驗(yàn)、頻譜分析、地震數(shù)據(jù)的互相關(guān)計(jì)算和各種地震數(shù)據(jù)成圖等功能，有助于充分利用Matlab平臺(tái)開(kāi)展水庫(kù)地震的研究工作。水庫(kù)地震研究的基礎(chǔ)工作就是要首先確定地震時(shí)、空、強(qiáng)的基本要素，本文利用Matlab對(duì)三峽庫(kù)區(qū)2995個(gè)地震事件的45999條震相記錄數(shù)據(jù)進(jìn)行了根據(jù)走時(shí)曲線的數(shù)據(jù)挑選工作，并將挑選后的數(shù)據(jù)利用Matlab生成定位程序的輸入文件，進(jìn)而利用Matlab對(duì)定位結(jié)果進(jìn)行初步成圖。本文在Matlab的數(shù)據(jù)批量化處理和數(shù)據(jù)圖形可視化上進(jìn)行了相關(guān)的工作，使得水庫(kù)地震研究的處理效率得到了極大的提高。

（2）作為一種面向科學(xué)與工程計(jì)算的高級(jí)語(yǔ)言，Matlab允許使用數(shù)學(xué)形式的語(yǔ)言編寫(xiě)程序，這讓本文在許多數(shù)學(xué)計(jì)算的代碼編寫(xiě)上變得更加簡(jiǎn)單直觀。Matlab還含有豐富的庫(kù)函數(shù)，可直接調(diào)用，開(kāi)展復(fù)雜的數(shù)學(xué)運(yùn)算，在編制代碼過(guò)程中，直接調(diào)用了如：fft函數(shù)、xcorr函數(shù)以及文本處理函數(shù)等一些庫(kù)函數(shù)，提高了編程效率，充分展現(xiàn)了Matlab其優(yōu)勢(shì)所在。

（3）本文利用Matlab對(duì)地震觀測(cè)報(bào)告數(shù)據(jù)利用走時(shí)曲線進(jìn)行了數(shù)據(jù)挑選和校驗(yàn)工作，并生成了HYPODD程序的輸入文件進(jìn)行定位，定位后又利用Matlab進(jìn)行成圖工作。開(kāi)展了地質(zhì)構(gòu)造條件分析。結(jié)果表明，庫(kù)區(qū)地震活動(dòng)頻次與庫(kù)水位升降過(guò)程正相關(guān)，都屬于水庫(kù)誘發(fā)地震。在仙女山斷裂過(guò)江段、九灣溪斷裂附近的地震屬于該斷層因蓄水引發(fā)的水庫(kù)誘發(fā)地震；巴東庫(kù)區(qū)神龍溪兩岸地震明顯呈現(xiàn)出3條線性分布，通過(guò)對(duì)比該地區(qū)碳酸鹽巖的分布特征，發(fā)現(xiàn)是由于水庫(kù)蓄水后，庫(kù)水從神龍溪兩岸等地下暗河滲入而誘發(fā)地震；在長(zhǎng)江三峽庫(kù)區(qū)泄灘鄉(xiāng)以西的兩岸存在著一些與蓄水相關(guān)的塌陷地震。

廖武林，張麗芬，姚運(yùn)生，2009．三峽水庫(kù)地震活動(dòng)特征研究．地震地質(zhì)，31（4）：707—714．

馬文濤，徐長(zhǎng)朋，李海鷗等，2010．長(zhǎng)江三峽水庫(kù)誘發(fā)地震加密觀測(cè)及地震成因初步分析．地震地質(zhì)，32（4）：552—562．

牟壘育，趙仲和，張偉等，2006．用INGLADA和GEIGER方法實(shí)現(xiàn)近震精定位．中國(guó)地震，22（3）：294—302.

譚雨文，劉國(guó)明，2011．Matlab在地震信號(hào)處理中的應(yīng)用實(shí)例．防災(zāi)減災(zāi)學(xué)報(bào)，27（3）：61—66.

萬(wàn)永革，2007．?dāng)?shù)字信號(hào)處理的Matlab實(shí)現(xiàn)．北京：科學(xué)出版社.

萬(wàn)永革，盛書(shū)中，程萬(wàn)正等，2012．考慮到時(shí)誤差的地震定位算法及其在四川地區(qū)2001—2008年地震定位的應(yīng)用．地震地質(zhì)，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的三峽水庫(kù)地震數(shù)據(jù)處理與分析．震災(zāi)防御技術(shù)，9（3）：447—453．

10.11899/zzfy20140311

中國(guó)地震局課題（2200408）和國(guó)家科技支撐課題（2008BAC38B04）資助

2014-02-12

藺永，男，生于1984年。中國(guó)地震局地質(zhì)研究所碩士研究生。專業(yè)方向：水庫(kù)誘發(fā)地震。E-mail：yong-89 @qq.com。

馬文濤，男，生于1958年。中國(guó)地震局地質(zhì)研究所副研究員。主要從事水庫(kù)誘發(fā)地震研究。E-mail：wentaoma @sina.com。