石永祥，蔣一然，寧杰遠(yuǎn),2*，李春來

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

0 引言

近年來，以機器學(xué)習(xí)與人工智能為基礎(chǔ)的地震自動拾取技術(shù)發(fā)展很快，可以從連續(xù)波形記錄中檢測到大量微震，具有了形成較為完備微震目錄的能力[1-4]。理論上，大部分的微震應(yīng)該分布在斷層面上，通過觀察微震的分布就可以獲得該地區(qū)的斷層空間結(jié)構(gòu)。但是，微震目錄具有微震數(shù)目龐大、信號質(zhì)量參差不齊的特點。同時，一個地區(qū)各種斷層結(jié)構(gòu)產(chǎn)生的微震在空間上會混雜在一起，難以通過人工分割開屬于不同斷層的微震序列。

機器學(xué)習(xí)方法在處理這類大數(shù)據(jù)的問題上具有天然優(yōu)勢，可以在一定程度上規(guī)避掉一些噪聲的影響，提高分選效率。由于獲得的微震目錄是無標(biāo)簽樣本，不能采取SVM、CNN 等目前流行的監(jiān)督類機器學(xué)習(xí)方法，只能采用非監(jiān)督學(xué)習(xí)方法。作為典型的非監(jiān)督方法，基于空間距離的K-Means 聚類方法可以自動分離不同空間位置的地震，因而被用于研究地震的空間分布，并劃分地震活動帶[5-6]。

真實的微震序列大都發(fā)生在斷層面上。因而，同一個斷層產(chǎn)生的微震序列是面狀分布，而不是球形分布。此時，K-Means 方法無法有效地分開同一地區(qū)不同斷層產(chǎn)生的地震序列。非斷層區(qū)域的地震密度一般較低，不同斷層產(chǎn)生的微震之間會存在低密度區(qū)。當(dāng)微震定位足夠精確時，可以通過微震低密度區(qū)來分割不同斷層，從而刻畫出這個地區(qū)的斷層空間分布。

本文使用一種基于密度聚類的方法（Density-Based Spatial Clustering of Applications with Noise,DBSCAN）進(jìn)行微震的聚類分析。在初始的微震目錄中，由于臺站密度有差異，不同區(qū)域的平均微震密度差別很大，難以通過一個標(biāo)準(zhǔn)來分離出不同斷層的微震序列。因此，對微震目錄進(jìn)行初次聚類，首先獲得較小范圍內(nèi)微震分布。之后，對分離后的小區(qū)域內(nèi)的微震再次聚類，從而獲得該區(qū)域微震的精細(xì)分布。為了更加精確地區(qū)分不同斷層產(chǎn)生的地震，在再次聚類時嘗試在聚類中加入波形的相關(guān)信息。在二次聚類之后，得到的微震分布會更加符合當(dāng)?shù)氐臄鄬臃植?，能夠從一些類中發(fā)現(xiàn)可能的斷層結(jié)構(gòu)。

1 DBSCAN 聚類方法

DBSCAN 方法[7-8]是一種經(jīng)典的密度聚類方法。該方法使用2 個控制參量來決定這個屬性：半徑?和最小樣本數(shù)目MinPts。；1 所示。DBSCAN 會對樣本中的每一個樣本進(jìn)行遍歷，當(dāng)一個樣本半徑?內(nèi)存在MinPts個樣本點時，這個樣本點被認(rèn)為處于稠密區(qū)域，滿足稠密條件。這個點與半徑?內(nèi)所有的點組成一個樣本點集合S，之后在這個集合S 內(nèi)尋找所有同樣滿足稠密條件的樣本點，并將它們的集合包含在S 內(nèi)并不斷循環(huán)，直到S“邊緣”的所有樣本點都無法滿足稠密條件，從而認(rèn)為S 中的所有點組成同一個類。這種聚類方法不需要迭代，所需要的運算時間與運算量都較小，能夠滿足對數(shù)萬樣本點的快速聚類需求。

