徐金榮， 王有學(xué)， 熊 彬， 蔣嬋君， 乃振龍， 張 琪， 王心宇， 曾 成， 胡錦鋒

(桂林理工大學(xué) 地球科學(xué)學(xué)院，桂林 541006)

基于IASP91地球模型的天然地震人機(jī)交互數(shù)據(jù)處理系統(tǒng)

徐金榮， 王有學(xué)， 熊 彬， 蔣嬋君， 乃振龍， 張 琪， 王心宇， 曾 成， 胡錦鋒

(桂林理工大學(xué) 地球科學(xué)學(xué)院，桂林 541006)

提出了一種基于IASP91地球模型確定各種地震波震相及其旅行時間拾取的人機(jī)交互方法，通過Intel Fortran平臺設(shè)計并實(shí)現(xiàn)了天然地震記錄可視化數(shù)據(jù)處理及多震相旅行時間拾取的交互式處理系統(tǒng)。該處理系統(tǒng)可以快速高效地處理天然地震數(shù)據(jù)、確定及拾取天然地震各震相的旅行時間。實(shí)際應(yīng)用表明，該系統(tǒng)操作簡單，準(zhǔn)確快速，靈活方便，是處理天然地震數(shù)據(jù)的有效工具。

天然地震； 人機(jī)交互； 數(shù)據(jù)處理

0 引言

地震學(xué)的任務(wù)之一是通過對天然地震的觀測及地震數(shù)據(jù)的分析，獲得各個地震事件所產(chǎn)生的不同地震波震相的傳播特征及旅行時間，進(jìn)而研究地球內(nèi)部的結(jié)構(gòu)特征[1]。然而對于一個流動觀測臺陣的觀測數(shù)據(jù)，其觀測數(shù)據(jù)量是驚人的。在海量的天然地震觀測數(shù)據(jù)中尋找出天然地震信號，并且確定該次事件地震記錄中的各個震相，是一件單一、重復(fù)性的且非常耗時的工作。同時，由于不同寬頻帶地震儀的觀測數(shù)據(jù)格式各不相同，諸如SAC、Miniseed、CSS、DMX等，使得天然地震觀測數(shù)據(jù)的處理與分析變得更為復(fù)雜。此時高效的天然地震數(shù)據(jù)處理系統(tǒng)的開發(fā)就變得非常必要。

震相地拾取是處理天然地震數(shù)據(jù)較為重要的一個環(huán)節(jié)，準(zhǔn)確的震相拾取是進(jìn)行地震波速度成像等地震數(shù)據(jù)處理方法的前提[2-3]。目前，在處理天然地震記錄的過程中，地震波震相的拾取主要有2種方式：①手動拾??；②計算機(jī)自動拾取。

手動拾取主要有：①直觀檢驗(yàn)法；②綜合檢驗(yàn)法；③達(dá)曲線檢驗(yàn)法；④合成地震圖檢驗(yàn)法等[4]，后來還有對地震記錄做一些放大、仿真、濾波之類的簡單處理的半自動識別方法。這些方法都是靠識別者的經(jīng)驗(yàn)來判斷、識別，不僅費(fèi)時，而且會產(chǎn)生較大的人為讀數(shù)誤差；另一方面，即使對于同一個地震事件記錄，不同的處理人員處理結(jié)果也會有所差別。

計算機(jī)自動拾取震相的方法主要有：①時域分析法(STA/LTA法、能量比法、振幅比法等)；②頻域分析法(Fourier變換、Walsh變換、最大熵譜等)；③模式識別法等[5-8]，其主要特點(diǎn)就是工作效率相對較高。但是它的各種誤判也不容忽視，比如在地震記錄中存在較強(qiáng)的干擾，采用自動拾取方法會產(chǎn)生誤判，難以獲取準(zhǔn)確的震相波至。因此，在地震記錄數(shù)據(jù)量非常大的情況下，只靠人工手動拾取不僅效率低，而且是不切合實(shí)際的；自動拾取方法雖然工作效率高，但對于存在較強(qiáng)噪音背景的地震事件的震相會產(chǎn)生誤判。因此，研究和開發(fā)一個穩(wěn)定、實(shí)用和高效的人機(jī)交互式天然地震數(shù)據(jù)處理系統(tǒng),具有一定的現(xiàn)實(shí)意義和應(yīng)用價值。

