利用模板匹配技術識別大崗山水庫庫區爆破事件

杜瑤 阮祥 戴仕貴 邵玉平 胥珍 余洋洋 羅勇

（中國成都 610041 四川省地震局）

0 引言

地震事件的識別與拾取在地震數據處理過程中非常重要，根據識別手段的不同，可劃分為人工經驗分析和計算機自動識別。隨著地震臺網的日益密集及專用地震臺網的不斷布設，對水庫庫區、頁巖氣開采區、石油鉆探區等工程場地的地震監測逐步開展起來，監測數據越來越多，人工拾取地震事件的效率和精度已遠不能滿足地震數據處理的需求。同時，隨著微小震級地震事件研究的不斷深入，對如何高效率、高精度、低錯誤率地識別微小地震事件提出了更高的要求。在專業地震臺網監測的水庫庫區及附近，地震類型較復雜，構造地震及人工地震活動等均可能會在監測區內出現。由于人工爆破波形與天然地震波形具有極高的相似性（尹欣欣等，2020），因此，如何有效地將不同類型的地震加以區分，從波形上分辨天然地震、人工爆破，提高地震定位、編目的質量，對于蓄水前后水庫庫區及附近地震活動研究，以及分析預報工作中正確判定人工爆破等異常都具有較重要的實際意義。

國內外針對地震類別的識別已開展了一些相關研究。張春賀等（2006）在小爆破識別與研究的基礎上，通過周期—頻度譜、波形特征等對爆破事件的自動識別進行了研究，自動識別概率達到了77%以上；崔鑫等（2016）通過對天然地震和人工爆破的頻譜分析認為，天然地震和人工爆破的記錄具有明顯不同的時頻特征；尹欣欣等（2020）利用波形互相關的模板匹配方法對爆破及天然地震事件進行震相識別，對于人工爆破，準確率為80%，對于天然地震，準確率為100%。在日常編目工作中發現，人工識別時震級較小的近震及爆破極易混淆，因此，建立一種能夠及時準確識別天然地震和人工爆破的方法，把人工爆破與震級較小的地方震區別開來，以得到真實的區域地震活動變化很有必要（鄭秀芬等，2006；包淑嫻，2011）。

大崗山水電站是在四川境內大渡河干流近期開發的大型水電工程之一。電站正常蓄水位 1 130 m，最大壩高210 m，總庫容7.42 億m3，電站裝機容量2 600 MW。電站于2005 年開工建設，2015 年建成。大崗山水庫于2015 年5 月28 日導流底孔下閘蓄水，這標志著水庫正式開始蓄水發電，該水庫屬典型的高山峽谷型高壩大庫容水庫。

大崗山水庫庫區位于青藏高原東南緣，地屬川西高原區。庫區所在區域地處以鮮水河—安寧河—小江斷裂為界的川滇斷塊與大涼山斷塊結合部位SN 向構造帶北段，該處為SN 向與NW、NE 向等多組構造的交匯復合部位，區域構造上位于由金坪斷裂和磨西斷裂所圍限的黃草山斷塊的西側邊緣，庫區內的主要構造形跡以近SN、NNW、NW 向的褶皺和斷裂為主。主要斷裂有大渡河斷裂、磨西斷裂、大發斷裂、安寧河斷裂等及規模較小的斷層，與庫水直接接觸的巖性主要為花崗巖，占庫區的大部分庫段，在水庫田灣支流有泥盆系、二疊系的一套碎屑巖出露，在庫中段有三疊系碎屑巖出露。從工程地震條件來看，庫區屬于斷裂塊狀巖體亞型，可能誘震強度為中等或強烈地震（《中國水力發電工程》編審委員會，2000）。區域內有較強構造地震背景，其中，1975 年1 月15 日大壩西邊的康定6.2 級地震震中距壩址約31 km，是歷史上距庫區最近且影響最大的地震。與低烈度地區水庫誘發地震不同，高烈度地區水庫在蓄水前就存在明顯地震活動（程萬正，2013），因此，在水庫蓄水后如何判斷是否庫水加卸載誘發的地震活動則需要開展更細致的工作。2014 年1 月、2015 年4 月、2015 年10 月等大崗山水庫蓄水前后庫區及附近均發生過大量微小地震，并且地震類型較復雜，有天然地震、工業爆破及蓄水后誘發地震等。

1 方法及數據

2012 年12 月大崗山水庫地震監測臺網軟硬件設備安裝完成，2013 年2 月完成調試運行，同年3 月開始考核運行，同時開始監測及觀測數據的處理。庫區共布設了8 個高增益短周期速度記錄地震儀，數據采集器為EDAS-24GN3 型，地震計為RSFS-1A 型，儀器觀測頻帶為2 s—40 Hz，布局孔徑NS20 km×EW40 km，臺站沿庫區均勻展布，監測范圍包括可能誘發地震的重點監視區段，監測震級下限為ML0.5，震中定位誤差≤1 km（阮祥等，2017）。根據大崗山水庫地震監測臺網分布，并結合其外圍地區的區域地質構造環境和歷史地震情況，將大崗山水庫庫區地震活動監測研究區域定為29.2°—29.9°N、101.9°—102.5°E。

