基于改進(jìn)的DBSCAN算法自動(dòng)識(shí)別斷層方法研究及其在唐山地區(qū)的應(yīng)用

張?zhí)K祥 盛書(shū)中* 席 彪 房立華 呂 堅(jiān) 王甘嬌 張 瀟

1)東華理工大學(xué)，地球物理與測(cè)控技術(shù)學(xué)院，南昌 330013 2)長(zhǎng)江大學(xué)，油氣資源勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430100 3)中國(guó)地震局地球物理研究所，北京 100081 4)江西省地震局，南昌 330039

0 引言

近年來(lái)，隨著地震觀測(cè)臺(tái)網(wǎng)的不斷加密和地震儀觀測(cè)能力的提高，觀測(cè)到的地震事件數(shù)劇增，定位的準(zhǔn)確性也不斷增加，為精細(xì)刻畫(huà)斷層構(gòu)造提供了基礎(chǔ)。大量基于震源空間位置分布獲取深部斷層構(gòu)造的研究工作陸續(xù)發(fā)表(Ouillonetal.，2008； 萬(wàn)永革等，2008； Hayesetal.，2010； Wangetal.，2013； 盛書(shū)中等，2014，2015； Skoumaletal.，2019； Kameretal.，2020； Brunsviketal.，2021)。這些研究表明，從海量的地震數(shù)據(jù)中獲取斷層形態(tài)及其參數(shù)是發(fā)震構(gòu)造研究的一個(gè)重要方法。當(dāng)前基于地震資料獲取斷層及其參數(shù)的研究中，選取數(shù)據(jù)的方法主要有2類(lèi)： 1)基于對(duì)斷裂構(gòu)造的認(rèn)識(shí)以及地震數(shù)據(jù)空間分布情況經(jīng)驗(yàn)性地選取地震數(shù)據(jù)，再由這些數(shù)據(jù)擬合出斷層面(萬(wàn)永革等，2008； Hayesetal.，2010； 盛書(shū)中等，2014； 高彬等，2016； 胡曉輝等，2019)。相關(guān)成果極大地推動(dòng)了由余震空間分布確定斷層參數(shù)的研究，但它們依賴(lài)于先驗(yàn)信息——須了解已有斷裂的構(gòu)造，或?qū)Φ卣鹂臻g分布的線性情況有所要求，因此這類(lèi)方法難以處理線性趨勢(shì)相對(duì)較差的地區(qū)。2)基于地震數(shù)據(jù)的空間叢集性，采用機(jī)器學(xué)習(xí)方法中的無(wú)監(jiān)督聚類(lèi)技術(shù)選取數(shù)據(jù)，該方法避免了對(duì)經(jīng)驗(yàn)的依賴(lài)，且更加適應(yīng)從海量地震數(shù)據(jù)中自動(dòng)獲取斷層。如Ouillon等(2008)基于K-means聚類(lèi)法提出的各向異性動(dòng)態(tài)聚類(lèi)方法(OADC)，該方法改進(jìn)了K-means算法，使其不需要輸入K值，并將其應(yīng)用于美國(guó)蘭德斯地震序列，獲得了該地區(qū)的已知斷層和隱伏斷層。但該方法在選取數(shù)據(jù)時(shí)，默認(rèn)將所有地震事件均分配給斷層面，忽略了部分離散事件對(duì)斷層參數(shù)擬合的影響。Skoumal等(2019)基于DBSCAN(Density Based Spatial Clustering of Application with Noise)算法提出了“斷層識(shí)別”方法(FaultID)，并據(jù)此識(shí)別出美國(guó)俄克拉荷馬州約2500個(gè)斷層，其中多數(shù)斷層為先前未被研究的斷層，但該方法在選取數(shù)據(jù)時(shí)需經(jīng)驗(yàn)性指定DBSCAN的2個(gè)全局參數(shù)。Kamer等(2020)提出“生長(zhǎng)聚類(lèi)”算法，選取數(shù)據(jù)并給出斷層結(jié)構(gòu)，由于該方法采用“自下而上”的模式，與其他方法相比對(duì)小規(guī)模、 復(fù)雜的斷層結(jié)構(gòu)的識(shí)別能力更強(qiáng)，但該方法仍需嘗試不同的參數(shù)，以便獲得最優(yōu)結(jié)果。Brunsvik等(2021)首先采用譜聚類(lèi)對(duì)地震數(shù)據(jù)進(jìn)行初步挑選，其次采用DBSCAN選取Paganica斷裂系統(tǒng)中各斷層的數(shù)據(jù)，并最終給出L’Aquila地震序列相關(guān)斷層的三維形態(tài)，但該方法在譜聚類(lèi)以及DBSCAN聚類(lèi)過(guò)程中仍然需要經(jīng)驗(yàn)性地設(shè)置參數(shù)。因此，本研究基于以上方法，進(jìn)一步提出改進(jìn)的DBSCAN聚類(lèi)算法，以自動(dòng)識(shí)別斷層，為快速、 自動(dòng)反演斷層的三維形態(tài)等提供新的思路和方法。