筆者使用Intel Fortran平臺進(jìn)行程序語言編譯[9]，基于windows系統(tǒng)搭建人機(jī)交互界面[10]，以便更簡潔、高效地實(shí)現(xiàn)天然地震數(shù)據(jù)的可視化顯示，更直觀通暢地實(shí)現(xiàn)人機(jī)交互[11-12]。同樣，通過人機(jī)交互界面實(shí)現(xiàn)地震震相的準(zhǔn)確識別及拾取等功能。

1 人機(jī)交互系統(tǒng)的總體設(shè)計

為了提高寬頻帶地震數(shù)據(jù)的處理速度及工作效率，筆者擬設(shè)計一套基于IASP91地球模型的天然地震人機(jī)交互數(shù)據(jù)處理系統(tǒng)，該系統(tǒng)采用人機(jī)交互方式及數(shù)據(jù)可視化技術(shù)，不僅能夠?qū)掝l帶地震數(shù)據(jù)進(jìn)行數(shù)據(jù)可視化、預(yù)處理、震相走時計算及拾取、數(shù)據(jù)的旋轉(zhuǎn)及重采樣等，而且能夠?qū)⑻幚斫Y(jié)果進(jìn)行保存、輸出，以便進(jìn)一步應(yīng)用于地球內(nèi)部結(jié)構(gòu)的研究。天然地震人機(jī)交互數(shù)據(jù)處理系統(tǒng)的總體設(shè)計如圖1所示。

圖1 天然地震人機(jī)交互數(shù)據(jù)處理系統(tǒng)總設(shè)計圖Fig.1 Site-plan of the interactive processing system

2 數(shù)據(jù)可視化及預(yù)處理

要實(shí)現(xiàn)天然地震數(shù)據(jù)的震相拾取及特征分析

等，首先需要將地震數(shù)據(jù)導(dǎo)入計算機(jī)內(nèi)存，然后根據(jù)其波形特征進(jìn)行人機(jī)交互的數(shù)據(jù)處理，完成數(shù)據(jù)分析；同時，天然地震數(shù)據(jù)記錄中可能存在壞道及嚴(yán)重的干擾等，需要對實(shí)際觀測的地震數(shù)據(jù)進(jìn)行濾波、壞道剔除等預(yù)處理。

2.1 地震數(shù)據(jù)的可視化

地震數(shù)據(jù)的高效可視化，可以提高地震數(shù)據(jù)處理工作的效率。筆者介紹的數(shù)據(jù)處理系統(tǒng)可以讀取寬頻地震儀記錄的多種格式(如SAC、DMX、CSS)等，同時能夠?qū)υ摂?shù)據(jù)進(jìn)行二維可視化。在進(jìn)行地震數(shù)據(jù)二維可視化過程中，首先要對地震數(shù)據(jù)進(jìn)行歸一化處理(道歸一)[13]，以避免各記錄地相互干擾，然后將地震數(shù)據(jù)進(jìn)行顯示。數(shù)據(jù)可視化的形式包括：①波形圖；②變面積圖。

1)波形圖是用曲線的形式表示一系列離散的地震數(shù)據(jù)。儀器采集的地震數(shù)據(jù)是一組離散的值，在屏幕坐標(biāo)中是一個個的離散點(diǎn)。在進(jìn)行地震波的波形圖繪制前常需要根據(jù)屏幕分辨率和采樣點(diǎn)的個數(shù)關(guān)系對地震數(shù)據(jù)進(jìn)行處理。如果采樣點(diǎn)較少就會出現(xiàn)折線，因此，在進(jìn)行地震曲線的繪制時需要對數(shù)據(jù)進(jìn)行插值使其圓滑。常用的插值方法有：①拉格朗日插值；②三次樣條插值；③二維插值等[15].文中的數(shù)據(jù)處理系統(tǒng)利用三次樣條算法進(jìn)行插值，在算法中只需要提供分段樣條曲線的控制點(diǎn)坐標(biāo)值，由節(jié)點(diǎn)內(nèi)部封閉的算法來完成樣條曲線的繪制。

