2017年3月渤海地震序列微震檢測與發震構造分析

馬婷 鄧莉 王曉山 宋程 譚毅培

1)天津市地震局，中國地震局地震工程綜合模擬與韌性抗震重點實驗室，天津 300201 2)河北省地震局，石家莊 050021

0 引言

地震序列指發震時間和震中位置相對集中的地震叢集。在華北地區，有數字化地震記錄的大震相對缺乏，而中小地震波形數據較為豐富，其數量優勢能夠在地震學研究中發揮重要作用(Brodsky，2019)。對中小地震序列發震構造的分析，成為華北地區地震活動性研究和地震危險性分析的重要基礎。

渤海海域地處郯城-廬江斷裂帶(簡稱郯廬帶)與張家口-渤海地震帶(簡稱張渤帶)的交匯部位，新生代主要發育張性和剪切性的斷裂，斷裂主要分NNE-NE向、NW向和近EW向3組(侯貴廷，2014)，地震地質構造較為復雜。渤海海域歷史上強震活動頻繁，史料記載1548年以來發生7級以上強震4次(國家地震局震害防御司，1995)。而近十多年華北地區有數字化地震波形記錄以來，渤海海域地震活動較弱。2017年3月24日渤海海域發生ML4.4地震，為十年來本區域發生的最大地震。根據河北省測震臺網地震目錄，3月24—28日，渤海海域發生地震18次，除3月24日ML4.4地震外，3月27日還發生1次ML3.5地震，地震震中位置如圖1 所示，ML2.0以上地震的震源參數見表1。研究此次地震序列的發震構造，有利于深入了解現今渤海海域地震活動特點，為區域地震危險性分析提供基礎資料。

圖 1 渤海地震序列2次ML3.0以上地震震中及本文所用臺站分布 紅色和藍色的圓圈分別對應ML4.4和ML3.5地震；黑色三角形表示本文使用的臺站；斷層數據引自鄧起東等(2007)

表 1渤海地震序列中ML2.0以上地震的震源參數

地震精定位和震源機制解是分析發震構造的重要手段，能夠較為精細地刻畫發震構造幾何特征(Long et al，2015；Yoon et al，2019)和地震破裂的錯動方式(易桂喜等，2019；王小娜等，2019)。前人對渤海地區歷史強震、1969年7.4級地震和區域中等地震的震源參數進行了復核，加深了對渤海海域地震活動與區域地質構造關系的認識(環文林等，1989；趙燕來等，1993；謝卓娟等，2008)。然而渤海海域連續運行的固定地震監測臺站稀少，區域臺網距離本次地震序列震中最近的臺站震中距約為80km(圖1)，對分析渤海地震活動的研究造成影響。如何利用現有的數字化地震波形，獲得可靠的渤海海域地震精定位和震源機制計算結果，值得深入研究。

隨著信息技術的發展，微震檢測和波形互相關技術能夠深入挖掘數字地震波形中攜帶的信息，逐漸成為地震序列分析中不可或缺的一環。常用的微震識別方法主要有3類，第一類基于地震波形與背景噪聲的振幅差別，多利用長時窗與短時窗的振幅或頻率等特征的變化(Allen，1978；Earle et al，1994)以及結合赤池信息準則(AIC)進行計算(劉希強等，2009)，此類方法已應用于中國地震臺網日常地震自動速報業務中；第二類為模板匹配方法，其利用地震波形相似性(Peng et al，2009)或音頻指紋算法(Yoon et al，2015)識別微震，近年來此方法已廣泛應用于余震序列(Wu et al，2017；Li et al，2018)、前震序列(Kato et al，2012；Ruiz et al，2014；Warren-Smith et al，2018)和震群(Frank et al，2018；Xue et al，2018)的遺漏地震檢測工作中；第三類基于人工智能的方法，以大量地震記錄作為訓練集學習地震波形的特征，進而在連續波形上搜尋這種特征從而檢測微震信號。國內外學者提出了多種人工智能地震和震相檢測算法(Perol et al，2018；Ross et al，2019；于子葉等，2018；趙明等，2019a、2019b；蔣一然等，2019；李健等，2020)，在地震序列分析和區域地震活動性研究中得到了應用(Ross et al，2020)。