本研究改進(jìn)了無(wú)監(jiān)督聚類(lèi)技術(shù)，實(shí)現(xiàn)了基于震源空間位置分布自動(dòng)識(shí)別斷層段并計(jì)算其參數(shù)。首先應(yīng)用改進(jìn)的DBSCAN聚類(lèi)算法自動(dòng)識(shí)別斷層段，其次對(duì)識(shí)別出的斷層段采用模擬退火全局搜索-高斯牛頓局部搜索結(jié)合法計(jì)算其斷層參數(shù)，再對(duì)鄰近的相似斷層段進(jìn)行合并，最終給出基于地震空間分布識(shí)別出的各斷層段及其參數(shù)。為了驗(yàn)證該方法的可行性，我們首先將其應(yīng)用于合成數(shù)據(jù)，其次將其應(yīng)用于構(gòu)造研究較為深入的唐山地區(qū)。

1 方法

本研究的目標(biāo)是基于震源空間位置分布自動(dòng)識(shí)別震源所刻畫(huà)的斷裂段，因此，提出改進(jìn)的DBSCAN聚類(lèi)算法以實(shí)現(xiàn)自動(dòng)識(shí)別斷層段。鑒于以下2方面考慮： 1)不同斷層以及同一條斷層不同段上的地震密集程度不同； 2)斷層具有任意形狀，且噪聲將影響斷層參數(shù)的擬合。我們選擇對(duì)可以發(fā)現(xiàn)任意形狀聚類(lèi)、 能挖掘數(shù)據(jù)中密度變化且對(duì)噪聲數(shù)據(jù)不敏感的考慮噪聲的空間密度聚類(lèi)算法(DBSCAN)進(jìn)行改進(jìn)(Skoumaletal.，2019； Brunsviketal.，2021)，以實(shí)現(xiàn)斷層段的自動(dòng)識(shí)別。

傳統(tǒng)的DBSCAN算法由Ester等(1996)提出，該算法認(rèn)為斷層上的地震事件由“核心震源”和“邊緣震源”組成，其余地震事件為噪聲(“核心震源”為在給定半徑eps的鄰域內(nèi)震源數(shù)量不小于minPts的震源； “邊緣震源”為在eps鄰域內(nèi)震源數(shù)量小于minPts，但其自身處于核心震源的鄰域內(nèi)的震源)。在DBSCAN算法中，聚類(lèi)是“密度相連”的集合，能夠?qū)⒆銐蚋呙芏鹊臄?shù)據(jù)聯(lián)系在一起，并且在具有噪聲的數(shù)據(jù)中能夠發(fā)現(xiàn)任意形狀的簇。前人經(jīng)研究發(fā)現(xiàn)，傳統(tǒng)的DBSCAN算法存在以下2方面缺陷： 1)聚類(lèi)結(jié)果對(duì)全局參數(shù)eps和minPts比較敏感，如增大minPts會(huì)增加將斷層上的地震事件誤認(rèn)為是噪聲的概率，增大eps會(huì)增加將相鄰的斷層合并在一起的趨勢(shì)及將非斷層上的地震事件誤認(rèn)為是斷層上的地震事件(Brunsviketal.，2021)； 2)單一的eps和minPts值無(wú)法同時(shí)識(shí)別活躍程度差異較大的斷層段(侯雄文，2017a)?？紤]到實(shí)際斷層的復(fù)雜性，固定輸入?yún)?shù)將導(dǎo)致在識(shí)別過(guò)程中丟失大量斷層。因此，本研究針對(duì)以上缺陷對(duì)DBSCAN算法進(jìn)行改進(jìn)，使其能夠在不同密度層次(侯雄文，2017a)的數(shù)據(jù)中自動(dòng)識(shí)別斷層段。

為了解決以上問(wèn)題，本文在付澤強(qiáng)等(2018)和李文杰等(2019)研究的基礎(chǔ)上，進(jìn)一步改進(jìn)DBSCAN算法。李文杰等(2019)提出利用K-平均最近鄰算法(K-Average Nearest Neighbor，K-ANN)和數(shù)學(xué)期望法生成eps和minPts閾值參數(shù)的候選集合，再基于密度層次穩(wěn)定的情況選取最優(yōu)參數(shù)(侯雄文，2017b)。K-平均最近鄰算法是K近鄰法(K-Nearest Neighbor，KNN)和平均最近鄰法(Guyonetal.，2002； Changetal.，2011)的延伸，該算法的基本思想是通過(guò)計(jì)算地震目錄中每個(gè)地震事件與其第K個(gè)最近鄰事件之間的K-最近鄰距離，并對(duì)所有數(shù)據(jù)點(diǎn)對(duì)的K-最近鄰距離求平均值，得到數(shù)據(jù)集的K-平均最近鄰距離。之后，對(duì)所有K值進(jìn)行計(jì)算，得到K-平均最近鄰距離向量，即eps參數(shù)候選集合。

