基于分形盒維算法的天然地震波初至拾取方法研究

張 琪， 王有學(xué)， 熊 彬， 蔣嬋君， 王心宇， 曾 成

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

0 引言

在天然地震數(shù)據(jù)的處理中，拾取初至信號(hào)是地震波走時(shí)成像、研究地球內(nèi)部結(jié)構(gòu)的基礎(chǔ)。實(shí)際工作中，根據(jù)地震波的振幅和波形特征，利用人工拾取初至信號(hào)。然而，天然地震波含有各種干擾，人工提取存在的主觀偏見(jiàn)給提取初至信號(hào)帶來(lái)了極大的不唯一性。因此，以計(jì)算機(jī)為基礎(chǔ)的自動(dòng)拾取方法毋容置疑地可以給我們提供更為準(zhǔn)確的地震波旅行時(shí)間信息。

對(duì)于自動(dòng)拾取地震波初至到時(shí)，地震科學(xué)家們提出了許多方法，諸如相關(guān)法、神經(jīng)網(wǎng)絡(luò)法、分形法等[1]。Gelchinsky等[2]提出了一種基于地震信號(hào)的特定相關(guān)屬性來(lái)確定初至的方法;王繼等[3]用多道互相關(guān)和最小二乘方法確定遠(yuǎn)震震相初至?xí)r間;潘樹(shù)林等[4]利用優(yōu)化相關(guān)法拾取地震事件的初至信號(hào)，這些方法均基于地震道的整體特征，對(duì)噪聲有較好的抑制作用，但其拾取結(jié)果仍然受到續(xù)至波的影響較大。Murat等[5]提出了在低信噪比的信號(hào)中反向追蹤神經(jīng)網(wǎng)絡(luò)技術(shù)能在噪聲干擾背景下識(shí)別初至，但神經(jīng)網(wǎng)絡(luò)拾取初至學(xué)習(xí)速度慢，實(shí)現(xiàn)復(fù)雜，需要綜合考慮網(wǎng)格結(jié)構(gòu)、訓(xùn)練方法等問(wèn)題[6];美籍?dāng)?shù)學(xué)家曼德布羅特[7]提出分形思想并于1975年創(chuàng)建分形理論，分形盒維數(shù)(Fractal Box)作為分形維數(shù)的一種定義，它的大小及變化可以反映信號(hào)的不規(guī)則度和復(fù)雜度。這種方法容易理解，易于數(shù)學(xué)計(jì)算且直觀形象，所以分形理論在解決自動(dòng)提取地震波初至到時(shí)的問(wèn)題上發(fā)揮了長(zhǎng)足的作用[8-10]。

對(duì)地震信號(hào)的抽樣序列來(lái)說(shuō)，將其作無(wú)量綱圖形化處理可得到一條自仿射曲線，具有不規(guī)則性和近似的自相似性,可以用分?jǐn)?shù)維來(lái)描述其不規(guī)則程度[11]。研究發(fā)現(xiàn)，在地震波初至到達(dá)之前，地震監(jiān)測(cè)系統(tǒng)監(jiān)測(cè)到的地震信號(hào)時(shí)間序列為隨機(jī)噪聲的響應(yīng)，地震波初至到達(dá)之后，地震信號(hào)時(shí)間序列為隨機(jī)噪聲與地震波動(dòng)的復(fù)合響應(yīng)[12]。由于地震信號(hào)的分段特征，在地震波初至到達(dá)前后地震波時(shí)間序列的分形特征會(huì)發(fā)生突變，根據(jù)其特征的變化可以定量地確定初至波的到達(dá)時(shí)間[13]。

筆者基于分形盒維算法(FBM- Fractal Box Method)提出了一種天然地震波初至拾取的方法，該方法能夠較好地解決自動(dòng)拾取地震波初至的問(wèn)題。

1 分形盒維算法的基本原理