本文以2017年3月渤海地震序列為研究對象，首先檢測序列中遺漏的微震事件，其次基于地震精定位和震源機制計算結果，對序列發震構造進行分析。

1 微震檢測

由于渤海地震序列地震波形相似度較高，同時序列發生在幾天之內，使得波形互相關計算量不大。因而，本文采用模板匹配識別方法(Peng et al，2009)，檢測渤海地震序列目錄遺漏的地震事件，從而提高地震序列的完整性。

模板匹配方法以較大地震記錄波形作為模板，在連續波形上進行互相關掃描，互相關達到一定閾值即為檢測到地震事件。遺漏地震事件檢測的計算步驟為：①選取地震目錄中所有ML2.0以上地震事件及其震相報告，選擇距離震中最近的3個寬頻帶臺站CLI、BDH和LUX的三分量數字化地震波形，截取地震事件S波到時前2s至后2s的波形作為模板；②將模板波形和連續波形通過2～8Hz的帶通濾波，再由每秒100個采樣點降為20個采樣點；③將模板波形在連續波形上進行滑動窗互相關掃描，滑動步長為1個采樣點(降采樣后采樣點間隔為0.05s)，計算3個臺站波形互相關系數的平均數；④取平均互相關系數大于9倍絕對離差中位數(MAD)，以及至少有一個臺站單臺互相關系數大于0.65作為閾值，挑選大于閾值的點作為檢測到的地震事件，其中臺網目錄中未記錄的地震作為遺漏地震事件；⑤利用遺漏事件水平向波形S波段最大振幅與模板地震的振幅比計算二者的震級差，從而估計遺漏事件的震級。

使用表1 中的6個ML2.0以上地震事件的波形作為模板，對3月24—28日時間段內的連續波形進行互相關掃描。將3月27日23時37分46.7秒的ML2.6地震作為模板，其在3月27—28日進行互相關掃描的結果，如圖2 所示。波形互相關值為3個臺的平均互相關系數，9倍MAD作為閾值(圖2 中紅色虛線)。圖2(a)結果顯示，在27日23時至28日4時檢測到較多的地震事件，顯示出在此時間段內多次發生地震，且其波形相似度較高。

圖 2 模板匹配微震檢測計算結果(a)掃描得到的3個臺的平均波形互相關系數，紅色虛線表示通過計算9倍MAD得到的地震檢測閾值；(b)互相關系數分布統計，縱坐標與圖(a)相同，統計區間步長為0.01，橫坐標為每個統計區間內的互相關系數個數，取對數坐標

圖 3給出了一個檢測到的目錄遺漏地震波形對比實例，所使用模板為3月28日3時19分55秒的ML2.3地震，截取模板S波到時前2s至后2s波形，經過2～8Hz帶通濾波(圖3 中紅色波形)。檢測到目錄遺漏地震為3月27日23時37分22秒的ML1.8地震。由于其后約25s發生一次ML2.6地震，振幅遠大于ML1.8地震，且震中距較遠臺站(如DOY、SUZ、DL2)的2個地震的記錄波形有所交疊，造成人工識別ML1.8地震較為困難。

圖 3 檢測到的目錄遺漏地震與模板地震波形對比紅色波形為模板地震；黑色波形為檢測到遺漏地震時間段的波形；紅色波形對比出的位置為檢測到遺漏地震波形的位置；由于需要清晰顯示ML1.8地震波形，ML2.6地震波形有所限幅

經檢測后較完整序列目錄的M-t圖如圖 4(a)所示，共檢測到目錄遺漏的地震事件32個，約為臺網目錄給出地震數量的1.8倍，其中檢測到ML2.0以上的遺漏地震2個。檢測到的遺漏地震較為集中地發生在3月24日ML4.4地震后15個小時內，以及3月27日23時至28日凌晨4時的時間段內(圖4(b))。