對(duì)于給定的eps參數(shù)候選集合，依次求出每個(gè)eps參數(shù)對(duì)應(yīng)鄰域的地震事件數(shù)量，并計(jì)算所有對(duì)象的eps鄰域數(shù)量的數(shù)學(xué)期望值，作為minPts參數(shù)候選集合。計(jì)算地震目錄D的minPts參數(shù)候選集合的具體公式為

(1)

其中，Pi為第i個(gè)對(duì)象的eps鄰域?qū)ο髷?shù)量，n為地震目錄D中的對(duì)象總數(shù)。

考慮到獲得可靠的斷層參數(shù)所需的地震數(shù)量不能過(guò)少，本研究采用來(lái)自于數(shù)據(jù)自身特征的簇內(nèi)最小事件數(shù)(M)作為斷層段識(shí)別的終止條件，并將其稱(chēng)為有效事件數(shù)閾值(M＿effect)，下面簡(jiǎn)述該值的確定過(guò)程。我們用矩形框?qū)Φ卣饠?shù)據(jù)中呈條帶狀且較為稀疏的區(qū)域進(jìn)行框選，統(tǒng)計(jì)框內(nèi)發(fā)生的地震事件數(shù)M。 考慮到M值的確定存在較大的主觀因素，該數(shù)值精度較低，因此，本研究中將M向下取整到十位數(shù)得到M＿effect。若算法所識(shí)別斷層內(nèi)的地震事件數(shù)不小于M＿effect，則稱(chēng)其為有效斷層。

為了盡可能地識(shí)別斷層，本研究采用由高到低的逐層密度聚類(lèi)法。斷層上的地震密度差異很大，使得在斷層識(shí)別時(shí)，地震稀疏段相對(duì)于地震密集段更容易被誤識(shí)為“噪聲”。因此，在本研究中，為了解決地震事件空間密度差異帶來(lái)的問(wèn)題并盡可能地識(shí)別出研究區(qū)的斷層，我們從高密度層次到低密度層次對(duì)地震目錄進(jìn)行多次不同密度層次的聚類(lèi)(付澤強(qiáng)等，2018)。在對(duì)下一密度層次斷層聚類(lèi)前，將上一密度層次識(shí)別的有效斷層數(shù)據(jù)刪除，保證高密度斷層的數(shù)據(jù)不影響低密度斷層的識(shí)別(Skoumaletal.，2019)。下面介紹各密度層次最優(yōu)參數(shù)eps和minPts的確定方法。從小到大選用不同K值(K=1，2，…，n)所對(duì)應(yīng)的eps和minPts，對(duì)地震目錄進(jìn)行聚類(lèi)，分別得到在不同K值下生成的有效斷層數(shù)目。在上述過(guò)程中，當(dāng)有效斷層數(shù)目連續(xù)3次出現(xiàn)相同值時(shí)，定義其密度層次達(dá)到穩(wěn)定，記該有效斷層數(shù)目為最優(yōu)數(shù)目N。 直到某K值所對(duì)應(yīng)的有效斷層數(shù)目不再為N，此時(shí)密度層次發(fā)生改變。將有效斷層數(shù)目N對(duì)應(yīng)的最大K值視為最優(yōu)K值，該K值對(duì)應(yīng)的eps和minPts為當(dāng)前密度層次的最優(yōu)參數(shù)(李文杰等，2019)。

基于地震叢集性地發(fā)生在斷層面上的研究，萬(wàn)永革等(2008)提出用模擬退火全局搜索-高斯牛頓局部搜索相結(jié)合的方法擬合斷層面參數(shù)。該方法將全局搜索和局部搜索相結(jié)合搜索一個(gè)平面，使得該平面到叢集性小地震的距離最小。該方法可以有效地避免對(duì)初始解的依賴(lài)，從而獲得全局最優(yōu)解。本研究將基于上述斷層段自動(dòng)識(shí)別結(jié)果，應(yīng)用萬(wàn)永革等(2008)的方法計(jì)算各段斷層面的參數(shù)及其誤差。

考慮到同一斷層不同段的活躍程度不同，即密度層次不同，在自動(dòng)識(shí)別斷層段的過(guò)程中，同一斷層可能被識(shí)別為2段。因此，在本研究中，將2個(gè)斷層面參數(shù)相近且緊鄰的斷層段合成為一段，并計(jì)算最終的斷層面參數(shù)。

上文簡(jiǎn)述了本研究提出的基于改進(jìn)的DBSCAN聚類(lèi)算法自動(dòng)識(shí)別斷層的方法，其具體實(shí)現(xiàn)步驟如下：

步驟1： 根據(jù)研究數(shù)據(jù)確定M＿effect。

步驟2： 計(jì)算地震目錄D的距離矩陣Dn×n。本研究采用半正矢公式(Markouetal.，2010)計(jì)算距離，具體公式為

(2)

其中，R為地球半徑；φ1、φ2表示2個(gè)地震事件的經(jīng)度；λ1、λ2表示2個(gè)地震事件的緯度；d為半正矢距離。

Dn×n={d(i，j)|1≤i≤n，1≤j≤n}

(3)