維的定義基于δ規(guī)模下的測(cè)量法的思想。設(shè)X是Rn中任一非空有界子集，對(duì)于一段離散的數(shù)字信號(hào)x(i) (i=1，2，…，N)，x(n)∈X，將一系列尺度為δ(i) (i=1，2，…，M)的方網(wǎng)格對(duì)X進(jìn)行覆蓋(圖1)，得到有效覆蓋網(wǎng)格數(shù)Nδ(i)(X) (i=1，2，…，M)，則X的上、下盒維定義[14]分別為

(1)

(2)

如果X的上、下盒維相等，則X盒維數(shù)為

(3)

離散的數(shù)字信號(hào)無(wú)法求得式(3)中的極限值，實(shí)際計(jì)算中，通過(guò)最小二乘法得到lnNδ(i)(X)～ lnδ(i)的擬合直線，其斜率就是的分形盒維數(shù)p。

圖1 數(shù)字信號(hào)分形維數(shù)的計(jì)算原理示意圖Fig.1 Calculation principle of fractal dimension on digital signal

對(duì)于地震波記錄來(lái)說(shuō)，它是一個(gè)離散序列x(i) (i=1，2，…，N)，我們?cè)谝粋€(gè)相對(duì)固定的時(shí)窗Wk(k=1，2，…，N)(長(zhǎng)度為l)內(nèi)，將時(shí)窗M等分，一般M取2的冪次方。以信號(hào)零值為零線，以相同的等分間隔畫(huà)平行于零線的平行線，使得最外側(cè)的平行線包括了整個(gè)時(shí)窗內(nèi)的地震信號(hào)，此時(shí)，在時(shí)窗內(nèi)形成一個(gè)邊長(zhǎng)為l/(M-1)的正方形陣列。然后對(duì)時(shí)窗內(nèi)的地震信號(hào)以Δt=l/(M-1)進(jìn)行重采樣，并對(duì)時(shí)窗內(nèi)重新采樣的地震信號(hào)計(jì)算與一系列尺度為δ(i)=l/(M-1)*(i-1)對(duì)應(yīng)的Nδ(i)，最后求得該時(shí)窗的分形盒維數(shù)pk。按照地震記錄的采樣間隔移動(dòng)時(shí)窗，重復(fù)上述步驟，得到每一個(gè)時(shí)窗所對(duì)應(yīng)的分形盒維數(shù)，直到覆蓋整個(gè)信號(hào)序列。

2 時(shí)窗大小對(duì)分形盒維的影響

通過(guò)對(duì)每個(gè)時(shí)間窗口下進(jìn)行分形盒維數(shù)的計(jì)算，可得到地震波信號(hào)的分形盒維數(shù)軌跡。但是，時(shí)窗長(zhǎng)度過(guò)小，計(jì)算量增大；時(shí)窗長(zhǎng)度過(guò)大，降低了信號(hào)識(shí)別能力。針對(duì)時(shí)窗長(zhǎng)度選取及它對(duì)分形維的影響,又分別選取時(shí)窗長(zhǎng)度l=0.3 s、0.6 s、1 s、2 s這四個(gè)取值進(jìn)行測(cè)試，其結(jié)果如圖2所示。

由圖2可以看出，當(dāng)時(shí)窗長(zhǎng)度取0.6 s時(shí)，分形盒維數(shù)曲線在有效信號(hào)附近下降趨勢(shì)明顯(圖2(c))，計(jì)算速度合適。

3 分形盒維在地震波初至拾取中應(yīng)用

根據(jù)分形盒維方法，我們對(duì)廣西地區(qū)寬頻帶地震流動(dòng)臺(tái)網(wǎng)記錄的地震事件 (表1、圖3(b)、圖4(b)、圖5(b)、圖6(b))進(jìn)行盒維分析。地震波采樣間隔10 ms，時(shí)窗長(zhǎng)度取0.6 s。為了求取分形盒維數(shù)的變化率，對(duì)分形盒維數(shù)曲線進(jìn)行了平滑[15](圖3(c)、圖4(c)、圖5(c)、圖6(c))，其計(jì)算結(jié)果如圖3～圖6(b)～圖6(d)所示。

