蔣一然，寧杰遠,2*，李春來

（1.北京大學地球與空間科學學院，北京 100871；2.河北紅山地球物理國家野外科學觀測研究站，河北 邢臺 054000；3.中國地震局地球物理研究所，北京 100081）

0 引言

模板匹配方法是檢測微小地震、慢滑移等活動的重要方法[1-4]。對于位置接近、震源機制解相似的2 個地震，其到同一臺站的路徑相近，記錄也相似。模板匹配方法就是使用其中一個地震的波形作為模板，與連續(xù)記錄互相關，當互相關值較高時則認為可能發(fā)生了相似地震。單臺的結果容易受隨機噪聲的干擾產生誤識別，因而將多個臺站的互相關結果疊加起來以減少誤識別的比例。同時，由于噪聲的影響，能量較小的震相與模板波形的互相關值雖然會略高于平均水平，但并不顯著。借助多臺疊加的方法，隨機噪聲引起的互相關值波動被平均，來自微小地震的高互相關值因為時間上的一致性被加強，使得小能量的微小地震信號變得顯著。這也使得模板匹配方法可以檢測到單臺上能量與噪音水平相當，甚至略低于噪音水平的微小地震。

與模板位置不同的微小地震的識別具有重要的意義，不僅能夠識別到更大范圍內的地震，還能用于斷層微結構研究。但是，當微小地震與模板位置有一定差距時，即使波形相似度較高，但由于波場到達臺站的走時與模板地震不同，不同臺站的互相關值峰值時間不一致，也不能疊加起來，因此無法被檢測到。一種改進的思路[5]是在將模板地震的位置在原始位置附近進行格點式的挪動，根據(jù)與原始位置的位置差和參考的速度結構對時間移動的大小做微小調整，從而提高模板檢測微震的空間范圍。在這種方法中，格點移動的范圍、移動的步長都將顯著地影響計算量，阻礙高精度定位；參考模型選取得合適與否，也將對檢測能力造成較大影響。另一種思路則是擴寬單臺上高互相關值峰的寬度，使得與模板走時差在一定范圍內的互相關峰值也能夠疊加起來[1,6]。這種擴寬單臺互相關峰值后疊加的方法被稱為弱模板匹配方法（weak matched filter technique，WMFT）。弱模板匹配方法在不顯著增加計算量的情況下提高了單個模板探測微小地震的半徑，因而是一種較好平衡計算代價和檢測能力的微小地震識別方法。

對于模板匹配而言，模板多可以提高檢測微小地震的能力。對于缺乏人工目錄和到時標準的流動臺網，模板相對缺乏。而在有大量地震目錄和波形積累的固定臺網中，人工標注的部分也往往是信噪比較高的部分，對應的地震震級較高。一般來講，震級差距較大的地震，震源時間函數(shù)差距也較大，臺站的波形記錄相似度也較低。震級較高的地震所能檢測到的微小地震的震級也會偏高。近年來，大量以深度學習方法為基礎的地震自動檢測和震相拾取方法被提出[7-11]，大大減少了獲取地震目錄和震相走時的成本，并提高了檢測結果的一致性。借助深度學習的震相自動拾取算法，可以獲得震級下限更小的地震目錄作為模板，從而探測震級更小的微小地震。

由于每個模板都需要和連續(xù)記錄進行互相關，即使是已經大大減少運算代價的弱模板匹配，仍然需要進行大量的運算?；ハ嚓P運算對于每段波形的掃描結果是獨立的，可以很好地并行化。使用torch[12]軟件包實現(xiàn)了弱模板匹配中互相關、歸一化、拓展互相關峰值3 個主要計算部分的GPU 并行，極大地提高了運算速度，使得在相同算力的情況下，可以使用更多的模板，從而檢測出更多的微小地震。