其中，Dn×n為n×n的實(shí)對(duì)稱(chēng)矩陣；n為地震目錄D中地震事件數(shù)；d(i，j)為地震目錄D中第i和j個(gè)事件間的半正矢距離。

步驟3： 對(duì)距離矩陣Dn×n的每行元素進(jìn)行升序排序，則第1列元素所組成的距離向量D0表示對(duì)象到自身的距離全為0。第K列元素構(gòu)成所有地震事件的K-最近鄰距離向量DK。

(4)

步驟5： 根據(jù)式(1)計(jì)算Deps對(duì)應(yīng)的minPts，得到候選minPts集合DminPts。

步驟6： 依次用不同K值(K=1，2，…，n)所對(duì)應(yīng)的參數(shù)對(duì)地震目錄進(jìn)行聚類(lèi)，得到不同K值聚類(lèi)的有效斷層數(shù)目。當(dāng)連續(xù)3次K值所對(duì)應(yīng)的有效斷層數(shù)目為N時(shí)，稱(chēng)N為最優(yōu)數(shù)目。

步驟7： 繼續(xù)執(zhí)行步驟6，直至有效斷層數(shù)目不再為N，則有效斷層數(shù)目為N時(shí)的最大K值為最優(yōu)K值。最優(yōu)K值對(duì)應(yīng)的eps和minPts則為最優(yōu)參數(shù)，基于該最優(yōu)參數(shù)進(jìn)行斷層段識(shí)別。

步驟8： 對(duì)識(shí)別出的有效斷層數(shù)據(jù)進(jìn)行斷層面參數(shù)計(jì)算，并將其從地震目錄D中刪除，對(duì)剩余地震數(shù)據(jù)執(zhí)行步驟2—8，重復(fù)上述步驟直至其識(shí)別斷層內(nèi)的地震事件數(shù)均小于M＿effect時(shí)，斷層識(shí)別終止。

步驟9： 對(duì)相鄰且參數(shù)相近的斷層段進(jìn)行合并，計(jì)算斷層參數(shù)。

2 理論測(cè)試

2.1 生成合成數(shù)據(jù)集

為了檢驗(yàn)本研究提出的方法，我們生成密集分布于斷層面上的地震事件和遍布整個(gè)空間的隨機(jī)地震事件(即噪聲)，由地震事件和噪聲事件組成合成數(shù)據(jù)集。首先，我們?cè)诳臻g{0≤x≤1，0≤y≤1，-10≤z≤0}(x、y和z分別表示經(jīng)度(°)、 緯度(°)和深度(km))上生成分布于3個(gè)斷層面上的地震數(shù)據(jù)，且假定震源點(diǎn)圍繞3個(gè)子斷層面服從正態(tài)分布，其分布區(qū)間為(-0.02°，0.02°)。斷層面方程分別為斷層1：x=0.5，斷層2：y=0.3和斷層3：y=0.7。斷層1、 斷層2和斷層3上的地震事件數(shù)分別為800、 1000和1200個(gè)(圖1a)。其次，在研究區(qū)域內(nèi)添加隨機(jī)地震事件，添加的地震數(shù)分別為總地震事件數(shù)的5%、 10%和20%(圖1b—d)。

圖 1 a 人工合成地震數(shù)據(jù)； b—d 分別添加總地震事件5%、 10%和20%隨機(jī)地震后的地震空間分布圖Fig. 1 The spatial distribution of synthetic seismic data(a) and adding random earthquakes which account for 5%，10% and 20% of total seismic events(b—d)，respectively.黑色圓圈表示隨機(jī)地震(即噪聲)，其余為斷層上的地震事件

2.2 自動(dòng)識(shí)別斷層面方法的測(cè)試

下面我們使用合成數(shù)據(jù)對(duì)上文所述的斷層面識(shí)別方法進(jìn)行測(cè)試，測(cè)試過(guò)程以添加5%噪聲事件的合成數(shù)據(jù)為例。首先，選擇其數(shù)據(jù)中呈條帶狀且較為稀疏的矩形區(qū)域計(jì)算M＿effect，得到的M＿effect為80(圖4a)。其次，按照上文所述步驟和方法進(jìn)行第1密度層次聚類(lèi)。由生成聚類(lèi)K與eps和minPts的關(guān)系圖(圖2a，b)以及K與有效斷層數(shù)目的關(guān)系(圖3a)確定最優(yōu)K值為12，其對(duì)應(yīng)的最優(yōu)聚類(lèi)參數(shù)eps和minPts為2.23和18(圖2a，b)?；谝陨献顑?yōu)參數(shù)進(jìn)行聚類(lèi)，獲得6個(gè)有效斷層和8個(gè)無(wú)效斷層。將已識(shí)別的有效斷層數(shù)據(jù)刪除，對(duì)剩余數(shù)據(jù)重復(fù)上述步驟，進(jìn)行第2密度層次聚類(lèi)。由圖2c、 2d和3b可見(jiàn)，聚類(lèi)最優(yōu)K值為19，其所對(duì)應(yīng)的最優(yōu)聚類(lèi)參數(shù)eps和minPts為6.39和45?；谝陨献顑?yōu)參數(shù)進(jìn)行聚類(lèi)，獲得4個(gè)有效斷層和1個(gè)無(wú)效斷層。再將已識(shí)別的有效斷層數(shù)據(jù)刪除，對(duì)剩余數(shù)據(jù)進(jìn)行第3密度層次聚類(lèi)，本次聚類(lèi)的最優(yōu)K值為26，所識(shí)別簇內(nèi)事件數(shù)均小于M＿effect，停止聚類(lèi)。

