文璽翔, 沈旭章*, 周啟明

1 中山大學地球科學與工程學院, 廣東省地球動力作用與地質災害重點實驗室, 廣州 510275 2 南方海洋科學與工程廣東省實驗室(珠海), 廣東珠海 519082

0 引言

2019年10月12日22時55分,廣西壯族自治區(qū)北流市與廣東省化州市交界發(fā)生MS5.2地震,震中位于22.16°N、110.52°E,震源深度10 km,震中烈度VI度. 此次地震是廣西地區(qū)自2016年蒼梧地震以來又一次5級以上的地震,致使震源區(qū)5個鄉(xiāng)鎮(zhèn)受災,經濟損失超過500萬元.

震中地區(qū)自新生代以來處在喜馬拉雅構造域與濱太平洋構造域的復合部位(如圖1所示),新構造運動主要包括間歇性整體抬升為主的運動以及沿NE向斷裂的斷塊差異運動(李冰溯等, 2019).李細光等(2007)和張培震等(2013)結合目前的震源機制解、震區(qū)發(fā)震構造以及GPS速度場等研究成果,得到該區(qū)域構造應力場的主壓應力方向為NW-SE向.研究區(qū)域主要分布北東向蕉林斷裂、北西向石窩斷裂及近南北向的新豐斷裂. 根據野外地質資料(潘建雄和黃日恒,1995),新豐斷裂為逆斷性質,傾向南東; 石窩斷裂為左旋正斷層,傾向南西,處于北西向巴馬—博白斷裂帶的分支斷裂.自2016年蒼梧地震后,該地區(qū)進入另一個地震活躍期(任鎮(zhèn)寰和羅振暖,1998).深入了解該區(qū)域發(fā)震構造對于華南地區(qū)地震活動、抗震設防以及地震預測預報工作都有重要意義.

圖1 地震及臺站分布圖(a) 臺站分布及研究區(qū)主要斷層. 主圖中三角形表示短周期密集臺陣位置,圓圈表示余震分布,顏色深淺表示地震發(fā)生時間; 圖中顯示了主震和前震的震源機制解 (數據來自于中國地震局地震預測研究所和廣東省地震局),并且展示了該區(qū)域三條主要斷裂: 石窩斷裂,新豐斷裂,蕉林斷裂; (b) 震源區(qū)域放大圖,其中圓圈的大小對應于地震震級(震級大小分布于-0.4～2.2之間),顏色深淺代表地震發(fā)生的時間.Fig.1 Map of earthquakes and stations(a) Map of stations and the major faults in the studying area; Triangles represent the locations of the high dense short-period seismic array, circles represent the distribution of aftershocks, and different colors represent the time of the earthquake; The focal mechanism solutions are also presented (These data came from the Institute of Earthquake Forecasting (IEF) of the China Earthquake Administration (CEA) and the Guangdong Earthquake Agency); The figure shows three major faults in this region, the NW Shiwo fault, the SE Xinfeng fault, and the NE Jiaolin fault; (b) Zoomed window for the region enclosed by the box in (a), where the size of the circle indicates the magnitude of the earthquake (the magnitudes of these events ranged from -0.4 to 2.2) and the different colors represent the time of the earthquake.

目前不同研究機構結果均顯示(廣西壯族自治區(qū)地震局,蘇珊; 中國地震局地震預測研究所; 中國地震局地球物理研究所; 中國地震臺網中心,趙博等; 防災科技學院,Seismology小組): 此次北流地震震源機制解為走滑性質,節(jié)面I走向為NWW,傾角約為70°,節(jié)面II走向NNE,傾角近乎垂直,主壓應力P軸方向為NW向,同該地區(qū)的構造應力場方向一致(許忠淮等,1989; 張培震等,2013).盡管根據目前已有的震源機制解結果與野外地質勘探調查、地震烈度分布及余震空間展布可以基本確定NWW向節(jié)面與真實的發(fā)震斷層石窩斷裂相對應(李冰溯等,2019; 王小娜等,2020; 閻春恒等,2019).但是受觀測資料所限,目前對北流地震的認識還面臨以下問題: (1) 震源機制解得出該節(jié)面運動性質為右旋走滑,與野外地質調查得到的石窩斷裂為左旋走滑的運動性質不同; (2) 根據定位和震源機制解結果,5.2級主震前2 s發(fā)生了4.2級前震,震中位于22.16N°,110.53E°,震源深度9.5 km,約在主震北西1 km處.這兩次地震的震源機制解雖均為走滑型,節(jié)面走向卻不同,因此可能是由于不同斷裂作用而導致的結果(王小娜等,2020; 周斌,2019),但是缺乏足夠觀測資料證明; (3) 此外,后續(xù)得到的部分余震震源機制解結果與主震前震都不同(郭培蘭等,2019).