表1 DBSCAN 方法流程偽代碼

DBSCAN 方法中沒有K-Means 中的中心點的概念，只需要度量任意兩點之間的距離即可。此處在初次使用震源位置信息聚類時，使用的距離度量只包含了位置信息：

其中：Eq1 與Eq2 代表不同的地震；X為定位得到的震源位置，單位為km。在二次聚類時，距離度量加入了不同事件在同一臺站的波形的信息：

式中：R1,i為i號臺站對事件1的記錄；R2,i是i號臺站對事件2 的記錄；coor(R1,i,R2,i)是2個記錄之間的點乘，也即零移互相關(guān)函數(shù)，取值歸一化到[0，1]，用于計算2 個記錄之間的差別；NS是臺站數(shù)目；這里的 α權(quán)重在不同小區(qū)域會有所變化，以適應(yīng)不同區(qū)域地震與臺站密度的差異。

2 聚類結(jié)果與分析

本文所用微震目錄是利用APP 方法[4]、從川滇地震臺網(wǎng)連續(xù)數(shù)據(jù)中獲取到P 波與S 波的到時之后，使用IASP91 地球模型對地震位置進(jìn)行定位得到的結(jié)果[9]。本文選用這個目錄中60 000 個地震記錄進(jìn)行分析，先使用地震的三維位置信息進(jìn)行粗略的聚類，之后再對幾個大類加入波形信息后進(jìn)行更加細(xì)致的分析。

2.1 使用位置信息的初次聚類

在第一次聚類中，DBSCAN 方法的2 個控制參量被設(shè)置為半徑?=10 km和MinPts=20，生成了震事件，微震分布與當(dāng)?shù)氐臄鄬幼呦蛴幸欢P(guān)系。88 個類（包括一個包含9 778 個事件的“噪點”類），結(jié)果如圖1 所示。在圖1 中，不同的類使用不同的顏色標(biāo)注，K=-1 對應(yīng)黑色的點為噪點，彩色的是有效類，并按照數(shù)目排序。聚類的結(jié)果中包含最多地震的類是在（26.5°N，102.8°E）附近（圖1 粉色橢圓內(nèi)），包含了23 575 個地震。這是因為這個區(qū)域加入了一個密集臺陣的數(shù)據(jù)，所以記錄到的微震較多。在圖2a 中展示了這一類的微震分布，同時在圖2b 與圖2c 中展示了AA′與BB′兩個切線上微震的深度分布，范圍為線兩側(cè)10 km。第二類（圖1 綠色橢圓）對應(yīng)在龍門山斷層（32°N，104°E）附近NESW 走向的微震序列，包含了8 856 個地震，與該地區(qū)主要的斷層走向相符（圖3）。第三類在（28°N，101.3°E）附近（圖1 淺藍(lán)色橢圓），包含了2 231 個地

圖1 利用位置信息進(jìn)行初次聚類分析的結(jié)果

圖2 初次聚類分析后微震最多的類

圖3 初次聚類分析后龍門山附近一類

圖1 結(jié)果中分離了9 778 個噪點類微震，分布見圖4。噪點類幾乎均勻分布在所有的地區(qū)，沒有明顯規(guī)律可循。噪點類的成因較多，可能與微震目錄中的誤識別或者定位過程中由于信號較弱導(dǎo)致的偏差有關(guān)。噪點類可以使用更低的稠密性條件進(jìn)行進(jìn)一步的聚類，從中抽取出更多的信息，但是總體利用價值不是很大。而對于圖1 中的第一類和第二類，由于區(qū)域較廣，同時地震密度大，無法從這次聚類的結(jié)果中直接觀察到斷層的信息，需要進(jìn)一步聚類分析。類似地，對于（26.3°N，100°E）附近的一系列類（圖1 深藍(lán)色橢圓內(nèi)），可以采用同樣的二次聚類方法從中抽取出更加精確的微震聚類信息，從而反映出更精細(xì)的斷層位置。本文將以第一類與第二類的二次聚類結(jié)果為例進(jìn)行研究。

