藺 永 馬文濤
(中國地震局地質研究所,北京 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的水庫地震數據處理與分析和圖形可視化工作,其目的是推進水庫地震的研究水平。
長江三峽水庫是現今世界上發電量最大的發電站,壩高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軟件平臺的后續處理軟件,開展地震數據讀取、處理與分析。
在將地震數據調入內存之后,可以利用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對地震波形數據進行頻譜分析
頻譜分析就是通過快速傅里葉變換實現地震信號從時間域到頻率域的轉化,從而獲得地震事件的優勢頻率及最大震幅值、拐角頻率等參數,進而對其波形進行濾波或者對事件類型進行分析。因此,地震數據的頻譜分析在地震類型識別與爆破分析中有著非常廣泛的應用?;贛atlab的頻譜分析可以對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條線性分布,通過對比該地區碳酸鹽巖的分布特征,發現是由于水庫蓄水后,庫水從神龍溪兩岸等地下暗河滲入而誘發地震,在長江三峽庫區泄灘鄉以西的兩岸存在著一些與蓄水相關的塌陷地震。
(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。