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

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

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

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

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

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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

H(f)=

(1)

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

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

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

(2)

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

3 系統(tǒng)功能

3.1 震相的拾取

天然地震信號中常見的地震震相有P、S、PKP、PKS、SKP、SKS、PKIKP、PKIKS、SKIKP、SKIKS以及LR等。基于IASP91模型[14-15]，在已知震源點及接收臺站的情況下，利用打靶法這一傳統(tǒng)的射線追蹤方法[16]，即可獲得各震相的理論旅行時間，為震相人機交互識別和拾取提供參考。基于此思想，我們的處理系統(tǒng)可快速計算某一地震對應震相到達各臺站的時間，同時軟件對該時間點進行標識并顯示該標識點于屏幕上，這也是地震記錄重新排列等功能的基礎。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)確定的前提下，利用臺站相對于震中的方位角β，可以將實際觀測的南北分量N(t)、東西分量E(t)進行旋轉(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ù)以臺站為單位存儲觀測的特點，使得數(shù)據(jù)的處理及分析，尤其是利用臺站間的互相關進行地震波的震相確定造成很大的不便。故在設計該地震數(shù)據(jù)處理系統(tǒng)時，使得經(jīng)過系統(tǒng)的預處理、地震事件確定、震相識別等，將系統(tǒng)窗口內(nèi)顯示的全部數(shù)據(jù)合并存儲于DAT格式的文件中，并以該事件發(fā)生的年月日+時分秒+設定震相+擴展名(.dat)進行命名。該功能使得后期對該地震記錄的處理大大提高了效率。

4 系統(tǒng)的人機交互設計

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

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

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

5 應用實例

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

5.1 參數(shù)設定

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

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

5.2 處理效果

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

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

圖7 參數(shù)設置及事件選擇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 結論

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

[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] 劉衛(wèi)國,蔡旭輝.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]葛文庚.數(shù)據(jù)可視化及其在地學中的應用[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].北京：地質(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

國家自然科學基金項目(41574030,41174077，41674075)；廣西自然科學基金重點基金項目(2016GXNSFDA380014);大陸構造與動力學國家重點實驗室開放課題(k201505)

徐金榮(1989-)，男，碩士，研究方向為固體地球物理，E-mail：10360928@qq.com。

王有學(1961-)，男，博士，教授，主要研究方向為地球物理、地球內(nèi)部結構及動力學，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