大量余震事件空間分布和震源機制解是深入認識和了解發(fā)震構造的最重要信息.近年來,機器學習方法在識別微震事件中取得了較大突破,此類方法已經被應用于多個區(qū)域并取得了較有意義的豐碩成果.如Ross等(2018)使用卷積神經網絡基于南加利福尼亞區(qū)域的1820000條地震波形所訓練的模型,通過交叉驗證,在P波到時拾取上與人工拾取平均誤差在0.023 s,并且在測定初動極性方面,與人工測定精度相比有95%的準確率; 趙明等(2019)利用U型卷積神經網絡在處理首都圈地震臺網以及四川地震臺網的數據時,Pg、Sg拾取到時上表現(xiàn)出了相較于傳統(tǒng)震相自動識別方法更優(yōu)的均方根誤差,準確率和查全率; Wang等(2019)提出基于深度卷積神經網絡(DCNNs)VGG-16結構的PickNet算法,在處理日本地區(qū)300個地震事件并與日本氣象廳(JMA)提供的震相目錄進行對比后,表現(xiàn)出了速度快、精度高以及P、S震相拾取數量更多的特點; Zhou等(2019)提出的基于卷積神經網絡的DetNet與循環(huán)神經網絡的PpkNet算法,在處理汶川余震的數據上表現(xiàn)出了地震信號檢測的完備性以及震相拾取的穩(wěn)定性與精確性,且在S波到時拾取上較好的表現(xiàn);Mousavi等(2020)提出的EQTransfomer方法,在處理日本2016年鳥取地震的波形數據后,得到了相比人工識別震級分布更廣,小震數量更多的地震目錄,且P波和S波到時拾取精度相對于人工拾取誤差在0.05 s左右.Zhou等(2021b)、Jiang等(2021)、廖詩榮等(2021)和蘇金波等(2021)將機器學習方法應用于2021年云南漾濞6.4地震研究,獲得了高精度、高完備性的余震序列目錄.Wong等(2021)使用機器學習方法和波形互相關方法對四川盆地地震進行研究,構建了一個高精度地震目錄, 進一步完善了該目錄的震級的完整性.Zhou等(2021a)將基于深度學習的PhaseNet方法應用于四川威遠頁巖氣區(qū)域臨時臺網記錄的連續(xù)地震數據中,獲得了相比于中國地震臺網約60倍的地震事件,證明了機器學習震相拾取方法在密集臺網中的適用性.

本研究將基于2019年10月12日廣西北流5.2級地震震源區(qū)震后約30天的短周期密集臺陣觀測資料,使用EQTransformer模型(Mousavi et al.,2020),對余震進行識別及拾取,進而使用Hypoinverse和Hypodd方法對結果進行精定位,挑取信噪比較高的地震事件進行震源機制解反演.根據余震空間分布及震源機制解特征,對該區(qū)域中強地震發(fā)震構造進行探討.

圖2 地震目錄中4次地震事件波形發(fā)震時刻標注在每幅圖的上方,臺站號位于每幅圖左側,按震中距排列.Fig.2 Waveforms of 4 earthquake events in regional earthquake catalogThe origin time is marked at the top of each subgraph, and the station numbers are located at the left, arranged according to epicentral distance.

圖3 959臺站24小時連續(xù)波形檢測結果黃色星表示識別出來的地震事件,紅色星表示多個臺站共同識別到的地震事件.Fig.3 24 hours continuous records and detection results of earthquakesYellow stars indicate detected events, and red stars indicate events recorded by multi-stations.