圖4 初次聚類分析后噪點類

2.2 包含P 波和S 波互相關(guān)信息的二次聚類

在二次聚類分析中，先考慮對初次聚類中數(shù)目最大、微震最為密集的第一類（圖2）進(jìn)行分析。為了增強不同地震事件的區(qū)分度，除了考慮對應(yīng)的位置信息外，還加入了不同事件之間波形的互相關(guān)信息。由于是微震，密集臺陣中只有少數(shù)臺站可以記錄到，而不同臺站記錄的不同事件的波形，有差異較大的方位角與距離，無法比較，故加入限制條件。只有當(dāng)2 個事件相距10 km 以內(nèi)，且有相同臺站記錄到這2 個事件的時候，才將互相關(guān)的信息加入到2 個事件的“距離”度量上，這樣就可以最大限度地避免臺站與事件之間的距離與方位角因素的影響。用公式表示為：

由于這一類包含的地震事件太多，所以沒有使用地震記錄的全波形做互相關(guān)，而是只使用P 波和S 波前后2 s 的記錄做零移互相關(guān)，以節(jié)省運算時間，提高運算效率?；ハ嚓P(guān)函數(shù)選擇所有同時記錄到這2 個事件的臺站，在相同臺站上對不同事件的P 波與S 波做互相關(guān)，之后疊加所有臺站互相關(guān)值，從而更好地衡量波形之間的相似度。

考慮到這次研究的區(qū)域變小，地震數(shù)目也下降，故選擇聚類參數(shù)半徑?=5 km 和MinPts=10，結(jié)果如圖5 所示。DBSCAN 方法一共生成了55 個類，其中黑色點是噪點類。聚類結(jié)果的第一類（圖5綠色橢圓）包含了13 236 個地震（圖6）。從圖6 可以發(fā)現(xiàn)，第一類在南端有很清晰的沿斷層走向的微震分布，并在（26.3°N,103°E）處發(fā)生了轉(zhuǎn)彎。這與當(dāng)?shù)氐臉?gòu)造聯(lián)系密切，預(yù)示著這個區(qū)域可能存在一個不在本文采用的斷層數(shù)據(jù)集內(nèi)的斷層。在圖5 中粉色橢圓對應(yīng)的一類中，也可以看到一個類似的地震事件分布，展示了較為清晰的傾角與走向（圖7）。

圖5 加入P 波和S 波互相關(guān)信息后二次聚類的結(jié)果

圖6 二次聚類（圖5）中綠色橢圓對應(yīng)的類

圖7 二次聚類（圖5）中的粉色橢圓對應(yīng)的類

從本次分析可以發(fā)現(xiàn)，二次聚類的結(jié)果相對于沒有聚類的結(jié)果，明顯帶有了更多的斷層信息，因而仍需對數(shù)目較多的類進(jìn)行第三次聚類，進(jìn)一步分析斷層信息。

2.3 包含全波形互相關(guān)信息的二次聚類

對圖1 中在龍門山附近區(qū)域的第二類（圖4）進(jìn)行二次聚類分析。與圖3 中的第一類不同，這個類只包含8 856 個地震，互相關(guān)所需要的運算量較少。因而采用P 波到達(dá)前10 s 到S 波后10 s 的波形進(jìn)行零移互相關(guān)，之后將互相關(guān)值加入到事件之間的距離度量上，設(shè)定

聚類參數(shù)設(shè)置為?=5 km和MinPts=10，DBSCAN 對龍門山地震序列附近再次聚類的結(jié)果，生成了18 個地震類（圖8）。其中數(shù)目最多的第一類（圖9）包含4 947 個地震。與第一次聚類得到的結(jié)果（圖4）相比，此時微震沿斷層分布的特征更加清晰。第一類的分布對應(yīng)著龍門山斷裂的一段，南端的NW 向翹起則對應(yīng)此處的NW 向分支斷層。