圖 4 經模板匹配檢測后較完整的序列目錄M-t圖(a)3月24—28日的M-t圖，其中紅色表示第一組地震序列，藍色表示第二組地震序列，黑色表示第三組地震，實線表示臺網目錄中給出的地震，虛線表示模板匹配方法檢測到的遺漏地震；(b)為圖(a)中3月27日23時至28日4時的放大圖，其中藍色實線表示臺網目錄中給出的地震，藍色虛線表示模板匹配方法檢測到的遺漏地震，均屬于第二組地震序列

通過模板匹配計算，我們觀測到此次渤海地震序列的地震可以分為2組，一組為主余序列，一組為震群。圖4(a)中紅色所示的第一組地震序列為ML4.4地震的余震序列，共檢測到14次地震，其地震頻度和強度均隨時間而衰減；藍色所示的第二組地震序列，集中發生在3月27日23時至28日4時之間，共檢測到35次地震，其中最大地震震級為ML3.5，次大地震震級為ML2.6，符合震群的定義；另外，臺網目錄中給出的3月27日9時7分22.5秒的ML1.6地震(圖4(a)中黑色實線所示)無法被其他任何模板地震檢測出，因而不屬于上述2組序列，判斷為一次孤立地震。此外，震群活動結束后的3月28日9時38分3.3秒仍有一個屬于第一組的地震發生，因而這2組地震并非先后發生，而是震群發生在余震序列活動期間。

第一組與第二組地震每組內地震的波形互相關系數較高，而2組地震間的波形互相關系數較低，其中一組的地震無法被屬于另一組的模板地震檢測出。選取記錄信噪比最好的BDH臺垂直向波形，經過2～8Hz帶通濾波后計算互相關系數，比較兩組地震序列內ML1.0以上地震間波形互相關系數占比。如圖5 所示，第一組序列除了模板ML4.4地震外，共有10個ML1.0以上地震，可計算得到10個互相關系數，根據互相關系數cc值分檔，其中滿足0.4≤cc<0.5、0.5≤cc<0.6和0.6≤cc<0.7的各1個，占比均為10%，滿足0.7≤cc<0.8的5個，占比50%，滿足0.8≤cc<0.9的2個，占比20%，10個互相關系數平均值為0.7033；第二組地震序列互相關系數平均值為0.8580。震群地震間的波形互相關系數平均值(0.858)高于余震序列地震間的平均系數(0.7033)，顯示出震群的波形一致性相對較高。

圖 5 第一組(a)和第二組(b)地震序列地震間波形互相關系數占比

2 震源機制計算與地震精定位

首先使用CAP方法(Zhu et al，1996)計算序列中2個震級最大地震的震源機制，選擇渤海周邊河北、山東和遼寧地震臺網的寬頻帶地震計記錄波形，體波濾波頻段為0.05～0.2Hz，面波濾波頻段為0.05～0.1Hz。參考基于深反射剖面對渤海灣地殼結果的研究(張成科等，2002)設定速度結構(表2)，vP/vS波速比設為1.73。

表2地震精定位和震源機制計算所用速度結構模型

ML4.4地震的震源機制計算結果如圖6 所示，節面I走向51°、傾角37°、滑動角-109°，節面Ⅱ走向254°、傾角55°、滑動角-76°；最佳擬合深度約為19km，擬合矩震級為MW4.08；結果顯示ML4.4地震震源機制以拉張為主。ML3.5地震的震源機制計算結果如圖7 所示，節面I走向41°、傾角90°、滑動角-18°，節面Ⅱ走向131°、傾角72°、滑動角180°；最佳擬合深度約為13km，擬合矩震級為MW3.49；結果顯示ML3.5地震震源機制以左旋走滑為主。震源機制計算結果給出2條正交的節面，判斷哪個節面為實際發震斷裂，需要結合地震序列精定位結果。

圖 6 ML4.4地震震源機制計算結果(a)波形擬合；(b)震源深度擬合

圖 7 ML3.5地震震源機制計算結果(a)波形擬合；(b)震源深度擬合

