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

張 琪， 王有學， 熊 彬， 蔣嬋君， 王心宇， 曾 成

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

0 引言

在天然地震數據的處理中，拾取初至信號是地震波走時成像、研究地球內部結構的基礎。實際工作中，根據地震波的振幅和波形特征，利用人工拾取初至信號。然而，天然地震波含有各種干擾，人工提取存在的主觀偏見給提取初至信號帶來了極大的不唯一性。因此，以計算機為基礎的自動拾取方法毋容置疑地可以給我們提供更為準確的地震波旅行時間信息。

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

對地震信號的抽樣序列來說，將其作無量綱圖形化處理可得到一條自仿射曲線，具有不規則性和近似的自相似性,可以用分數維來描述其不規則程度[11]。研究發現，在地震波初至到達之前，地震監測系統監測到的地震信號時間序列為隨機噪聲的響應，地震波初至到達之后，地震信號時間序列為隨機噪聲與地震波動的復合響應[12]。由于地震信號的分段特征，在地震波初至到達前后地震波時間序列的分形特征會發生突變，根據其特征的變化可以定量地確定初至波的到達時間[13]。

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

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

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

(1)

(2)

如果X的上、下盒維相等，則X盒維數為

(3)

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

圖1 數字信號分形維數的計算原理示意圖Fig.1 Calculation principle of fractal dimension on digital signal

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

2 時窗大小對分形盒維的影響

通過對每個時間窗口下進行分形盒維數的計算，可得到地震波信號的分形盒維數軌跡。但是，時窗長度過小，計算量增大；時窗長度過大，降低了信號識別能力。針對時窗長度選取及它對分形維的影響,又分別選取時窗長度l=0.3 s、0.6 s、1 s、2 s這四個取值進行測試，其結果如圖2所示。

由圖2可以看出，當時窗長度取0.6 s時，分形盒維數曲線在有效信號附近下降趨勢明顯(圖2(c))，計算速度合適。

3 分形盒維在地震波初至拾取中應用

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

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

地震事件時間/UTC經度/°E緯度/°N震源深度/km震級震中距/°馬里亞納群島地震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 馬里亞納群島地震記錄及其分形計算結果Fig.3 Seismic data from earthquake in Mariana islands and computed results by FBM(a)地震信號；(b)地震信號的分形盒維數；(c)平滑后的分形盒維數；(d)分形盒維數的變化率

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

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

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

4 結論

根據分形理論，對基于分形理論的地震波初至拾取的方法進行了研究，主要結論如下：

1)分形維理論可以根據地震波初至到達前后地震波時間序列的分形特征突變，定量地確定初至波的到達時間，能夠較好地解決自動拾取地震波初至的問題。

2)基于分形理論的地震波初至時間的拾取方法易于數學計算，其拾取精度可以達到0.01 s；拾取有唯一點，解決了手動拾取帶來的不唯一性。對于采樣間隔10.0 ms，記錄長度在30 s之內的地震信號，耗時為3.42 s，以后可以通過改進算法更好的縮短拾取耗時。

3)試驗結果表明，當窗口大小為0.6 s時，計算速度合適，可以保證地震波波至時間的拾取精度。