采用Zhang 等（2015）對微震同時進行檢測及定位的方法，即匹配定位方法（Matching and location，簡稱M&L 方法），并在GeoTaos 軟件平臺下進行模板掃描與定位。該方法需要精確的模板地震事件，同時更注重微弱小震事件的拾取與定位。M&L 方法計算流程（圖2）及原理如下 （Zhang et al，2015）。

圖2 M&L 方法模板掃描流程Fig.2 Flow diagram of M&L method

圖1 研究區域位置及地震臺站分布陰影部分為研究區Fig.1 Map showing study area and stations

根據震群序列重新定位后將震級較大的地震作為模板事件，以模板事件位置為中心在經度、緯度、深度等3 個方向上進行網格化，計算每個可能的微震位置（格點）與模板事件位置之間的參考震相（S 波）在同一地震臺站上的走時差，然后按此走時差對所有臺站分量記錄的模板事件參考震相（S 震相）與相應臺站分量的連續數據進行滑動互相關，并計算互相關波形疊加后的平均相關系數CC和信噪比SNR。當平均相關系數和信噪比大于閾值時認為檢測到了1 個微震，并把其位置確定在最大相關系數的格點位置（曾憲偉等，2017）。其中，由模板事件與檢測出的微震事件之間的位置差異所導致的走時差計算公式如下（Wen，2006；Wen et al，2010；Zhang，2013）

其中，dDk為模板事件與可能的微震事件之間相對差異所導致的臺站k震中距的差異；dh為2個事件間深度的相對變化；分別為地震相位P 的走時分別對模板震中距D（水平慢度）和模板深度h（垂直慢度）的導數，利用一維速度模型，分別計算每個臺站及其相關地震相位的。走時差Δt(k,P)的計算對參考模型依賴較小。

研究中獲取了部分爆破作業的準確時間數據，并根據這些資料在地震臺網目錄獲取事件，從連續波形中挑選了信噪比高、震相清晰的5 次爆破事件（圖3、表1）。利用這5 次已知的爆破事件作為模板事件，進行去均值、去線性趨勢等預處理以及1—10 Hz的帶通濾波，連續波形重采樣率至50 Hz，因爆破事件樣本及記錄到爆破的臺站數量較少，故直接使用觀測報告中定位的震相文件，沒有再進行重新定位。由于震中距越大的臺站記錄波形特征越不明顯，近臺震相波形特征明顯，因此模板選取近臺記錄的爆破事件（圖4）。

圖3 研究區域爆破事件分布Fig.3 Distribution of epicenters of blasting events in the study area

圖4 爆破模板事件波形（a）HCP 臺2014-03-10 T 18:44:56.1；（b）CAK 臺2014-03-12 T 18:18:36.4；（c）BLG 臺2014-03-24 T 18:01:43.5；（d）CAK 臺2014-04-09 T 18:10:39.5；（e）CAK 臺2014-04-28 T 18:07:50.6Fig.4 Template waveforms of blasting events

表1 爆破事件統計Table 1 Statistics of blasting events

選取大崗山水庫蓄水前后2014 年3 月1 日至2015 年12 月31 日庫區地震活動監測研究區域（29.2°—29.9°N，101.9°—102.5°E）8 個臺站（XIX、AJW、ZMC、NTW、HCP、BLG、CAK、XMC）三分量地震儀記錄的連續波形，利用M&L 方法檢測定位。

由于近震爆破與天然地震在波形上有差異，前者P 波能量強于S 波，更容易從噪聲中檢測出來，因此選取P 波及后5 s 的時間作為匹配相位和匹配窗口，水平、垂直搜索距離分別為2 km，水平搜索網格數為20，垂直搜索網格數為10，模板相位信噪比閾值為10，綜合考慮匹配的爆破數量及質量，確定3 個及以上臺站同時記錄到時檢測相關系數閾值為0.3，只有2 個臺站記錄到時為0.8，僅1 個臺站記錄到時為0.95，所采用的四川地區一維平均速度模型見表2。圖5 分別選取了與模板事件2、4、5 同一臺站（CAK 臺）記錄的事件波形匹配出的事件，2 個不同波形之間的相關系數分別為0.628 5、0.773 4、0.409 7。從圖5 還可見，掃描出的事件波形與模板事件均有震中距為數十千米的小噸位爆破的特征：①縱波在垂直分向的初動向上；②垂直分向上呈振幅較大的脈沖型；③有較明顯的反射和折射波等。同一臺站記錄的波形相似度極高，可確定為同類型事件。