圖4 不同臺站同一地震事件的波形(a)、(b) 為臺站959(SP424)和臺站983(SP327)記錄到的同一個地震事件; (c)、(d) 為臺站1025(SP209)和臺站1011(SP309)記錄到的同一個地震事件.Fig.4 Waveforms recorded by different stations of the same earthquake(a) and (b) indicate the same seismic events recorded by stations 959 (SP424) and 983 (SP327); (c) and (d) indicate the same seismic event recorded by stations 1025 (SP209) and 1011 (SP309).

圖5 單臺站余震識別結果三角形表示臺站位置,不同顏色表示識別出的余震事件數目.(a) 每個臺站單獨識別出的微震數目; (b) 精定位之后每個臺站的識別結果; (c) 精定位與地震臺網共有的11個余震事件的臺站識別結果.Fig.5 The detection results of single stationTriangles indicate station locations, and different colors indicate the number of aftershocks. (a) the number of events detected individually from each station; (b) the number of events detected by each station after relocation; (c) single station about the 11 events in the regional catalog during the same period.

1 數據

北流地震之后一周,中山大學地球科學與工程學院在震源區(qū)布設了由120個EPS-M6Q短周期地震儀組成的密集地震臺陣,臺間距3～5 km,自2019年10月18日起進行了為期約一個月的觀測,在觀測期間記錄到了大量震級較小的余震.在密集臺陣觀測同期,區(qū)域地震目錄共記錄到了余震13次,這些地震事件在本研究布設的密集地震臺陣上都有質量較好的地震事件波形記錄.如圖2為已有地震目錄中4次震級較小的地震事件波形記錄,在波形記錄中可以看到清晰的P波和、S波.由于本研究中臺陣觀測密度遠遠高于區(qū)域地震臺網,因此也記錄到了非常多的區(qū)域地震臺網沒有記錄到的余震事件.使用機器學習方法,可以識別出更多震級較小的余震.

2 機器學習方法及微震識別

本文使用最新深度學習神經網絡算法(Mousavi et al.,2020)對北流的短周期地震數據進行處理.該方法包含由多任務結構組成的一個主編碼器和三個獨立的解碼器,解碼器和編碼器主要由一維卷積層(1D convolutions)、雙向和單向長短期記憶門(long-short-term memories,LSTM)、NIN網絡層(Network-in-Network)、殘差連接層(residual connections)、前饋層(feed-forward layers)、轉換層(transformer)和自注意力層(self-attentive)組成.編碼器在時域中分解地震信號,并通過在它們的時間依賴性上提取高維特征.然后,解碼器利用這些信息將高級特征映射到與地震信號相關的三個概率序列: 即每個時間點是否存在地震信號、P波和S波.

相比其他方法,該方法主要優(yōu)點為: (1) 效率更高且降低了人為主觀性的影響; (2) 所受信噪比與檢測區(qū)間的影響更小; (3) 對于內存和計算的需求更低; (4) 具有更好的泛化性.

圖3展示了利用連續(xù)記錄波形確定余震事件的過程.首先對輸入數據進行基于三分量數據標準差的振幅歸一化,并將地震事件、P波和S波的識別閾值分別設置為0.4、0.5和0.5,這樣做是在保證每個臺站識別及拾取結果相對可靠的情況下,盡可能增加關聯(lián)地震事件的數目.本研究所使用的震相關聯(lián)算法基于每個事件(包括其P波和S波到時)的識別時間,將不同臺站之間識別誤差不超過15 s的事件歸為同一事件(Mousavi et al.,2020).在SP424臺站(儀器編號959)連續(xù)24小時的觀測記錄中,使用機器學習方法識別出了8個疑似余震事件.進而根據其他臺站的記錄結果對這些疑似事件進行了核實和確定,如SP327臺站(儀器編號983)同樣在20時34分識別出一例事件,當某一地震事件由至少三個(如3個P或者兩個P加1個S,等等)臺站記錄到時才會將其關聯(lián)到初始地震目錄中.圖4為短周期密集臺陣不同臺站確定的余震事件波形.

使用該方法,本研究最終識別出了地震數據中記錄到的大量余震事件波形,并對結果進行關聯(lián),得到有441個地震事件的地震目錄(如圖5所示),大約是同時間段的廣東和廣西聯(lián)合測定得到地震事件數目的34倍.