為了提高地震精定位結果可信度和精度，基于波形互相關對地震序列震相相對到時進行校正。本次地震序列發生在渤海中部(圖1)，100km內僅有2個短周期井下臺，且其記錄波形質量較差。最近的寬頻帶基巖臺(CLI臺)距震中約110km，該震中距范圍極易出現Pg、PmP、Pb震相交疊的現象，導致準確識別震相較為困難，因而臺網給出的震相報告存在人工識別震相誤差。我們采用波形互相關震相檢測技術精確拾取地震間同臺震相到時差，對所有目錄給出的地震和檢測到的遺漏地震的震相重新標定了到時。

震相校正的計算步驟為：①人工校核地震序列2組地震中最大的ML4.4和ML3.5地震的P波和S波震相到時，作為2組地震計算波形互相關的模板地震；②對于序列中的其他小地震，以微震檢測計算中得到的發震時刻為基準，根據模板地震的走時計算小地震震相到時的初始值；③將模板地震和小地震記錄波形經過2～8Hz帶通濾波，截取模板地震各臺站P波和S波到時前0.5s至后1.5s波形，以及小地震各臺站P波和S波到時初始值前0.5s至后1.5s波形，做波形互相關計算；④若波形互相關函數出現明顯的尖峰，則認為檢測到相應震相，作為校正后的小地震震相到時，若未出現明顯尖峰，則認為未檢測到震相；⑤將所有基于波形互相關檢測到的震相到時數據作為校正后的新震相報告。

基于校正后的新震相報告，利用雙差定位法(Waldhauser et al，2000)對序列進行精定位，使用的地殼速度結構見表2。ML4.4余震序列地震震中初始位置統一設置為臺網給出的ML4.4地震震中位置，初始深度設為CAP方法擬合的ML4.4地震深度19km；ML3.5震群地震震中初始位置統一設置為臺網給出的ML3.5地震震中位置，初始深度設為CAP方法擬合的ML3.5地震深度13km；反演計算選用奇異值分解(SVD)方法。

最終得到的精定位結果如圖 8 所示，ML4.4余震序列有5個地震可給出精定位結果(圖8(a)中紅色圓圈所示)，ML3.5震群有25個地震可給出精定位結果(圖8(a)中藍色圓圈所示)。其他地震因可用波形互相關檢測到的震相數量較少，無法得到高精度的定位結果，故在計算中被舍去。雙差定位法給出了地震序列的反演計算誤差，其中ML3.5震群精定位結果震中平均誤差為16.1m，最大誤差為30.1m，深度平均誤差為26.3m，最大誤差為37.9m；ML4.4余震序列精定位結果震中平均誤差為63.3m，最大誤差為328.1m，深度平均誤差為101.1m，最大誤差為386.9m。由于反演使用奇異值分解方法，雙差定位法給出的結果誤差是有意義的。

圖 8 渤海地震序列精定位結果(a)震中分布；(b)ML3.5震群在AA′剖面上的投影；(c)ML3.5震群在BB′剖面上的投影

結果顯示，2組地震發震構造走向均為NE向，與震源機制解中節面I走向基本一致。圖8(b)給出了ML3.5震群在與斷層走向正交的AA′剖面上的投影，顯示其發震構造為近直立的斷層，與震源機制解節面I傾角90°的結果一致。ML4.4余震序列精定位的地震僅有5個，無法刻畫斷層幾何特征，故未給出在剖面上的投影圖。

3 討論

2017年渤海地震序列發生在張渤帶與郯廬帶的交匯部位。張渤帶是一條晚第四紀、具有相當規模和較強活動性的NW向活動構造帶(徐錫偉等，2002)。郯廬帶是我國東部規模最大的一條巨型走滑斷裂帶，該斷裂走向為NNE-NE向。地震精定位結果顯示渤海序列走向為NE向，其發震構造與郯廬帶走向接近。郯廬帶渤海灣段歷史上發生過多次強震，其走向在山東段和渤海灣南段為NNE向(曹筠等，2018；顧勤平等，2020)，渤海灣北段為NE向(陳瑛等，2020)。地震序列震中所在的渤海中部正處于郯廬帶走向由NNE向NE轉變的過渡地帶。震源機制解給出的2次最大地震NE向節面走向分別為51°和41°，其更接近于渤海灣北段的走向。

