張?zhí)K祥 盛書中* 席 彪 房立華 呂 堅 王甘嬌 張 瀟

1)東華理工大學，地球物理與測控技術學院，南昌 330013 2)長江大學，油氣資源勘探技術教育部重點實驗室，武漢 430100 3)中國地震局地球物理研究所，北京 100081 4)江西省地震局，南昌 330039

0 引言

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

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

1 方法

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

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

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

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

(1)

其中，Pi為第i個對象的eps鄰域對象數(shù)量，n為地震目錄D中的對象總數(shù)。

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

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

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

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

上文簡述了本研究提出的基于改進的DBSCAN聚類算法自動識別斷層的方法，其具體實現(xiàn)步驟如下：

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

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

(2)

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

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

(3)

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

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

(4)

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

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

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

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

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

2 理論測試

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

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

圖 1 a 人工合成地震數(shù)據(jù)； b—d 分別添加總地震事件5%、 10%和20%隨機地震后的地震空間分布圖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.黑色圓圈表示隨機地震(即噪聲)，其余為斷層上的地震事件

2.2 自動識別斷層面方法的測試

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

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

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

圖 3 添加5%噪聲的合成數(shù)據(jù)第1次聚類(a)和第2次聚類(b)的K值與有效斷層數(shù)目關系圖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ù)自動識別斷層面的結果Fig. 4 Automatic fault plane recognition results for synthetic data.a 添加5%噪聲的合成數(shù)據(jù)未合并的識別結果，矩形框為最小事件數(shù)閾值選取區(qū)域； b 添加5%噪聲的合成數(shù)據(jù)斷層面合并后的識別結果； c 添加10%噪聲的合成數(shù)據(jù)斷層面合并后的識別結果； d 添加20%噪聲的合成數(shù)據(jù)斷層面合并后的識別結果； 用不同顏色標示識別出的斷層段

2.3 測試結果分析

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

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

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

3.1 數(shù)據(jù)

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

圖 5 唐山地震序列的震中(a)、 深度(b)和震級(c)分布圖Fig. 5 Distribution of epicenter(a)，depth(b)and magnitude(c) of the Tangshan earthquake sequence，China.矩形框為最小事件數(shù)閾值的選取區(qū)域

3.2 斷層自動識別結果

重復與理論測試部分相同的數(shù)據(jù)處理步驟，自動識別出唐山地區(qū)的斷層段?；趫D 5 框內(nèi)的97個地震事件數(shù)據(jù)，確定M＿effect為90。其次，對數(shù)據(jù)進行不同密度層次聚類，基于上述終止原則進行了4次聚類，其最優(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。最終自動識別出不同密度層次的37個聚類(圖7a—d)，滿足條件的有效斷層有11個。

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

圖 6 唐山地區(qū)4層密度聚類的K值與有效斷層數(shù)目的關系圖(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次聚類結果圖(a—d)和最終斷層識別結果圖(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)的小地震擬合結果Fig. 8 The fitting result by using small earthquakes of Weishan-Fengnan Fault(f1).a 小地震的水平面投影； b 小地震的斷層面投影； c 小地震垂直于斷層面的橫斷面投影； d 小地震與斷層面的距離

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

3.3 斷層識別結果分析

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

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

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

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

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

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

表 2 震源機制解目錄Table2 Catalog of focal mechanism

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

4 結論與討論

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

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

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

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