例如已知控制點(diǎn)p1、p2、p3、p4，利用三次樣條函數(shù)，p1為起始點(diǎn)，p2 和p3為對應(yīng)控制點(diǎn)，p4為結(jié)束點(diǎn)。經(jīng)過三次樣條插值，既能使數(shù)據(jù)曲線圓滑，又可實(shí)現(xiàn)較小的插值誤差(圖2(a)和圖2(b))。

2)變面積圖是將地震波形數(shù)據(jù)的正值或負(fù)值區(qū)域用藍(lán)色筆刷進(jìn)行填充。在進(jìn)行地震數(shù)據(jù)顯示時，為了讓波形及同相軸更加清晰，易于分辨，常將地震數(shù)據(jù)的正值或者負(fù)值區(qū)域進(jìn)行填充(圖2(c))。

2.2 地震數(shù)據(jù)的預(yù)處理

一方面由于故障可能引起的地震儀器的工作異常，地震記錄形成壞道；另一方面，較強(qiáng)的干擾也會影響地震信號的信噪比，影響地震信號的識別與分析。這就需要對數(shù)據(jù)進(jìn)行必要的預(yù)處理，包括壞道剔除、濾波、增益控制等。

圖2 數(shù)據(jù)的可視化Fig.2 Visudization of seismic data(a)連接p1、p2、p3、p4的折線圖；(b)連接p1、p2、p3、p4三次樣條圖；(c)波形的變面積顯示

1)壞道剔除。壞道的存在不僅影響處理結(jié)果的準(zhǔn)確度，同時也使得系統(tǒng)的處理速度變慢。文中的數(shù)據(jù)處理系統(tǒng)使用人機(jī)交互的方式進(jìn)行壞道剔除，在地震波形圖中鼠標(biāo)左鍵點(diǎn)擊進(jìn)行壞道選擇，再執(zhí)行系統(tǒng)中的壞道剔除動作，即可完成壞道剔除；如果用戶發(fā)現(xiàn)錯選了良好的壞道，也可交互取消選中的道集，重新進(jìn)行壞道的拾取、剔除。值得注意的是，在選擇壞道過程中，即使每個臺站選擇任意一個分量，在剔除壞道過程中，不管屏幕顯示的是哪個分量，被選擇的壞道所在臺站的所有數(shù)據(jù)將從該次讀入中移除，不再參加處理。

2) 濾波處理。在濾波處理過程中，本文數(shù)據(jù)處理系統(tǒng)采用帶通濾波器(圖3)，并且利用快速傅里葉變換(FFT)對數(shù)據(jù)進(jìn)行濾波處理，以提高地震記錄的信噪比[16]。帶通濾波器的頻率響應(yīng)H(f)的數(shù)學(xué)表達(dá)式為式(1)。

H(f)=

(1)

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

圖3 帶通濾波器示意圖Fig.3 Schematic of band-pass filter

3)增益處理。為了使有效的弱信號更為明顯，使強(qiáng)弱不一的有效信號均能在屏幕上均勻地顯示，常需要對地震信號進(jìn)行增益控制。增益是對地震記錄的一種時變比例均衡，這種比例函數(shù)是根據(jù)所期望的規(guī)則確定的。如果地震數(shù)據(jù)的時間序列為f(t)，增益函數(shù)為c(t)，那么，經(jīng)增益處理后的地震記錄F(t)為式(2)。

(2)

本系統(tǒng)所使用的增益函數(shù)c(t)為一常數(shù)，并且用戶可以通過人機(jī)交互界面的地震記錄的干擾情況進(jìn)行分析，通過設(shè)置項(xiàng)選擇合適的值進(jìn)行設(shè)定。

3 系統(tǒng)功能

3.1 震相的拾取