圖8 龍門山斷層附近一類二次聚類的結(jié)果

進(jìn)一步分析聚類結(jié)果。從圖9 事件的深度以及在經(jīng)度與緯度上限的分布中可以看到，在（31.6°N，103.8°E）附近地震震源深度大多大于10 km，淺層地震很少，而兩端存在著大量的淺源地震，形成了明顯的地震空區(qū)。這種分布特征與“5·12”汶川地震余震的空區(qū)分布十分相似[10]，兩端大量的淺源地震主要與龍門山斷裂在兩端發(fā)展出的大量不同走向的分支斷層有關(guān)（如該類南端的NW 向分支）。而中間分支斷層較少，地震主要發(fā)生在龍門山主斷裂上，震源的位置較深，形成了這個地震空區(qū)。同時，圖8 分類結(jié)果中龍門山斷裂的北端（104.8°E 以東）被劃分成了單獨的一類，這可能是因為北端處于斷裂帶的末端，發(fā)震結(jié)構(gòu)、活動性與主斷裂不相同，因而被分離出。整個分段結(jié)果與Tian et al.[11]的研究結(jié)果形成了很好的對應(yīng)關(guān)系，從地震空間分布特征的角度驗證了Tian et al.關(guān)于震源機制解分區(qū)分布的結(jié)果。尤其是清晰地從地震空間分布的角度驗證了關(guān)于存在兩種逆斷層地震及其空間分布的結(jié)果，并從地震空間分布的角度表明了第二類逆斷層地震的活動范圍，證實了Tian et al.關(guān)于現(xiàn)存斷層釋放應(yīng)力場不同分量的結(jié)論。

圖9 龍門山斷層附近二次聚類后的第一類

在聚類結(jié)果中，還在斷裂帶的南端分離出了一些其他的類，反映了當(dāng)?shù)貜?fù)雜的斷層分支的存在。原可以進(jìn)一步分離出其中多個斷層，但是受限于定位的精度與地震數(shù)目，進(jìn)一步的聚類存在困難。

3 結(jié)論

通過使用基于密度的聚類方法分析了自動拾取技術(shù)獲得的川滇地區(qū)微震目錄，并進(jìn)一步地利用空間位置信息分離出了不同地震活動帶內(nèi)的微震序列。對于使用APP 方法定位出的60 000 個地震，先使用位置信息進(jìn)行初步聚類，之后再對數(shù)目較多的類調(diào)整聚類參數(shù)后進(jìn)行二次聚類；二次聚類中同時加入波形互相關(guān)信息，以加強不同斷層類之間地震事件的區(qū)分度。一次聚類的結(jié)果大概展示了該地區(qū)微震分布的區(qū)域性，但是受到臺站密度等多方面的影響，并不能給出相對清晰的斷層結(jié)構(gòu)。對一次聚類結(jié)果進(jìn)行二次聚類，則更好地展示出了地震序列沿斷層的分布，并給出了一些斷層所對應(yīng)的地震序列。

同時，二次聚類后依然存在一些包含大量微震的類可以進(jìn)一步分離。但是想要獲得更可靠的結(jié)果，需要掃描更長時間的連續(xù)記錄，獲得更多的微震序列。并且，也可以使用Tomo-DD 等反演方法得到精確的地下結(jié)構(gòu)對地震進(jìn)行定位，從而獲取更加精確的地震分布信息，提高聚類方法結(jié)果的精確度與可靠程度。

初步結(jié)果表明，基于密度的聚類方法具有推廣價值。

致謝中國地震局地震預(yù)測研究所為本研究提供地震波形數(shù)據(jù);研究工作得到北京大學(xué)高性能計算校級公共平臺支持；感謝參與討論并給出建議的溫景充、鮑鐵釗、伍晗與殷常陽同學(xué)。