圖 2 添加5%噪聲的合成斷層數(shù)據(jù)第1次聚類(lèi)(a、 b)和第2次聚類(lèi)(c、 d)的K-eps、 K-minPts關(guān)系圖Fig. 2 K-eps，K-minPts diagrams of first clustering(a，b) and second clustering(c，d) for synthetic data with 5% noise.標(biāo)記處為最優(yōu)參數(shù)

對(duì)合成數(shù)據(jù)集進(jìn)行了3個(gè)密度層次聚類(lèi)，自動(dòng)識(shí)別出19個(gè)聚類(lèi)，其中有效聚類(lèi)為10個(gè)(圖4a)。圖4a 中4和5斷層段，7、 8和9斷層段的走向相近且相鄰，因此，將其分別合并為同一斷層段。1、 3、 6和10斷層段走向相近但不相鄰，考慮到識(shí)別共軛斷層處時(shí)往往將同一斷層劃分為多段，因此將其合并為同一斷層段。最終在研究區(qū)獲得3個(gè)斷層面(圖4b)。重復(fù)上述方法分別處理噪聲水平為10%和20%的合成數(shù)據(jù)，最終識(shí)別出的斷層面情況見(jiàn)圖4c、 4d。

圖 3 添加5%噪聲的合成數(shù)據(jù)第1次聚類(lèi)(a)和第2次聚類(lèi)(b)的K值與有效斷層數(shù)目關(guān)系圖Fig. 3 K-number of effective faults diagram of first clustering(a) and second clustering(b) for synthetic data with 5% noise.

圖 4 添加噪聲數(shù)據(jù)自動(dòng)識(shí)別斷層面的結(jié)果Fig. 4 Automatic fault plane recognition results for synthetic data.a 添加5%噪聲的合成數(shù)據(jù)未合并的識(shí)別結(jié)果，矩形框?yàn)樽钚∈录?shù)閾值選取區(qū)域； b 添加5%噪聲的合成數(shù)據(jù)斷層面合并后的識(shí)別結(jié)果； c 添加10%噪聲的合成數(shù)據(jù)斷層面合并后的識(shí)別結(jié)果； d 添加20%噪聲的合成數(shù)據(jù)斷層面合并后的識(shí)別結(jié)果； 用不同顏色標(biāo)示識(shí)別出的斷層段

2.3 測(cè)試結(jié)果分析

基于改進(jìn)的DBSCAN算法，對(duì)添加不同噪聲水平的合成數(shù)據(jù)進(jìn)行識(shí)別，得到的斷層結(jié)果如圖 4 所示。由圖 4 可見(jiàn)，本文的算法可以較好地識(shí)別出理論斷層，隨著噪聲的增加，斷層面識(shí)別的效果稍有降低。

由于假定的斷層面是共軛斷層，斷層上的地震分布存在交集，而本研究提出的算法將交集部分?jǐn)?shù)據(jù)歸于某條斷層上，會(huì)出現(xiàn)將共軛斷層中的一條斷層識(shí)別為不連續(xù)的斷層段的現(xiàn)象(圖 4)，這種不連續(xù)總體上對(duì)斷層認(rèn)識(shí)的影響較小。

3 在唐山地區(qū)的應(yīng)用

3.1 數(shù)據(jù)

改進(jìn)的DBSCAN算法理論測(cè)試結(jié)果較好。下面我們將用地震活躍且研究程度較高的唐山地區(qū)的地震數(shù)據(jù)對(duì)其做進(jìn)一步測(cè)試。研究區(qū)的范圍為(39°～40°N，117.2°～119.2°E)，數(shù)據(jù)為經(jīng)雙差定位后的4256次地震事件，時(shí)間范圍為2000年1月1日—2011年4月3日，震源深度分布在0～20km，震級(jí)主要為ML<2.7(圖 5)。地震分布明顯呈條帶狀，有利于本文斷層識(shí)別方法的應(yīng)用。

圖 5 唐山地震序列的震中(a)、 深度(b)和震級(jí)(c)分布圖Fig. 5 Distribution of epicenter(a)，depth(b)and magnitude(c) of the Tangshan earthquake sequence，China.矩形框?yàn)樽钚∈录?shù)閾值的選取區(qū)域

3.2 斷層自動(dòng)識(shí)別結(jié)果

