一種基于模板匹配的震相關聯和提取技術及其初步應用

龍鋒 祁玉萍 趙敏 芮雪蓮

四川省地震局，成都 610041

0 引言

隨著地震學的發展和計算機能力的提升，以及對地震學研究結果精細化的需求越來越高，實測地震學的發展趨勢將是廣泛采用海量臺站、小孔徑、寬頻帶、高動態范圍及高采樣率的設備進行觀測。在此基礎上產生的海量數據除了對存儲造成壓力外，地震事件的識別、震相的讀取以及地震定位等傳統依賴人工判別的步驟也會成為極為耗時耗力的工作，從而使后續的地震學研究受到影響。

近幾年發展了眾多計算機自動識別及定位的方法。按其方法大體可以分為兩類，其中一類是以已知地震為模板，采用互相關等信號增強技術來識別較小的地震事件。這種方法實現較為簡單，運行快速，但缺點在于其僅能識別出模板地震的重復/相似事件。由于其操作簡便，被大量應用于遺漏地震檢測中(Gibbons et al，2006；Peng et al，2006、2009；Shelly et al，2007；譚毅培等，2014；尹欣欣等，2020)。第二類為自主識別法，其通過一定的算法識別不同震相具有的特定形態，從而進行地震的識別與定位(Kao et al，2004、2007；Gharti et al，2010；Liao et al，2012)，但這種方法約束性不強，容易出現虛報與漏報。近年來出現的神經網絡方法可通過事前的大樣本量學習，提取出地震信號和噪聲的特有表征，從而避免了不斷增長的模板波形庫。在實際運算過程中，地震信號的識別被轉化為機器學習中的監督二元分類問題(Perol et al，2018)，這種方法不受模板的制約，是一種相對“普適”的地震自動識別方法。然而大樣本量的機器學習依賴高昂的硬件設備，模板匹配仍是目前微小地震識別所采用的主要方法。

常規模板匹配基于重復/相似地震，即假定他們在位置上是相同的，不同臺站所記錄到的與模板具有高度相關的事件的時間偏移量被歸算為待識別地震的發震時間(Peng et al，2006)；也有算法考慮到微小地震與模板之間可能存在空間上的差異，在確定發震時刻的同時，通過類似于三維格點搜索等方法來確定待識別地震的空間位置(Zhang et al，2015)。但目前這些方法存在一個共同的缺點，即不能提供屬于每個地震事件的震相到時信息，而到時信息是地震精確定位和體波成像的基礎資料。理論上，當模板與待識別地震之間存在一一對應關系時，可以根據達到給定閾值的相關系數的時間偏移量來給出各臺站的震相到時。然而現實情況是，一方面，待識別地震往往震級偏小，頻發的小震會使得觸發信號短時內在同一個臺站多次出現；另一方面，在通用的模板匹配識別微震的過程中，地震信號會表現為互相關信號的時間序列，而如何區分這些觸發的信號屬于同一個地震是難點。

21世紀初以來，龍門山斷裂帶上陸續發生了2008年汶川8.0級地震與2013年蘆山7.0級地震，2個地震序列之間所在的龍門山斷裂帶南段存在一條長約40km的地震空區(杜方等，2013；高原等，2013；何富君等，2017；梁春濤等，2018)，如圖1 所示。對于其未來的強震風險，學術界有眾多爭論： 庫侖應力結果顯示，空區位于2個強震所造成的應力擾動增強區，存在發生強震的風險(Li et al，2014；Liu et al，2014)；而介質成像結果顯示，空區內介質強度不高，難以積累強震所需的足夠應力(Pei et al，2014；顏照坤等，2014)。基于區域臺網的地震定位結果研究，認為空區內地震雖然偏少且弱，但基本上仍表現為雙石-大川斷裂等主要構造在活動(芮雪蓮等，2020)，也有基于密集臺陣觀測的模板匹配研究，結果認為空區存在明顯的微震“虧空”，否定了該區域通過微地震釋放累計應力的調節模式(王朝亮，2019)。本文以龍門山斷裂帶南段的地震空區為研究對象，發展一種提取單個地震震相到時的技術流程并對其進行驗證，而對于空區內所有模板地震的全面匹配，還有待于后續研究。

圖 1 龍門山斷裂帶南段構造及本文所使用模板和臺站分布

1 數據與方法