本文將并行優(yōu)化后的弱模板匹配方法應用到川滇地區(qū)，以基于深度學習的APP++震相拾取方法[11]拾取到的地震[13]作為模板，掃描了2019 年6 月17 日四川長寧地震前后各1 個月共60 天的連續(xù)波形，檢測到81 704 個微小地震，并篩選出7 618 個可靠的微小地震。在由雙差地震層析成像[14-15]確定的模板地震位置的基礎上，利用主事件法對模板探測到的微小地震位置進行定位，并分析了地震序列的時空分布特征。

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

1.1 GPU 并行下的弱模板匹配

包括弱模板匹配在內的模板匹配方法面臨的主要困難之一是龐大的計算量，弱模板方法因計算量更大，問題尤為突出。提高模板匹配的計算效率，就可以在相同的時間內，使用更多的模板掃描更多的連續(xù)記錄，從而檢測更多的微小地震，從而得到有意義的研究結果。為了盡可能地提高模板匹配方法的運行效率，本文具體分析了弱模板匹配方法的計算過程和主要的運算代價部分，借助深度學習軟件包torch 中的相關函數(shù)，實現(xiàn)了弱模板計算的GPU 并行化。

弱模板匹配方法的第一步是將模板地震的每一個震相波形與相應臺站的連續(xù)數(shù)據(jù)互相關，得到互相關值的時間序列。具體地講，對于一個模板，其第i個波形記錄wi(t)與相應臺站的聯(lián)系記錄ri(t)做 互相關得到的互相關時間序列corri(t)如下：

式中：t表示時間；dt是時間采樣間隔；Ti為用于互相關的模板波形的起始時刻；N為用到與互相關的采樣點的個數(shù)；To為模板地震的發(fā)震時刻。這里使用了(Ti-To)對不同走時的震相的互相關結果進行了對齊，分母用于將互相關結果歸一化到0～1 之間。

深度學習算法中，一般不關心卷積核中參數(shù)的時序，因而一般用互相關操作來進行等價卷積，以簡化程序的實現(xiàn)。具體到torch 軟件包中，用一維的卷積操作conv1d（實際為互相關）來實現(xiàn)在GPU中的并行計算。值得注意的是，這里的卷積操作并沒有對每一段進行歸一化。因此，作為歸一化因子的分母需要單獨計算。分母的計算由兩部分構成，一部分是模板波形的平方和，另外一部分是每時刻用于互相關的連續(xù)波形的平方和。第一部分對于每個時刻是不變的，計算量較小；第二部分是隨時間變化的，計算量與歸一化因子的分子相當。在此將第二部分的計算轉化為互相關計算：

式中：1(t)表示一個值恒為1 的時間序列，分母第二部分的計算相當于ri的平方與這個序列做互相關，也可以使用conv1d 函數(shù)予以實現(xiàn)。在北京大學實驗室的工作站上，對使用CPU（AMD 2990wx）和GPU（1 080 Ti）計算互相關值的平均用時進行了統(tǒng)計（圖1），結果顯示使用GPU 的平均耗時約為CPU的1/40。需要特別指出的是，如果使用頻率域互相關的方法，運算速度會更快。但是，由于連續(xù)記錄時間長度長、波形能量變化大，采用頻率域計算的方法，會面臨大數(shù)吃小數(shù)的數(shù)值問題，對于互相關值的計算精度影響很大。在測試過程中，conv1d 函數(shù)在CPU 上進行時域互相關時僅調用了2～3 個核心，并行度不高。如果提高CPU 上的并行化程度，能夠在一定程度上提高計算效率，但GPU 仍然有較大的效率優(yōu)勢。

圖1 CPU-GPU 互相關效率測試

在得到歸一化后的互相關值后，弱模板匹配方法要求需要對每一道互相關值的峰值進行拓寬，得到新的互相關值序列。具體地講，根據(jù)直接互相關序列的均值m和標準差s，設定單臺的互相關值閾值threshold為：

式中：mul為峰值需要大于平均值的最低標準差倍數(shù)。當互相關值大于閾值時，視為需要拓寬的峰值。這里倍數(shù)mul的選取，需要保證能夠正確拓寬微小地震的高互相關值，又不過多地拓寬隨機噪聲引起的互相關高值。經過實驗，一般取mul為4。