3 地震精定位及震源機制解反演

機器學習方法識別確定出的地震事件形成了關聯(lián)后的震相文件.為了對余震事件空間分布進行深入分析,我們使用Hypoinverse方法(Klein,2002)對地震事件先進行初定位.初定位使用的地殼模型如圖6所示.

圖6 初定位使用的地殼速度模型橙色線表示S波波速,藍色線表示P波波速.Fig.6 Crustal velocity model used for initial locatingThe orange line shows the S-wave velocity, and the blue line shows the P-wave velocity.

為消除所選速度模型對非震源區(qū)附近的影響,確認關聯(lián)后識別地震的可靠性,本研究進一步使用雙差定位法(Waldhauser,2001)對余震序列進行重定位.其原理是將同一臺站記錄到的兩個相鄰地震(事件對)的觀測走時差與理論走時差的殘差(即“雙差”) 表示為一個方程,利用最小二乘法求解所有臺站和事件對組成的方程組,來獲得地震事件的相對位置.該方法近年來在國內外地震序列精定位研究方面得到了較好的應用,如房立華等(2018),Yang等(2022).經過上述的精定位過程,最終得到了一個包含299個事件的地震目錄.

本研究中使用臺站比較密集,對很多余震事件觀測的方位覆蓋較好,為研究震源機制解提供了可靠保證.因此我們使用gCAP(general Cut And Paste)全矩張量反演方法求解了余震震源機制解及震源矩心深度.由于反演過程對波形記錄質量和方位覆蓋要求較高,為得到較為可靠穩(wěn)定的震源機制解,本研究選取發(fā)生在主震附近且信噪比較高的65個余震波形進行震源機制解反演.該方法通過對P波震相及其后續(xù)震相(Pnl)和S波(或面波)賦予不同的權重,依據理論與實際波形之間的擬合誤差函數,使用網格搜索法來反演地震矩張量的六個獨立分量,即走向、傾角、滑動角、標量地震矩M0, 以及ISO(各向同性)參數ζ和CLVD(補償線性偶極子)參數(Zhu and Ben-Zion,2013).

求解處理過程主要包括以下步驟: (1) 在SAC文件中寫入頭段信息,截取重定位后的地震事件,選用起震時刻前20 s至后100 s的數據; (2) 對截取波形得到的地震數據進行去均值、去尖以及去除儀器響應等操作,接著對數據進行重采樣,并進行分量旋轉; (3) 采用頻率-波數法(FK)來計算(Zhu and Rivera,2002) 擬合所需要的格林函數模型,使用全球地殼模型CRUST1.0,網格中心位于110.5°E,21.5°N.反演的震源深度范圍2～16 km,每0.5 km計算一個模型; (4) 將非雙力偶分量ISO和CLVD約束為0,擬合求解地震的最佳雙力偶節(jié)面解.反演中Pnl和S波截取波形窗長分別為25 s和40 s,二者的濾波范圍為0.05～0.1 Hz,走向,傾角和滑動角的搜索步長為5°,每隔0.5 km反演一次震源機制獲取最佳矩心深度.

4 結果與討論

4.1 余震精定位結果

初定位的結果(如圖7所示)顯示余震事件主要發(fā)生在研究區(qū)域的三個位置,即主震及前震附近(測線AA′-BB′)、蕉林斷裂附近(測線CC′-DD′),以及石窩斷裂南端(測線EE′-FF′).主震附近的余震平面有近NS和EW向兩個延展方向,深度多在3～11 km之間,震源位置隨深度有向東移動的趨勢,且大部分余震發(fā)生在較前震和主震更淺的位置,震中位置也更接近于前震; 蕉林斷裂附近與石窩斷裂南側的余震,震源深度集中在5 km左右,雖然這些余震在已有的固定臺站地震目錄中并沒有記錄,且余震數目較少,對于指示孕震斷層空間展布的作用較弱,但是卻可以證明受北流地震的影響,該地區(qū)在至少60 km的范圍內有多條次級斷層處于活動狀態(tài).