收集并挑選了2013年5月至2016年12月空區內35次1.5≤ML≤3.6高信噪比的地震事件波形，共計使用84個臺站(圖1)。按照Zhang等(2015)相同的操作，對模板及連續波形均進行了2～8Hz的Butterworth帶通濾波，并將波形的采樣率降至25sps，從而在奈奎斯特頻率范圍內保留足夠信息的同時，有效減少計算量。定義模板中P波和S波震相的前1s和后3s為模板窗，采用時間域互相關方法計算各模板窗在連續波形中的歸一化相關系數NCC。為了協同三分量處理，將同一個臺站的三分量截取相同的模板窗，并將相關系數疊加后平均(圖2)。從圖2(c)中可以看出，平均化處理后的NCC能更好地壓低噪聲，從而提高匹配時的信噪比。

圖 2 模板匹配處理流程(以20160215023958.20模板事件的JJS臺記錄為例)(a)模板窗；(b)三分量連續波形片段；(c)三分量NCC；(d)平均NCC

地震事件發生后多個臺站會在短時內記錄到震相，從而使得震相密度記錄在時間上發生變化，因此密度聚類算法可用于對疑似的同一個地震震相進行關聯。基于密度的聚類算法是數據挖掘技術中被廣泛應用的一類方法，此類算法假定聚類結構能通過樣本分布的緊密程度確定。通常情況下，密度聚類算法從樣本密度的角度來考察樣本之間的可連接性，并基于可連接性樣本不斷擴展聚類簇，以獲得最終的聚類結果。其優點在于既可以找出不同形狀的聚類，且事先不需要知道聚類的個數(Ester et al，1996；Schubert et al，2017；周志華，2016)。

震相到時僅涉及時間數據，因此可表達為一個簡單的一維聚類問題。我們選擇了“帶噪聲的基于密度的空間聚類算法”(Density-Based Spatial Clustering of Application with Noise，DBSCAN)，這是一種著名的密度聚類算法，其基于一組鄰域參數(eps，minPts)來刻畫樣本分布的緊密程度，其中eps定義了簇的尺度上限，即一簇之內任意樣本之間的距離不得大于eps。在震相聚類之前，我們統計發現35個模板中第一個與最后一個震相的到時差均限定在30s以內，因此eps設定為30s；minPts為定義一個獨立簇的最小樣本個數。考慮到至少需要4個臺站才能進行地震定位，而極限條件下每個臺站僅提供1個震相，因此minPts設定為4。需要注意的是，這里其實設定的是震相個數，而不是臺站個數，因為1個臺站可能會有2個震相，因此即便是被算法識別為單獨的簇(即1個地震)，但臺站數量會少于4個，所以其仍不能被定位。

震源位置的空間差異會在震相到時和NCC中體現，采用NCC=0.7作為疑似震相信號觸發閾值，這樣既可以保證疑似信號的準確性，也能使得待檢測信號能在模板事件周圍有足夠的搜尋空間。通過互相關的計算，共計得到資料時段內1300萬余次觸發。采用上述的eps=30、minPts=4為參數進行DBSCAN聚類計算，獲得39萬次聚類簇，其中770萬次觸發參與了聚類，未參加聚類的離群點可能是太過細小的地震達不到聚類標準或隨機噪聲。圖3(a)展示了20160215023958.20模板事件在2016年3月2日的聚類結果，當天存在2次離群點，而具體的聚類簇中的細節展示了各臺站各震相到時的先后順序(圖3(b))。

圖 3 基于DBSCAN的震相聚類(a)以20160215023958.20模板事件在2016年3月2日連續波形上的聚類結果為例，其中藍色為成功的聚類簇，灰色為離群點；(b)成功聚類簇時間段內的細節展示(圖(a)中的虛線框部分)，字符分別為臺站名、震相名和NCC

近一半的觸發未能構成聚類簇，說明離群點在時間上是廣泛存在的，不排除存在某些離群點與鄰近的觸發恰好構成了聚類簇。通過抽樣調查，部分簇中存在如下典型異常：①重復的臺站記錄；②定位誤差大。其中前者是由于短時內同一個臺站多次觸發所致，通過簇內遍歷并選取定位誤差最小所對應的到時可將其解決；而后者則反映了簇內一個或多個觸發并不屬于同一個地震事件，我們采用“自洽性檢驗”的方法來挑選簇內屬于同一地震的觸發。