由圖3可以看出，在10 s以前分形盒維數(shù)p在1.5左右，而在10 s以后，分形盒維數(shù)p急劇降至1.2以下，由分形盒維數(shù)p的變化率(圖3(d))自動(dòng)拾取初至?xí)r間為t=9.51 s；在圖4中，在3 s以前分形盒維數(shù)p在1.5左右，而在3 s以后，分形盒維數(shù)p明顯地降至1.4左右，由分形盒維數(shù)p的變化率(圖4(d))自動(dòng)拾取初至?xí)r間為t=3.08 s。在圖5和圖6中，分形盒維數(shù)p具有明顯變化，由分形盒維數(shù)p的變化率(圖5(d)與圖6(d))自動(dòng)拾取初至?xí)r間分別為t=7.70 s、t=18.91 s。

地震事件時(shí)間/UTC經(jīng)度/°E緯度/°N震源深度/km震級(jí)震中距/°馬里亞納群島地震2016/7/29 21:18:26145.5418.54207.627.733.47廣西梧州地震2016/7/31 9:18:13111.4724.1322.554.92.05俾斯麥群島地震2016/8/31 3:11:35152.79-3.69476.006.849.77日本本州地震2016/11/21 20:59:50141.3537.3914.696.929.57

圖3 馬里亞納群島地震記錄及其分形計(jì)算結(jié)果Fig.3 Seismic data from earthquake in Mariana islands and computed results by FBM(a)地震信號(hào)；(b)地震信號(hào)的分形盒維數(shù)；(c)平滑后的分形盒維數(shù)；(d)分形盒維數(shù)的變化率

圖4 廣西梧州地震記錄及其分形計(jì)算結(jié)果Fig.4 Record data from earthquake in Wuzhou and computed results by FBM(a)地震信號(hào)；(b)地震信號(hào)的分形盒維數(shù)；(c)平滑后的分形盒維數(shù)；(d)分形盒維數(shù)的變化率

圖5 俾斯麥群島地震記錄及其分形計(jì)算結(jié)果Fig.5 Seismic data from earthquake in Bismarck islands and computed results by FBM(a)地震信號(hào)；(b)地震信號(hào)的分形盒維數(shù)；(c)平滑后的分形盒維數(shù)；(d)分形盒維數(shù)的變化率

圖6 日本本州地震記錄及其分形計(jì)算結(jié)果Fig.6 Record data from earthquake in Japan's Honshu and computed results by FBM(a)地震信號(hào)；(b)地震信號(hào)的分形盒維數(shù)；(c)平滑后的分形盒維數(shù)；(d)分形盒維數(shù)的變化率

4 結(jié)論

根據(jù)分形理論，對(duì)基于分形理論的地震波初至拾取的方法進(jìn)行了研究，主要結(jié)論如下：

1)分形維理論可以根據(jù)地震波初至到達(dá)前后地震波時(shí)間序列的分形特征突變，定量地確定初至波的到達(dá)時(shí)間，能夠較好地解決自動(dòng)拾取地震波初至的問(wèn)題。

2)基于分形理論的地震波初至?xí)r間的拾取方法易于數(shù)學(xué)計(jì)算，其拾取精度可以達(dá)到0.01 s；拾取有唯一點(diǎn)，解決了手動(dòng)拾取帶來(lái)的不唯一性。對(duì)于采樣間隔10.0 ms，記錄長(zhǎng)度在30 s之內(nèi)的地震信號(hào)，耗時(shí)為3.42 s，以后可以通過(guò)改進(jìn)算法更好的縮短拾取耗時(shí)。

3)試驗(yàn)結(jié)果表明，當(dāng)窗口大小為0.6 s時(shí)，計(jì)算速度合適，可以保證地震波波至?xí)r間的拾取精度。