從實際操作上講，這個拓寬的過程和深度學習中的最大值池化過程maxPooling 類似。最大值池化是選取時間點附近一定時間窗內的最大值代替原有值，即：

式中：twin表示選取的用于池化的時間窗大小，等于需要拓寬峰的時間長短。時窗長度變長可以增加模板的檢測范圍，但同時也會增加誤識別比例，其選取需要平衡這二者。帶閾值的峰的拓展中與最大值池化略有不同的地方在于，當最大值小于閾值時，就不使用最大值對原有的值進行替代。因此，最終得到的拓寬后的互相關值序列corriw(t)為：

這樣，就可以借助torch 軟件包中的maxPooling函數(shù)實現(xiàn)并行的峰的拓寬。結合上面的處理方法，就可以將弱模板主要計算部分全部通過GPU 并行加以實現(xiàn)。在實際掃描過程中，一般先把一天的連續(xù)記錄加載到內存中，然后與多個模板進行互相關。考慮到數(shù)據(jù)在內存和顯存之間交換會耗費相當?shù)臅r間，預先將波形數(shù)據(jù)加載到顯存中。一般來講，因顯存的大小有限，無法同時將一天的連續(xù)數(shù)據(jù)和全部的模板波形加載到顯存中。由于每個臺站的連續(xù)記錄會和很多模板在這個臺站上的波形進行互相關，而每一次互相關中連續(xù)數(shù)據(jù)的數(shù)據(jù)量大于模板波形的數(shù)據(jù)量，所以可以優(yōu)先將連續(xù)記錄加載到顯存中。

之后，將不同模板波形的互相關結果進行平均，得到總的互相關疊加結果corrall(t)：

式中：I表示所有使用到的模板波形數(shù)量。根據(jù)corrall(t)的均值M和標準差S設定檢測微小地震的閾值Threshold為：

式中：Mul是需要給定的微小地震互相關值高于平均值的標準差倍數(shù)，一般取Mul等于6，得到初步的微小地震目錄后再進行分析和質量控制。

1.2 濾波頻段的選取

在進行互相關之前需要對數(shù)據(jù)進行濾波，對微小地震波形主要能量所處頻段外的其他信號進行壓制。在模板匹配中，濾波頻段的選擇還應該考慮使得相似地震的互相關值更容易具有較高的互相關值。有兩個完全相同的波場，都記為wi(t)，兩者之間做互相關。在（1）式的定義下，互相關值在T0處取得最大值1。但是由于離散采樣，實際記錄中兩個波場的采樣點會有一定差異，設其為dt′∈，dt為數(shù)據(jù)的采樣間隔。此時實際只能獲得T0+dt′處的互相關值：

將wi(Ti+n×dt)和wi(Ti+n×dt+dt′)作傅立葉展開：

在dt′的偏移下，同一頻率的信號會有一個相差，越高頻的部分，相差越大，越容易造成互相關值的降低。通過壓制高頻部分的能量，可以在一定程度上減少采樣點差異對互相關值的影響。

本文中使用數(shù)值試驗的方式來測試在采用不同的濾波頻段對互相關值的影響。先生成1 000 Hz采樣的隨機信號當作連續(xù)波場的近似，并將波場降采樣到50 Hz，生成兩段采樣點相差0.01 s 的波形，并且進行0.5～20 Hz、0.5～12 Hz、0.5～8 Hz 三個不同頻段的帶通濾波。將在同一頻段濾波但有采樣點偏差的記錄做互相關（圖2）。結果可以看出，當波形包括較高頻段的能量時，即使是相同連續(xù)波場下采樣的波形互相關值也顯著低于1。隨著帶通濾波最大頻率的下降，互相關值逐漸上升。當取0.5～8 Hz 時，互相關值趨于穩(wěn)定。因此我們選取0.5～8 Hz作為弱模板匹配的互相關頻率。

圖2 采樣點有偏差時相同波場的互相關結果在不同帶通濾波下的變化測試

1.3 微小地震的相對定位