圖5 M&L 方法檢測出的事件Fig.5 Comparison of the template seismograms (blakc traces) with signals detected in the continuous waveforms (red traces) by using M&L Method

表2 四川地區一維平均速度模型Table 2 One-dimensional average velocity model for Sichuan area

2 結果及分析

利用M&L 方法，剔除因斷記等因素造成的相關系數＞1 的誤匹配結果，再經過人工復合，最終獲得23 次疑爆記錄，其中，5 次疑爆事件與已知爆破事件波形位置吻合且相關系數為1.0，可以確定為爆破。圖6 為爆破模板事件及掃描出的爆破（疑爆）事件分布。在圖6 中劃分了A、B、C、D、E 等5 個事件集中區域，其中，A 區域有疑爆事件12 次（含爆破1 次）；B 區域有疑爆事件2 次（含爆破1 次）；C 區域有爆破事件1 次；D 區域有疑爆事件4 次（含爆破1 次）；E 區域有疑爆事件4 次（含爆破1 次）。匹配結果為已知爆破模板事件個數的4.6 倍，有效地提高了爆破事件檢出率，同時將匹配識別出的爆破事件形成目錄與原分析編目報告中的事件進行比對檢驗，區別并剔除原報告中人工誤判的爆破事件，提高報告準確性。將匹配出來的疑爆事件及距離最近的臺站（CAK）記錄到的NS 向波形分別在Google Earth 上投影（圖7），由于疑爆事件發生時間均為2014—2015 年，而Google Earth 的影像則是2020 年8 月，所以存在部分爆破環境與當時實際環境不一致的可能，并且可能受速度模型的影響，事件定位也會有一點偏離，但偏離量很小，基本在0.5 km 之內。由圖7 可見，掃描出的爆破（疑爆）事件附近存在采石場或蓄水前修建大壩的場地，這與爆破實際情況較符合，通過模板搜索出的爆破（疑爆）事件詳見表3。

圖6 爆破模板及掃描出的爆破（疑爆）事件分布Fig.6 Distribution of templates of blasting events and detected blasting or suspicious blasting events

圖7 匹配出來的疑爆事件及CAK 臺記錄到的NS 向波形在Google Earth 上的投影Fig.7 Projection of detected blasting or suspicious blasting events and seismograms recorded by CAK station in NS direction at satellite map

表3 通過模板掃描出的疑爆及爆破事件Table 3 Catalogue of suspicious blasting and blasting events detected by template

對表3 中的發震時刻在2014 年3 月1 日至2015 年12 月31 日人工編目觀測報告中進行搜索，找到并剔除了觀測報告中人工誤分的爆破及疑爆事件8 次（圖8），表明可以利用M&L 模板匹配方法來有效區別不同類型的地震事件。

圖8 人工編目觀測報告中被誤判為爆破事件的地震震相報告Fig.8 Seismic phase report of misjudged blasting event in the original observation reports

3 結論

選取大崗山水庫庫區及附近波形記錄較好的5 次爆破作為模板事件，采用M&L 方法對2014 年3 月至2015 年12 月連續波形進行檢測與定位，并進行人工復核，最終獲得23 次爆破（或疑爆）事件，這些事件均呈現出震中距為數十千米范圍內小噸位爆破的特征：①縱波在垂直分向的初動向上；②垂直分向上呈振幅較大的脈沖型；③有較明顯的反射和折射波等。掃描出的爆破（或疑爆）事件與模板事件有較高的相關性，相關系數達0.773 4。將部分掃描出的事件投影至Google Earth，從影像圖上可見附近存在采石場或蓄水前修建大壩的場地，這與爆破實際情況較符合，重定位偏離基本在0.5 km 之內。將識別出來的爆破（或疑爆）目錄與人工編目觀測報告進行對比，剔除了8 次人工誤判的爆破事件，有效地提高了觀測報告的精度。大崗山水庫庫區及附近在水庫蓄水前存在大量的人工活動，而真實的地震活動背景才對庫區蓄水前后地震活動的研究具有參考意義。微小地震與爆破在人工分析處理中容易混淆，而利用人工地震活動作為模板事件在連續波形中掃描定位的方法，可有效地剔除人工地震活動事件。建議建立非天然地震數據庫，利用匹配定位等方法區分人工活動地震事件，在地震監測工作中復檢工程地震臺網監測報告，落實因人工活動形成的震群異常。此外，因爆破等人工活動地震事件震源位置精確，可將其用于反演速度結構等研究工作中。

日本產業技術綜合研究所雷興林研究員提供Geotaos 軟件，并給予了指導與幫助，在此表示感謝。