重定位的結果如圖8所示,重定位后的地震目錄共包括299個地震事件,而同期的廣東地震臺網的余震目錄中只有13個地震事件.主震附近余震源深度主要在4～11 km之間,與已有地震目錄中的事件深度相近.但與之不同的是,該區(qū)域還有許多發(fā)生在淺部5 km以內的地震,且余震中分布僅僅表現(xiàn)出沿EW向的條帶特征.余震震中位置多集中在前震附近,位于主震西北1～3 km范圍內; 其它區(qū)域的余震,震源深度集中在3～7 km左右,水平和垂向上沒有明顯的延展趨勢.大部分余震發(fā)生在布設臺站后的15天以內,但是震中位置與震源深度同時間沒有明顯相關性.

由主震附近精定位結果的三維空間分布(如圖9所示)可以看到微震序列隨震源深度有沿NEE向移動的趨勢,在8～10 km的范圍內,微震分布集中,有指示斷層面參數的作用(萬永革等,2008),該斷層走向約為NWW-SEE,傾向約在SSW-SW之間,傾角約為70°,與目前多個機構得到的主震震源機制解節(jié)面產狀吻合(廣西壯族自治區(qū)地震局,蘇珊;中國地震局地震預測研究所;中國地震局地球物理研究所; 中國地震臺網中心,趙博等;防災科技學院,Seismology小組),因此確定該斷層對應主震的孕震斷層,即石窩斷裂.此外,還可以由微震分布推測該地區(qū)存在一條分布于地下4～12 km的斷層,走向約在E-NEE之間,傾角近90°,其產狀與前震震源機制解節(jié)面相吻合(王小娜等,2020),這條斷裂與石窩斷裂相切,且二者走向較為接近,因此很有可能是該斷層先發(fā)生錯動,而后觸發(fā)石窩斷裂活動.

圖7 初定位結果(a) 余震序列分布圖,藍色圓點表示識別的余震事件,紅色圓點表示區(qū)域地震目錄中的事件; (b)—(g) 測線AA′-FF′ 的余震序列縱剖面,綠色五角星表示主震和前震的位置.Fig.7 Initial locations of the aftershocks(a) Distribution of the aftershocks. The blue circle represents the picked events, the red circle represents the events in the catalog; (b)—(g) The longitudinal profile of aftershocks, and the green star represents the position of the mainshock and foreshock.

圖8 重定位結果(a) 余震序列分布圖,顏色表示識別的余震事件的時間,藍色圓點表示區(qū)域地震目錄中的事件; (b)—(g) 測線AA′-FF′的余震序列縱剖面,綠色五角星表示主震和前震的位置.Fig.8 The results of relocation(a) Distribution of the aftershocks, color shows the elapsed time of the recognized aftershocks relative to the main shock, the blue circle shows the events in the regional earthquake catalog; (b)—(g) The longitudinal profile of aftershocks, and the green star shows the position of the main earthquake and foreshock.

圖9 主震附近重定位結果(a) 余震序列分布三維圖, (b) 余震分布三視圖; 不同顏色的球表示不同的發(fā)震時刻,紫色平面表示石窩斷裂的空間展布,青色平面表示推測斷裂的空間展布.Fig.9 3D view of relocation results near main shock(a) 3D map of aftershock distribution; (b) depth distribution of aftershocks in latitude and longitude profiles; Different colors of balls show different seismogenic moments. The purple plane shows the spatial distribution of the Shiwo fault, while the cyan plane shows the spatial distribution of the presumed fault.

圖10 震源機制解結果(a) 震源機制解結果,不同顏色的震源球表示不同的震級范圍; (b)—(g) 測線AA′-FF′的余震震源機制解縱剖面.Fig.10 The results of focal mechanism solution(a) The result of focal mechanism solution. Different colors of the focal spheres show different magnitude; (b)—(g) The longitudinal profile of aftershock focal mechanism solution from survey line AA′ to FF′.

圖11 主震附近震源機制解結果(a) 震源機制解結果,不同顏色的震源球表示不同的震級范圍; (b)—(c) 測線AA′-BB′的余震震源機制解縱剖面.Fig.11 The results of focal mechanism solution near the mainshock(a) The result of focal mechanism solution. Different colors of the focal spheres show different magnitude; (b)—(c) The longitudinal profile of aftershock focal mechanism solution from survey line AA′ to BB′.