在檢測到微小地震后，在每個臺站的直接互相關結果對應時間點附近選取互相關值較高的點作為微小地震在該臺站上的震相到時，并以疊加的互相關值峰值位置的時刻作為對發(fā)震時刻的初步估計。之后采用主事件法，計算小地震相對于模板地震的位置。用x表示模板地震的位置，x′表示微小地震的位置，ti和τio分別表示模板地震和微小地震的第i個波形的走時。這里 τio是由微小地震初步估計的發(fā)震時刻計算的。當對微小地震的發(fā)震時刻調整 Δτ時，獲得新發(fā)震時刻下的走時 τi：

考慮到2 個地震事件的位置十分接近，于是根據(jù)模板地震走時及其對震源位置的一階導數(shù)來估計微地震的走時作為理論到時 τit：

式中：I表示使用到的震相波形的總數(shù)；wi是根據(jù)每個震相波形的互相關值設定的權重。在此設定為當互相關值大于0.4 時，權重wi為互相關值的平方，否則設定為千分之一。由于模板地震一般具有較高的能量，走時和絕對位置更容易確定。而微小地震能量較低，絕對位置確定較困難。這里依賴于模板地震，計算微小地震相對于模板地震的相對位置，則更為容易。

1.4 模板地震與連續(xù)數(shù)據(jù)

蔣一然等[11]基于U-net++[18]網絡改進了APP 方法，提出了基于臺陣策略APP++震相自動拾取算法，并使用該方法掃描了川滇地區(qū)的固定臺網（116 個）和小江斷裂帶附近的高密度流動臺陣（51 個）2014—2019 年間6 年的連續(xù)波形數(shù)據(jù)（臺站分布如圖3），共檢測到73 291 個地震，537 554 個P 波記錄，471 459 個S 波記錄。在此基礎上，蔣一然等[13]挑選了其中的21 160 個地震，利用雙差地震層析成像方法，獲得了該地區(qū)淺部的三維速度結構和17 240 個地震的精確空間位置。本文使用這17 240個地震作為研究模板地震，并依據(jù)其精定位結果對微小地震進行相對定位。需要特別指出的是，蔣一然等[13]在進行雙差地震層析成像的時候，使用了由互相關提取的走時差作為約束，而本文中對微小地震定位也是使用了互相關提取的走時差，兩者具有一致性，因而不同模板掃描得到的微小地震的位置間也是可比較的。

圖3 APP++方法同時間段內識別到的地震

下面選取與模板地震同一臺陣的2019 年6 月17 日長寧地震前后一個月共60 天的連續(xù)波形，使用GPU 并行下的弱模板匹配技術掃描其中的微小地震。

2 弱模板匹配結果

2.1 初步結果和質量控制

在50 Hz 采樣率下，設定帶通濾波頻帶為0.5～8 Hz，震相波形截取為震相到時前3 s 和后4 s；模板地震至少含有6 個震相波形，如果超過40 個震相波形，則取相應臺站震中距最小的40 個震相波形；對于單臺互相關結果，當某些時刻的互相關值大于整體平均值4 倍標準差時，將其后面0.45 的互相關值全部置為該互相關值；取疊加結果中高于平均值6 倍標準差的峰值為檢測到的微小地震。一共檢測到81 704 個地震，圖4 展示的是其中的2 個例子。

圖4 兩個弱模板匹配圖

進一步統(tǒng)計了相同時間段內，在該區(qū)域由弱模板匹配和APP++方法拾取到的地震震級頻次分布（圖5a）?？梢钥吹剑珹PP++的結果在完備震級之上基本滿足G-R關系[19]，而弱模板匹配的結果則在震級小于2 時，數(shù)量明顯增多。經分析，這是由于設計的閾值過低，導致一部分噪聲被誤識別為地震，從而使得地震目錄的G-R關系被破壞。通過統(tǒng)計不同Mul取值下，弱模板匹配檢測到地震數(shù)目的變化（圖5b），發(fā)現(xiàn)在8～12 之間，地震數(shù)目的變化和Mul取 值變化間滿足很好的指數(shù)關系。而當Mul小于8 時，誤識別比例增加，地震數(shù)目變化不再符合之前的指數(shù)關系。故將Mul設定為8，重新統(tǒng)計了地震震級頻次分布（圖5a），發(fā)現(xiàn)提高閾值后的震級分布基本滿足G-R關系。和APP++方法比較，質量控制后的弱模板匹配方法的地震檢測能力比APP++方法提高約0.5 個震級，數(shù)目約為APP++方法的3 倍（APP++拾取到2 779 個地震分布如圖3，弱模板拾取到7 618 個地震分布如圖6）。