天然地震信號中常見的地震震相有P、S、PKP、PKS、SKP、SKS、PKIKP、PKIKS、SKIKP、SKIKS以及LR等?；贗ASP91模型[14-15]，在已知震源點(diǎn)及接收臺站的情況下，利用打靶法這一傳統(tǒng)的射線追蹤方法[16]，即可獲得各震相的理論旅行時間，為震相人機(jī)交互識別和拾取提供參考?；诖怂枷?，我們的處理系統(tǒng)可快速計算某一地震對應(yīng)震相到達(dá)各臺站的時間，同時軟件對該時間點(diǎn)進(jìn)行標(biāo)識并顯示該標(biāo)識點(diǎn)于屏幕上，這也是地震記錄重新排列等功能的基礎(chǔ)。IASP91模型及地球內(nèi)部的地震震相如表1及圖4所示。

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

3.2 地震數(shù)據(jù)的旋轉(zhuǎn)

地震記錄有東西(EW)、南北(NS)、垂直(VZ)三個分量，在地震數(shù)據(jù)的研究工作(如接收函數(shù)、橫波分裂等)中，需要將地震波的兩個水平分量分解為徑向(SV波)及橫向(SH波)兩個分量。在地震事件已經(jīng)確定的前提下，利用臺站相對于震中的方位角β，可以將實(shí)際觀測的南北分量N(t)、東西分量E(t)進(jìn)行旋轉(zhuǎn)，得到徑向分量R(t)及橫向分量T(t)為式(3)。

表1 IASP91地球速度模型參數(shù)表

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

(3)

3.3 地震數(shù)據(jù)的合并及存儲

天然地震數(shù)據(jù)以臺站為單位存儲觀測的特點(diǎn)，使得數(shù)據(jù)的處理及分析，尤其是利用臺站間的互相關(guān)進(jìn)行地震波的震相確定造成很大的不便。故在設(shè)計該地震數(shù)據(jù)處理系統(tǒng)時，使得經(jīng)過系統(tǒng)的預(yù)處理、地震事件確定、震相識別等，將系統(tǒng)窗口內(nèi)顯示的全部數(shù)據(jù)合并存儲于DAT格式的文件中，并以該事件發(fā)生的年月日+時分秒+設(shè)定震相+擴(kuò)展名(.dat)進(jìn)行命名。該功能使得后期對該地震記錄的處理大大提高了效率。

4 系統(tǒng)的人機(jī)交互設(shè)計

系統(tǒng)的人機(jī)交互貫徹了數(shù)據(jù)輸入到數(shù)據(jù)處理完成的全過程(圖5)。

圖5 震相拾取過程的人機(jī)交互示意圖Fig.5 The interaction of phases picking

經(jīng)過了該系統(tǒng)對地震數(shù)據(jù)的導(dǎo)入、預(yù)處理，接著進(jìn)行人機(jī)交互的震相識別及拾取。在進(jìn)行各地震震相的識別和拾取時，需要響應(yīng)動作①：系統(tǒng)首先從準(zhǔn)備好的地震事件列表文件中，自動篩選出地震數(shù)據(jù)所在日期的所有地震事件，人工選擇地震記錄所對應(yīng)的地震事件；然后手動輸入需要識別的震相。這時，系統(tǒng)執(zhí)行對應(yīng)震相的識別計算，通過IASP91模型的速度參數(shù)，結(jié)合射線追蹤理論，便可確認(rèn)對應(yīng)震相在每個獨(dú)立臺站的到時。接著響應(yīng)動作②：機(jī)器拾取的震相需要人工判斷其震相識別正誤。如果識別正確，則進(jìn)行下一步的震相人機(jī)交互精確拾取，保存拾取結(jié)果，結(jié)束交互動作；如果識別錯誤，則返回動作①，從今選擇地震事件，直至系統(tǒng)識別到正確的震相。

5 應(yīng)用實(shí)例

選用南嶺地區(qū)天然地震流動臺站的觀測數(shù)據(jù)的一部分(圖6)，以日本福島外海2013年10月26日(北京東八區(qū)時間)7.1級地震為例，檢驗(yàn)該處理系統(tǒng)的震相識別及拾取效果。

5.1 參數(shù)設(shè)定