渤海海域內郯廬斷裂帶東西兩側的構造存在差異，斷層西側多為高速區，中強地震震源深度淺于斷層東側(汪晟等，2017)。謝卓娟等(2008)通過統計分析認為渤海海域內地震深度多在10～20km的地殼中，渤中拗陷內地震震源深度一般較深。渤海地震序列震中位置處于郯廬帶西側(圖1)，CAP方法給出的2個最大地震深度為19km和13km，精定位后序列地震的深度均不超過20km，與謝卓娟等(2008)給出地震深度分布基本一致。1969年渤海海域M7.4地震的深度尚存在爭議，束沛鎰等(1983)通過遠震波形反演得到震源深度為25km，而周蕙蘭等(1989)同樣使用遠場波形得到初始破裂點深度為8km，不在謝卓娟等(2008)得到的中小地震震源深度分布優勢區間內。隨著渤海海域地震監測能力提升，渤海小震、微震震源深度的測定精度也將不斷提升。

此次渤海地震序列震中位于渤中盆地內的渤中凹陷。渤中盆地為NNE向“郯廬深斷裂”與NWW向“北京-塘沽-蓬萊深斷裂”的交匯區域，盆地構造演化受這2條重要的區域構造帶影響。渤中凹陷是盆地新生代沉降中心，充填的新生代地層達10km以上(陸克政等，1997)。渤中凹陷區斷層NE、NW和EW三組方向，主要分為2類：一類以延展斷層為特征的控盆構造，另一類以花狀構造和產狀較陡為特征的與走滑作用有關的次級斷層(侯貴廷，2014)。渤西地區拗陷的主伸展斷層比較低緩，傾角多為30°(童亨茂，2003)。傾角較低的伸展斷裂與ML4.4地震震源機制解節面I給出的參數(傾角37°、滑動角-109°)較為吻合；傾角較陡走滑型次級斷層與ML3.5地震震源機制解節面I給出的參數(傾角90°、滑動角-18°)較為吻合。

綜合地震精定位和震源機制計算結果，我們推測ML4.4余震序列發震構造為渤中凹陷內NE向低傾角的伸展性正斷層，ML3.5震群發震構造為NE向傾角較陡的次級走滑斷層。根據陸克政等(1997)和周斌等(2000)給出的渤中盆地構造剖面圖，得到渤海地震序列發震構造示意圖，如圖9 所示。ML4.4余震序列發生在伸展性正斷層上，淺部傾角較大，深部傾角逐漸變小。ML3.5震群發生在次級走滑斷層上，傾角較陡。由于中小地震破裂尺度較小，很難將其發震構造與已知的斷層相互對應。雖然可以對發震構造的幾何特征做出一定判斷，但確定地震發生在哪條斷層上仍十分困難。

圖 9 渤海地震序列發震構造示意圖(a)渤中盆地構造剖面圖，N+Q：新第三系和第四系，Ed：東營組，Es1-2：沙一、二段，Es3：沙三段，Es4：沙四段，據陸克政等(1997)修改；(b)視角為E向，斷層走向為NE向；斷層1表示ML4.4余震序列發震斷層，紅色圓圈為精定位后的余震序列震源位置；斷層2表示ML3.5震群發震斷層，藍色圓圈為精定位后的震群震源位置

各類微震檢測方法均有其優勢與不足，在實際應用中應根據研究對象和研究目標選擇使用的方法。長短窗類方法在信噪比較高時能夠可靠準確地識別震相(劉希強等，2009)，計算速度較快，但會漏掉低信噪比的信號(趙明等，2019b)。模板匹配類方法適用于檢測具有相近震源位置和震源機制的地震，對于波形高度相似的地震序列較為有效(Schaff，2010)，但其檢測精度取決于使用的模板，計算量隨模板數量增大呈指數級增加(Skoumal et al，2014)，且無法檢測到不與任何模板地震波形高度相似的地震。人工智能類算法優勢在于訓練好的模型具有泛化能力，可以發現訓練集不包含的新特征(趙明等，2019b)，在識別準確率和效率等綜合性能方面有所提升。