重復(fù)與理論測(cè)試部分相同的數(shù)據(jù)處理步驟，自動(dòng)識(shí)別出唐山地區(qū)的斷層段。基于圖 5 框內(nèi)的97個(gè)地震事件數(shù)據(jù)，確定M＿effect為90。其次，對(duì)數(shù)據(jù)進(jìn)行不同密度層次聚類(lèi)，基于上述終止原則進(jìn)行了4次聚類(lèi)，其最優(yōu)K值分別為10、 4、 5和14(圖6a—d)； 最優(yōu)eps參數(shù)分別為1.6、 1.9、 2.7和5.5； 最優(yōu)minPts參數(shù)分別為59、 15、 14和27。最終自動(dòng)識(shí)別出不同密度層次的37個(gè)聚類(lèi)(圖7a—d)，滿足條件的有效斷層有11個(gè)。

通過(guò)對(duì)11段斷層面參數(shù)進(jìn)行對(duì)比分析及合并，最終獲得8個(gè)斷層(圖 8，表1)： 第1段和第4段(圖7a)、 第4段和第11段(圖7b)的走向差異分別為2.7°和2.0°，傾角差異分別為2.0°和2.4°，將上述3個(gè)斷層段合并為巍山-豐南斷裂段斷層(f1)； 第3段(圖7a)和第6段(圖7a)的走向差異為9.6°，傾角差異為2.0°，將其合并為灤縣斷裂北段(f3)。我們根據(jù)8個(gè)斷裂段附近的斷層和城市確定其名字，具體見(jiàn)圖7e。

圖 6 唐山地區(qū)4層密度聚類(lèi)的K值與有效斷層數(shù)目的關(guān)系圖(a—d)Fig. 6 Relationship between K value of four-level density clustering and the number of effective faults in Tangshan area(a—d).

圖 7 唐山地區(qū)第1、 2、 3、 4次聚類(lèi)結(jié)果圖(a—d)和最終斷層識(shí)別結(jié)果圖(e)Fig. 7 Graph of the first，second，third and fourth clustering results(a—d) and final fault identification results(e)in Tangshan area.

圖 8 巍山-豐南斷裂段(f1)的小地震擬合結(jié)果Fig. 8 The fitting result by using small earthquakes of Weishan-Fengnan Fault(f1).a 小地震的水平面投影； b 小地震的斷層面投影； c 小地震垂直于斷層面的橫斷面投影； d 小地震與斷層面的距離

表 1 利用本文方法求得的唐山地區(qū)各段斷層面的走向、 傾角和標(biāo)準(zhǔn)差Table1 Fault plane parameters of Tangshan region obtained by using the automatic fault identification method

3.3 斷層識(shí)別結(jié)果分析

為了分析本研究結(jié)果的可靠性，我們將研究結(jié)果與唐山斷裂帶的實(shí)際活動(dòng)斷裂(楊雅瓊等，2018)、 基于震源機(jī)制解的斷裂分段結(jié)果(楊雅瓊等，2016)和經(jīng)驗(yàn)性分段結(jié)果(萬(wàn)永革等，2008； 胡曉輝等，2019)進(jìn)行對(duì)比分析。萬(wàn)永革等(2008)經(jīng)驗(yàn)地劃分出陡河斷裂段、 巍山-豐南斷裂南段、 巍山-豐南斷裂北段、 灤縣-樂(lè)亭斷裂段和盧龍斷裂段； 楊雅瓊等(2016)采用震源機(jī)制解劃分出巍山-豐南斷裂南段、 巍山-豐南斷裂北段、 灤縣-樂(lè)亭斷裂段和盧龍斷裂段。本文方法的識(shí)別結(jié)果同樣包括陡河斷裂段、 巍山-豐南斷裂段、 灤縣-樂(lè)亭斷裂段和盧龍斷裂段，但將巍山-豐南斷裂南段與巍山-豐南斷裂北段合為一段。研究者在唐山斷裂帶分段方面存在一定爭(zhēng)議： 萬(wàn)永革等(2008)將其經(jīng)驗(yàn)性地分為3段； 楊雅瓊等(2016)基于震源機(jī)制解將其劃分為2段，分界點(diǎn)位于唐山市附近。 本研究將唐山斷裂帶劃分為2段，且分界點(diǎn)位于唐山市西南方向(圖7e)?？梢?jiàn)，本文結(jié)果與前人的研究結(jié)果存在異同，本研究結(jié)果與唐山斷裂的實(shí)際活動(dòng)更加一致(圖 10)。