以具有20個震相記錄的20160215023958.20模板事件為例，其具體方法為：①計算模板中各震相與其余震相的到時差，構建到時差矩陣(圖4(a))，模板的到時差矩陣定義了每一個震相到時在這個二維坐標中的位置與幅度，是一個相對時標框架。從圖4(a)中可以看出，該模板各震相到時基本限定在30s內；②計算與該模板相對應的聚類簇的到時差矩陣，并與模板的到時差矩陣做差，形成到時雙差矩陣(圖4(b))，其反映了簇時標與模板時標的差異，能直觀顯示出哪些震相到時與模板之間存在較大偏離。圖 4(b)顯示該聚類簇具有14個震相記錄(無記錄的震相用灰色表示)，各震相到時差要普遍小于模板，表明在不進行約束的情況下，相比于模板，該簇所反映的地震震中更接近于所用到的臺站的幾何中心，但這些臺站幾乎都偏向一側。而部分震相出現到時雙差偏移量超過10s的現象，說明疑似地震與模板時間之間存在幾十甚至上百千米的距離，但由于這二者是通過高相關匹配得到，理論上不會相距太遠，因此那些具有較大偏移量的到時雙差震相有很大概率是混入簇中的誤觸發，應該予以排除。在剔除掉超過2s的到時雙差震相后，到時雙差范圍限定在1s之間(圖4(c))， 且正向和反向偏移的震相數量相當，幅度較為一致，說明疑似地震與模板事件之間距離較近，滿足這樣關系的震相到時則說明他們是自洽的。

圖 4 模板到時差矩陣(a)模板到時差矩陣；(b)聚類簇與模板的到時雙差矩陣；(c)自洽性檢驗后的到時雙差矩陣；灰色部分表示無記錄和未挑選的震相

但到目前為止，仍不能確定被保留下來的震相是全部可用的，圖4(c)顯示1號震相可與13個震相構成自洽關系，其中就包括2號震相，但2號震相僅能與4個震相構成自洽關系，表明部分震相在自洽性上離其他震相更遠。因此，對圖4(c)中可構成自洽震相的數量，按由多到少的順序統計了其觀測數，結果如圖5 所示，圖中顯示所有的震相均能確保有4個以上臺站參與記錄，且具備7個以上自洽記錄占大多數。因此，我們選取了具備7個以上自洽記錄的震相并將其編號，按由多到少的順序求取其交集，最后得到的集合即為該疑似事件中被各震相廣泛自洽的到時信息。

圖 5 各震相的觀測統計柱狀圖

2 結果分析

按照以上規則流程，從39萬次聚類簇中獲得了118次疑似事件的震相到時信息，包含504條P波震相和984條S波震相。利用HYPOINVERSE2000(Klein，1989)對疑似事件進行了絕對定位，少量事件由于震相較少，約束不夠導致誤差較大，從定位結果中剔除了14次走時殘差大于1s或水平誤差大于3km的疑似事件。考慮到部分模板可能高度相似，因此會識別出重復的疑似事件。將發震時刻在3s以內且震中距小于5km的地震事件作為重復事件，在104個疑似事件中識別出10組共計11次重復事件，每一組僅保留水平誤差最小的結果。按照該流程，最終獲得93個確切的地震信息，為所使用模板的2.65倍，包含390條P波震相以及791條S波震相。

地方震震級可通過測量各記錄臺站的S波最大振幅獲得，已知對于某個臺站，震級計算公式為(中國地震局，2001；劉瑞豐等，2015)

Mo=lg(Ao)+R(do)

(1)

其中，Mo為模板地震震級，Ao為該臺站記錄到的S波最大振幅，R(do)為臺站震中距為do時的量規函數，同樣地，自動識別出的地震震級Mx可寫為

Mx=lg(Ax)+R(dx)

(2)

其中，Ax為與模板同一個臺站記錄到的S波最大振幅，dx為自動識別出的地震與該臺的震中距。兩式相減，并考慮到do≈dx，可得

Mx=lg(Ax/Ao)+Mo

(3)

S波最大振幅由2個水平向波形的最大振幅算術平均求得，利用式(3)計算同一個疑似地震每個臺站的震級，最后求其平均值作為整個疑似事件的震級。結果顯示，采用上述流程在該研究區的震級識別下限為0.5級左右，范圍在0.5～1.0級之間。原則上，如果降低判別標準(如在自洽性檢驗中降低所需觀測數)，可進一步降低識別的震級下限數值。由圖6 展示的M-t圖可以看出，識別出的地震震級普遍比模板震級低，主要是由于高震級地震本就可以通過臺網波形目視分析得到。

圖 6 研究區內的M-t圖