利用地震波形的相似性，除了能實現檢測微震的功能，還可以對序列中的地震進行聚類分析(王偉濤等，2012)。本文對渤海地震序列的分析，通過模板匹配算法檢測出序列可分為2組，并對2組地震分別進行精定位，發現2組地震的發震構造確實存在較大差別。若未經過波形互相關聚類，將2組地震合在一起分析，則可能無法得到可靠的發震構造分析結果。不可否認，基于人工智能的微震檢測技術在檢測能力和計算效率上顯示出了越來越大的優勢(趙明等，2019b；Perol et al，2018)，但在地震序列研究中，由于地震波性相似度較高，地震活動時間相對較短，模板匹配方法對波形相似的地震檢測能力較強的優勢體現較為明顯，同時能避免檢測出不屬于地震序列的其他地震的干擾。同時，與GPU加速技術(Beaucé et al，2017)、音頻指紋技術(Yoon et al，2015)等新技術相結合，在一定程度上彌補了模板匹配類方法計算效率較低的不足。通過給模板地震各臺站走時加入一個微小的變化值，M&L方法(Zhang et al，2015)或弱模板識別方法(吳夢羽等，2018)可以檢測到更多與模板地震“不太相似”的微震。因而，模板匹配方法仍然廣泛應用于地震序列微震檢測工作中(Ross et al，2019；Yao et al，2020；Sianipar，2020)。

對地震序列精定位結果誤差的估計對發震構造分析可靠性有重要作用。自助法(Bootstrap)是估計誤差的有力工具。對震相到時數據加入一定的隨機誤差(王清東等，2015)，或對震相到時、地震初始位置和速度結構同時加入誤差(Hardebeck，2013；姜金鐘等，2016)，通過重復定位幾百次確定精定位結果的誤差。Hardebeck(2013)分析了基于波形互相關校正后的相對到時誤差，發現其絕大多數在0.005s內。本文采用每秒100個采樣點的數據，0.005s為采樣間隔的一半。因而，可以認為經過波形互相關震相檢測后，每秒100個采樣點數據的相對到時誤差是非常微小的，對每一組地震相對定位結果的影響很小。本文2組地震的初始位置分別置于臺網給出的ML4.4和ML3.5震中位置，震源深度基于CAP擬合結果。鑒于人工識別的震相到時存在一定誤差，導致震中定位結果存在一定不確定性，ML4.4和ML3.5兩個模板地震初始位置的不確定性會造成本文精定位結果中2組地震的相對位置存在誤差。同時，速度結構的不確定性，尤其是震源區P波和S波速度的誤差，對雙差定位也有影響。利用Bootstrap方法定量估計模板地震初始位置和速度結果造成的相對定位結果誤差，是今后研究和改進的方向。

4 結論

本文對2017年3月渤海地震序列進行了微震檢測、震源機制計算和地震精定位計算，得到如下認識：

(1)通過模板匹配方法，檢測到目錄遺漏地震32個，約為臺網目錄給出地震數量的1.8倍，其中檢測到ML2.0以上遺漏地震2個；

(2)基于波形互相關特征，發現渤海地震序列分為2組，一組為ML4.4地震主余序列，一組為最大震級ML3.5的震群，另有一個ML1.6地震與其他地震波形相似度較低，可能為一個孤立的地震事件；

(3)地震精定位結果顯示2組地震均為NE走向，震源機制解計算得到ML4.4地震發生在低傾角正斷層，ML3.5地震發生在高傾角走滑斷層；

(4)結合研究區地質構造相關研究結果，推測ML4.4余震序列發震構造為渤中凹陷內NE向低傾角的伸展性正斷層，ML3.5震群發震構造為NE向傾角較陡的次級走滑斷層。

致謝：中國地震局地球物理研究所韓立波研究員對本文撰寫給予指導和建議，天津地震臺提供了數據資料支持，審稿專家提出了寶貴的修改意見，部分圖件采用GMT軟件包(Wessel et al，1995)繪制，在此一并表示感謝。