本文研究結(jié)果中的陡河斷裂段(f8)、 盧龍斷裂段(f5)和灤縣-樂(lè)亭斷裂段(f2)與萬(wàn)永革等(2008)、 胡曉輝等(2019)所給出的結(jié)果差異見(jiàn)圖 9。總體上看，本文結(jié)果與萬(wàn)永革等(2008)結(jié)果的走向差異大于胡曉輝等(2019)的結(jié)果，反映了地震定位誤差對(duì)斷層面參數(shù)擬合的影響。胡曉輝等(2019)搜集了18組由不同學(xué)者或機(jī)構(gòu)給出的同一地震斷層面的參數(shù)差異，獲得的走向和傾角差異范圍分別為2°～34°和2°～28°。以此差異作為參考，除陡河斷裂段(f8)與胡曉輝等(2019)的傾角差異外，本文得到的走向和傾角與前人研究結(jié)果的差異值均在統(tǒng)計(jì)范圍內(nèi)(圖 9)。由于胡曉輝等(2019)基于中國(guó)地震臺(tái)網(wǎng)統(tǒng)一地震目錄擬合參數(shù)，其數(shù)據(jù)誤差較大(黃文輝等，2016)，且本文陡河斷裂段的傾角與萬(wàn)永革等(2008)的差異仍在統(tǒng)計(jì)范圍內(nèi)，故可認(rèn)為本文研究所得結(jié)果是可靠的。灤縣-樂(lè)亭斷裂段(f2)的斷層面參數(shù)相比杜晨曉等(2010)根據(jù)三維有限差分?jǐn)鄬铀矐B(tài)破裂動(dòng)力學(xué)模型反演得到的走向偏大12.2°、 傾角偏大9.3°。本研究的結(jié)果表明，唐山地區(qū)除陡河斷裂段(f8)和雷莊斷裂段(f7)外均為近直立斷層(表1)，這與前人的研究結(jié)果相同(李欽祖等，1980； 張之立等，1980； 萬(wàn)永革等，2008； 劉亢等，2015； 胡曉輝等，2019)。

圖 9 本文結(jié)果與前人研究結(jié)果走向(a)和傾角(b)的差異值Fig. 9 Difference of fault strike(a)and dip(b)between the results of this paper and those of previous studies.線段為胡曉輝等(2019)總結(jié)的斷層參數(shù)差異范圍

本研究所得結(jié)果除與上述研究一致外，在其余4條斷層中，有2條斷層的參數(shù)與地質(zhì)資料給出的斷層參數(shù)相符，1條斷層參數(shù)與其上發(fā)生的震源機(jī)制解參數(shù)相近。劉亢等(2015)和楊雅瓊等(2016)的研究結(jié)果表明，在古冶小地震密集區(qū)的東側(cè)存在一條新生地震帶，這與徐家樓-王喜莊斷層(f4)相吻合。劉亢等(2015)得到其為近NNE向的近直立斷層，與本研究的結(jié)果基本一致(表1)。本文新識(shí)別的灤縣斷裂北段(f3)與唐山地質(zhì)圖中NE向的盧龍斷層(圖 10)較為一致； 本文新識(shí)別的雷莊斷裂段參數(shù)與其上發(fā)生的震源機(jī)制解參數(shù)一致； 本文新識(shí)別的陳官屯斷裂段(f6)仍有待進(jìn)一步證實(shí)。

圖 10 唐山地區(qū)的主要斷層分布圖(劉亢等，2015)Fig. 10 Distribution of major faults in Tangshan area(after LIU Kang et al.，2015).F1陡河斷層； F2巍山-豐南斷層； F3古冶-南湖斷層； F4豐臺(tái)-野雞坨斷層； F5薊運(yùn)河斷層； F6寧河-昌黎斷層； F7灤縣-樂(lè)亭斷層； F8盧龍斷層； F9建昌營(yíng)斷層； F10滄東斷層； F11徐家樓-王喜莊斷層

斷層上地震的震源機(jī)制解，其節(jié)面之一的參數(shù)和斷層參數(shù)相似，因此常被用于驗(yàn)證所識(shí)別斷層參數(shù)的可靠性(Ouillonetal.，2008； Wangetal.，2013； 盛書(shū)中，2015； 盛書(shū)中等，2015； Skoumaletal.，2019)。為了分析本研究所得結(jié)果的可靠性，我們將其與搜集到的斷層上的震源機(jī)制解(表2)進(jìn)行對(duì)比分析。本文得到的盧龍斷裂段的走向小于震源機(jī)制解的走向(40.3°)，從圖7e 可見(jiàn)，由于該斷層段聚類(lèi)的震中分布呈分段線性，而本研究中將其視為一段擬合，由此導(dǎo)致走向差異過(guò)大。本文所得雷莊斷裂段的走向與2次雷莊斷裂段上的震源機(jī)制解走向差異分別為2.5°和6.5°，傾角差異分別為6.2°和7.2°，證實(shí)了本文新識(shí)別的雷莊斷裂段的可靠性。本文識(shí)別的巍山-豐南斷裂段于2020年7月12日發(fā)生唐山古冶區(qū)5.1級(jí)地震，本文所得斷裂段的走向相比3所機(jī)構(gòu)測(cè)定的震源機(jī)制解(表2)走向分別偏小10.6°、 12.6°和3.6°； 傾角分別偏大15.4°、 13.4°和15.4°。本文識(shí)別的灤縣-樂(lè)亭斷裂段于2021年4月16日發(fā)生唐山灤州市 4.3級(jí)地震，本文所得斷裂段的走向相比震源機(jī)制解(表2)走向偏大10.2°、 傾角偏大14.3°?？紤]到線性擬合可能會(huì)與實(shí)際斷層形態(tài)不符(曾憲偉等，2019)，因此，本文的研究結(jié)果與斷層上震源機(jī)制解存在較小差異是合理的。綜上所述，與現(xiàn)有震源機(jī)制資料相比(表2)，除盧龍斷裂段因呈分段線性導(dǎo)致較大誤差外，其他斷層參數(shù)與斷層上震源機(jī)制解的差異均在合理范圍內(nèi)。