4.2 震源機制解反演結果

挑選的65個余震事件分布及其震源機制解如圖10所示,主震附近的余震(具體震源機制解結果如圖11所示) 矩震級基本在2.0以內,多為走滑型,震源機制解的兩個節(jié)面走向分別在NNW-NWW之間以及NNE-NEE之間,傾角接近直立,P軸方向基本與該地區(qū)主壓應力方向保持一致.震級較小的地震震源深度在5～8 km,震級較大的地震震源深度集中在9～11 km.該部分余震震源機制解與前震的震源機制解類似,節(jié)面Ⅰ走向為NNW-SSE,節(jié)面Ⅱ走向為NEE-SWW,震源深度在9～10 km,位于主震NW約1～3 km,進一步證明了上述推測斷層的存在.

蕉林斷裂附近的余震,震源機制解多表現(xiàn)為走滑型以及逆沖走滑型,兩組節(jié)面傾角在40°～60°,節(jié)面走向多為NW-NWW,震級多在MW1.3～1.9之間; 石窩斷裂南側地震,震源機制解多為過渡型和走滑型,節(jié)面Ⅰ傾角在70°～90°,節(jié)面Ⅱ傾角在40°～60°,兩組節(jié)面走向近SN、EW向,震級較小,都在1.3以下.

由于前震和主震的孕震過程并不相同(王小娜等,2020; 周斌,2019),暗示主震的產生可能是由于上述另一條斷層的錯動而觸發(fā),且在主震后,兩條斷層同時處于活動狀態(tài),隨之又產生了不同震源機制解的余震.此外,蕉林斷裂北東端以及石窩斷裂北端的孕震斷裂同樣受到此次北流地震的影響而處于活動狀態(tài),產生了諸多震級較小的地震.由于該部分地震數量較少,分布較為隨機,難以確定孕震斷層的具體位置及產狀.

5 結論

本文基于廣西北流5.2級地震震后布設在震源區(qū)120個短周期地震儀組成的密集臺陣資料,使用機器學習方法識別了441個余震事件,并通過Hypoinverse和Hypodd方法進行精定位.從精定位的299個地震事件中選擇65個波形清晰、有多臺站記錄的事件使用gCAP方法進行震源機制解反演.根據主震、前震以及余震精定位及震源機制解結果,結合震區(qū)地質構造情況和烈度情況,對廣西北流5.2級地震進行發(fā)震構造分析,得到如下主要結論:

(1)震源區(qū)位于石窩斷裂與新豐斷裂的交匯處,可能存在相互切割的斷裂.主震附近的余震主要發(fā)生在其北西約1～3 km的位置,震源深度在4～11 km以內,震源機制解多表現(xiàn)為走滑型,節(jié)面走向與前震更接近.結合精定位與震源機制解結果,確定主震區(qū)域主要存在兩條斷裂,即走向NWW-SEE向的石窩斷裂,傾角約70°,以及走向NEE-SWW,傾角近90°的一條斷裂,該斷裂活動性更強,石窩斷裂可能是伴隨該斷裂的影響而發(fā)生錯動.

(2)蕉林斷裂附近同樣識別出了較多余震,這些余震震源深度在3～8 km以內,震級相對于主震區(qū)域更小,震源機制解多表現(xiàn)為逆沖型與走滑型.表明蕉林斷裂在受北流地震的影響后處于活動狀態(tài).此外石窩斷裂以南也識別出了部分余震,震級約在1.0左右,多為走滑型與過渡型,該區(qū)域并沒有出露于地表的斷裂,但是處在石窩斷裂與蕉林斷裂的交匯處,同樣可能存在次級隱伏斷裂.

致謝本文在處理波形數據的過程中使用到了SAC以及Obspy(Beyreuther et al.,2010) 軟件,余震識別使用了S.Mostafa Mousavi等提供的模型和EQTransformer軟件,地震定位使用了Fred W.Klein提供的HYPOINVERSE和Waldhauser提供的HypoDD程序,作震源機制解使用了美國圣路易斯大學朱露培教授分享的gCAP程序,本文圖件繪制主要用到了 GMT(Wessel et al.,2019),在此一并表示感謝.