龍鋒 祁玉萍 趙敏 芮雪蓮

四川省地震局，成都 610041

0 引言

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

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

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

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

圖 1 龍門山斷裂帶南段構(gòu)造及本文所使用模板和臺(tái)站分布

1 數(shù)據(jù)與方法

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

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

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

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

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

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

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

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

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

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

圖 5 各震相的觀測統(tǒng)計(jì)柱狀圖

2 結(jié)果分析

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

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

Mo=lg(Ao)+R(do)

(1)

其中，Mo為模板地震震級(jí)，Ao為該臺(tái)站記錄到的S波最大振幅，R(do)為臺(tái)站震中距為do時(shí)的量規(guī)函數(shù)，同樣地，自動(dòng)識(shí)別出的地震震級(jí)Mx可寫為

Mx=lg(Ax)+R(dx)

(2)

其中，Ax為與模板同一個(gè)臺(tái)站記錄到的S波最大振幅，dx為自動(dòng)識(shí)別出的地震與該臺(tái)的震中距。兩式相減，并考慮到do≈dx，可得

Mx=lg(Ax/Ao)+Mo

(3)

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

圖 6 研究區(qū)內(nèi)的M-t圖

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

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

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

圖 7 研究區(qū)精定位后的震中分布

圖 8 不同方向的深度剖面

3 結(jié)論

本文設(shè)計(jì)了一套基于模板匹配的震相關(guān)聯(lián)和提取流程，不僅可對(duì)高度匹配的疑似地震進(jìn)行時(shí)空強(qiáng)參數(shù)的測算，也能提供每個(gè)疑似事件對(duì)應(yīng)的震相到時(shí)信息，可供后續(xù)的速度結(jié)構(gòu)成像、臺(tái)站校正等分析使用。

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

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

由龍門山斷裂帶南段地震空區(qū)內(nèi)的模板識(shí)別測試結(jié)果可以看出，本文設(shè)計(jì)的流程可提供包括震源參數(shù)和每個(gè)地震震相到時(shí)信息在內(nèi)的完整的震相報(bào)告，在模板識(shí)別領(lǐng)域可作為一套計(jì)算機(jī)自主處理的可行方案。