從所有地震的時間演化圖像(圖6)來看，2014年后地震發生率總體比較平穩，但2013年6月開始至2013年年底未識別到任何地震。這可能是由于該區域的確處于地震平靜期，也有可能是所選擇的模板不合適，恰好在這個時段未發生和模板匹配的地震。

由于通過模板匹配獲得了更多的地震樣本，使得采用雙差法(Waldhauser et al，2000)對地震進行相對定位成為可能。考慮到地震總體分布彌散，在定位過程中，將組隊的搜尋半徑設置為20km，并利用龍門山斷裂帶南段的速度模型(Long et al，2015)采用LSQR進行2次迭代后得到研究區最終的地震精定位結果，定位誤差為水平向0.6km，垂直向0.9km，走時殘差0.13s。

精定位后的震中分布(圖7)顯示，模板地震大多沿NE-SW走向的雙石-大川斷裂兩側分布，其中絕大多數模板并未在其周圍地區識別出小震。被識別出的小震主要集中分布在雙河南東和安順以西，其中后者由51個地震(3個模板，48個識別地震)組成，形成一個走向 NE-SW、 與雙石-大川斷裂平行的長約2km的條狀。為了研究地震在深部的分布特征，沿垂直和平行于區域斷裂走向的方向劃分了2個剖面。垂直于區域斷裂走向的AA′剖面(圖8(a))顯示地震密集分布在2個區域，其中雙河以西的一叢地震震級偏大，以2～3級地震為主，幾個稍大地震的剖面連線顯示該斷面傾向NW，最深處可達18km。盡管選擇的模板數量較少，按照這叢地震的傾向與地表斷裂出露點的位置關系判斷，其反映的應該是雙石-大川斷裂的位置和產狀。另一叢地震位于安順以西，分布密集，震級偏小，大部分為根據模板識別出的地震。剖面顯示深度小于10km時地震分布密集，而深于10km地震較少，但有2次2級左右地震。如果用曲線連接淺部的密集區和更深處的2次較大地震，則可勾勒出一個淺部直立、隨著深度增加傾向NW的“鏟狀”構造。但地表并無已知斷裂對應(圖7)，應該為某條未知的斷裂。平行于區域斷層走向的深度剖面(圖8(b))顯示震源深度有自SW(圖8(b)左側)向NE(圖8(b)右側)變深的趨勢，而2013年蘆山7.0級地震序列也展現了類似的分布特征(Long et al，2015)，可能反映了龍門山斷裂帶南段不同段落具有類似的構造屬性。

圖 7 研究區精定位后的震中分布

圖 8 不同方向的深度剖面

3 結論

本文設計了一套基于模板匹配的震相關聯和提取流程，不僅可對高度匹配的疑似地震進行時空強參數的測算，也能提供每個疑似事件對應的震相到時信息，可供后續的速度結構成像、臺站校正等分析使用。

以龍門山斷裂帶南段的地震空區為例，按照該流程，在研究區中挑選了35個模板，對2013年5月至2016年12月的連續波形進行了互相關計算，并設置NCC=0.7，共記錄到1300萬余次觸發。通過采用eps=30、minPts=4為參數進行DBSCAN聚類計算，共獲得39萬個聚類簇；在自洽性檢驗過程中，保留最少7個以上自洽記錄的震相，得到了118次疑似事件。去除重復地震后得到93個確切的地震信息，為所使用模板的2.65倍，包含390條P波震相以及791條S波震相。識別出的地震震級在0.5～1.0之間，2014年后的地震發生率比較均勻，而2013年下半年的平靜可能與研究區真實閉鎖或模板的挑選有關。

對所有地震進行精定位后的結果顯示，絕大部分模板并未識別出自身的重復/相似地震，而被識別出的地震是由少量模板匹配成功的，主要分布在雙河和安順附近。垂直于區域斷裂走向的深度剖面可以清晰顯示出傾向NW的雙石-大川斷裂的位置和規模，而安順以西一叢被識別出的地震呈現出淺部直立、深部傾向NW的“鏟狀”分布形態，反映了某條不知名斷裂的幾何特征。平行于區域斷裂走向的深度剖面顯示震源深度具有與蘆山7.0級地震序列類似的自SW向NE變深的趨勢，可能反映了龍門山斷裂帶南段不同段落具有類似的構造特征。

由龍門山斷裂帶南段地震空區內的模板識別測試結果可以看出，本文設計的流程可提供包括震源參數和每個地震震相到時信息在內的完整的震相報告，在模板識別領域可作為一套計算機自主處理的可行方案。