圖5 地震數(shù)目統(tǒng)計圖

圖6 篩選后弱模板匹配檢測到的地震分布圖

2.2 地震序列

圖6 中，地震的震中位置與當?shù)氐臄鄬游恢糜泻芎玫囊恢滦浴D中剖面水平線匯聚的位置為長寧地震主震及其余震發(fā)生的區(qū)域。選擇紅色矩形框內的地震，統(tǒng)計其發(fā)震時刻和震級分布（圖7）。圖7a～7b 中分別展示了APP++方法拾取到的地震序列和弱模板匹配方法拾取到的地震序列。APP++地震序列中，可以看到有4 個左右時間上的叢聚序列（近垂直的細條）；而在弱模板匹配目錄中，借助更多小震發(fā)現(xiàn)較多時間上的叢聚序列。另外，由于弱模板匹配識別到的地震事件較多，對于地震活動性的分辨力也會更高。由弱模板匹配方法的目錄中可以發(fā)現(xiàn)，在長寧地震發(fā)生前有一個較為明顯的平靜期，這可能與地震臨震狀態(tài)的動力學過程有關。借助弱模板匹配方法識別到的地震目錄可以更好地研究地震活動性及其相應的構造過程。

圖7 2019 年6 月17 日長寧地震前后地震震級隨時間的變化

2.3 地震的空間位置分布

以圖6 中的4 條線段為剖面，展示地震在垂直剖面上的位置（圖8）。A—A′剖面120～200 km 之間和C—C′剖面170～250 km 之間，地震叢集大致為南北兩個中心，主要集中在20 km 以內，更深部的地震沒有體現(xiàn)出較為明顯的取向；B—B′的120～200 km和D—D′剖面200～300 km，兩個叢集的垂直投影相互重疊不易分辨，深部的地震體現(xiàn)出一定的線性特征，可能與該位置的斷層結構有關；D—D′剖面150 km處出現(xiàn)的地震叢集體出現(xiàn)明顯的取向，反映地震在相應斷層發(fā)生的位置分布。獲得的微小地震空間位置分布可以進一步用于對斷層幾何形態(tài)、地震活動性及地震活動性變化的研究。

圖8 沿剖面的地震位置分布圖

3 結論

1）將弱模板匹配的互相關、歸一化、拓寬峰值3 個主要計算部分實現(xiàn)了GPU 并行，提高了運算效率。

2）進一步考慮了離散采樣對波形互相關結果的影響，討論不同濾波頻段下互相關結果的變化，并選取了合適的濾波頻段。

3）將GPU 并行下的弱模板匹配方法應用于川滇地區(qū)，掃描了2019 年6 月17 日長寧地震前后各1 個月共60 天的連續(xù)數(shù)據(jù)，并利用地震數(shù)量分布與互相關值的關系對檢測結果進行質量控制，并最終檢測到了基于深度學習方法自動檢測的地震目錄3 倍數(shù)量的地震。

4）微小地震構建了更為細致的地震序列，并能夠反映更多時間上地震的叢聚。

5）借助于經過模板匹配方法精定位過的模板位置信息，對微小地震利用主事件法定位，獲得了微小地震的空間位置分布，并可以進一步用于對斷層幾何形態(tài)、地震活動性及其變化的研究。

致謝中國地震局地震預測研究所為本研究提供地震波形數(shù)據(jù)；研究工作得到北京大學高性能計算校級公共平臺支持。