表 2 震源機(jī)制解目錄Table2 Catalog of focal mechanism

總體來(lái)看，除盧龍斷裂段外，本文所識(shí)別斷層的傾角均大于震源機(jī)制解傾角。我們認(rèn)為造成傾角偏大的主要原因是雙差精定位后地震分布高度集中，帶狀分布更明顯，雙差定位后震中寬度分布的壓縮大于其在深度上的壓縮，且在斷層參數(shù)擬合的過(guò)程中僅采用內(nèi)部90%的小地震所在區(qū)域確定斷層面的4個(gè)頂點(diǎn)位置(萬(wàn)永革等，2008)，最終造成本研究擬合傾角均偏大的結(jié)果(圖8b)。

4 結(jié)論與討論

本研究改進(jìn)了DBSCAN算法，實(shí)現(xiàn)了由震源位置空間分布數(shù)據(jù)自動(dòng)識(shí)別斷裂段，并將該方法應(yīng)用于合成數(shù)據(jù)和唐山地區(qū)。合成數(shù)據(jù)和唐山地區(qū)斷層的成功識(shí)別表明： 1)本研究所提出的改進(jìn)的DBSCAN算法可以實(shí)現(xiàn)斷層段的自動(dòng)識(shí)別； 2)基于精定位數(shù)據(jù)，在唐山地區(qū)自動(dòng)識(shí)別出8條斷層段，本研究所得結(jié)果與經(jīng)驗(yàn)分段(萬(wàn)永革等，2008； 胡曉輝等，2019)和震源機(jī)制解分段(楊雅瓊等，2016)較為一致，表明本研究提出方法的合理性和有效性。

本研究所述方法在自動(dòng)識(shí)別斷層過(guò)程中，采用了從高密度層次到低密度層次的多密度層次聚類(lèi)，因此，本研究方法是從高密度向低密度的有序遞進(jìn)，與其他方法相比其優(yōu)點(diǎn)在于提高了對(duì)于低密度地震數(shù)據(jù)中斷層的識(shí)別能力。本研究方法相比前人的斷層識(shí)別方法(Scitovski，2017； Shangetal.，2018； Skoumaletal.，2019； Kameretal.，2020； Brunsviketal.，2021)部分地解決了經(jīng)驗(yàn)性選擇參數(shù)問(wèn)題，但仍存在以下缺陷： 1)有效事件數(shù)閾值的確定具有主觀性。在后續(xù)研究中，我們將考慮如何避免該主觀因素的影響。2)將實(shí)際斷層面近似為平面，難以反映實(shí)際斷層的三維復(fù)雜形態(tài)，如汶川地震斷層的鏟狀構(gòu)造(張培震等，2012)、 俯沖帶弧狀構(gòu)造(張澤明等，2021)等復(fù)雜情況。為了提高自動(dòng)識(shí)別斷層的精度，在后續(xù)的研究中，我們將考慮實(shí)際斷層的三維形態(tài)，對(duì)空間進(jìn)行網(wǎng)格化，獲取每個(gè)網(wǎng)格單元內(nèi)的斷層元，再以空間插值的方式連接所有斷層元，給出符合實(shí)際的三維斷層。3)本研究中忽略了深度因素，僅使用震中距離進(jìn)行聚類(lèi)。唐山地區(qū)斷裂以近直立的走滑型斷裂為主，且噪聲水平相對(duì)較低，本研究所得結(jié)果較為合理和可靠，在后續(xù)研究中，將使用震源距離進(jìn)行聚類(lèi)，以適應(yīng)復(fù)雜斷裂構(gòu)造的識(shí)別。

基于上述分析可見(jiàn)，本研究改進(jìn)的DBSCAN算法可以較好地實(shí)現(xiàn)斷層段的自動(dòng)識(shí)別，但仍存在以上亟需改進(jìn)的問(wèn)題。在后續(xù)研究中，我們將繼續(xù)對(duì)斷層段自動(dòng)識(shí)別方法進(jìn)行改進(jìn)，提高其自動(dòng)識(shí)別斷層能力，為活動(dòng)斷層的應(yīng)變分析、 地震地質(zhì)災(zāi)害評(píng)估和斷層滑動(dòng)趨勢(shì)等相關(guān)研究提供更加準(zhǔn)確的斷層資料。

致謝防災(zāi)科技學(xué)院萬(wàn)永革研究員為本文提供了模擬退火全局搜索-高斯牛頓局部搜索結(jié)合法的程序； 東華理工大學(xué)劉茜茜碩士、 防災(zāi)科技學(xué)院許鑫碩士、 李梟碩士和馮淦碩士對(duì)本文繪圖提供了幫助與指導(dǎo)； 審稿專(zhuān)家對(duì)本文提出了寶貴意見(jiàn)； 部分圖件采用GMT(Wesseletal.，1998)軟件繪制。在此一并表示感謝！