在系統(tǒng)的參數(shù)設(shè)定界面Paramerters Setting 中，對所需要讀取的數(shù)據(jù)段進(jìn)行設(shè)置(圖7)。其中包括數(shù)據(jù)段開始時刻、時段長度、濾波系數(shù)、分量顯示、波形放大、變面積圖顯示、記錄排序等設(shè)置項(xiàng)。由于福島地震距離臺網(wǎng)最近臺站約為3 000 km，到達(dá)臺站的時刻大約在地震發(fā)生后的1 100 s附近，故在示例中參數(shù)設(shè)定數(shù)據(jù)段長度設(shè)為30 min,起始時間當(dāng)延后與發(fā)震時刻，保證初至波及后續(xù)波都被讀取進(jìn)來。完成了上述的參數(shù)設(shè)置后，通過菜單欄的file->open，并且選擇相應(yīng)的數(shù)據(jù)格式文件，即可導(dǎo)入與之對應(yīng)時段的各臺站地震數(shù)據(jù)，并在屏幕上進(jìn)行數(shù)據(jù)的可視化。圖8為進(jìn)過預(yù)處理后的地震數(shù)據(jù)顯示。

圖6 南嶺地區(qū)天然地震流動臺陣Fig.6 Portable seismic array in Nanling

5.2 處理效果

完成了數(shù)據(jù)的預(yù)處理后，通過參數(shù)設(shè)置窗口，并且完成響應(yīng)動作①及響應(yīng)動作②后(圖4)，即可進(jìn)一步進(jìn)去的震相精確拾取界面(圖9)。

圖9可以看出，經(jīng)過地震記錄的震相識別、波形放大、變面積圖顯示，可以清楚地確定該地震的P震相的同相軸，紅色十字線處的同相軸提供了清晰的拾取標(biāo)準(zhǔn)。為了保證拾取的準(zhǔn)確性，用戶可進(jìn)行多次的拾取選定，且以對同一臺站的最后一次拾取為準(zhǔn)。這樣既可快速處理海量的地震數(shù)據(jù)，又保證了震相到時人工拾取的精確度。

圖7 參數(shù)設(shè)置及事件選擇Fig.7 Parameters setting and seismic event selection

圖8 預(yù)處理后的地震記錄Fig.8 Seismogram after bad traces elimination and filtering

圖9 人機(jī)交互拾取P震相示意圖Fig.9 Interactive picking of Seismic phases

6 結(jié)論

根據(jù)天然地震數(shù)據(jù)處理的需求，筆者基于Intel Fortran 開發(fā)平臺，在Windows環(huán)境下設(shè)計并實(shí)現(xiàn)了一套天然地震人機(jī)交互式的數(shù)據(jù)處理系統(tǒng)，該系統(tǒng)能夠?qū)崿F(xiàn)多種格式的地震數(shù)據(jù)導(dǎo)入、預(yù)處理、震相識別及到時拾取等功能。實(shí)際測試結(jié)果表明，該人機(jī)交互系統(tǒng)操作方便，效果良好，為天然地震數(shù)據(jù)的處理及分析提供了一個高效、靈活、方便的工具。

[1] 傅淑芳,劉寶成,李文藝.地震學(xué)教程[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].地震學(xué)報,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] 高淑芳,李山有.一種改進(jìn)的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] 劉衛(wèi)國,蔡旭輝.Fortran 90程序設(shè)計教程[M].北京：北京郵電大學(xué)出版社,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]張宏林,孔艷,王哲．按實(shí)例學(xué) 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]李樂山.人機(jī)界面設(shè)計[M].北京：科學(xué)出版社.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]葛文庚.數(shù)據(jù)可視化及其在地學(xué)中的應(yīng)用[J].工程地球物理學(xué)報,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]王有學(xué).地震波理論基礎(chǔ)[M].北京：地質(zhì)出版社,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

國家自然科學(xué)基金項(xiàng)目(41574030,41174077，41674075)；廣西自然科學(xué)基金重點(diǎn)基金項(xiàng)目(2016GXNSFDA380014);大陸構(gòu)造與動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室開放課題(k201505)

徐金榮(1989-)，男，碩士，研究方向?yàn)楣腆w地球物理，E-mail：10360928@qq.com。

王有學(xué)(1961-)，男，博士，教授，主要研究方向?yàn)榈厍蛭锢怼⒌厍騼?nèi)部結(jié)構(gòu)及動力學(xué)，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