杜 航 楊 云 鄭江蓉 王 俊 張 揚(yáng) 宮 杰

（江蘇省地震局, 南京 210014）

引言

北京時間2021 年5 月22 日02：04，青海果洛州瑪多縣發(fā)生MS7.4 地震（34.59°N，98.34°E），震源深度17 km。震中位于瑪多黃河鄉(xiāng)，距縣城38 km，距西寧380 km，震中海拔高度4 258 m。地震表現(xiàn)出無明顯前震、余震較多的特征，主震后1 小時內(nèi)發(fā)生MS3.0 以上余震10 次，最大余震為MS4.3，一天內(nèi)發(fā)生MS3.0 以上余震36 次，最大余震為MS5.1，為主震-余震型地震。中國地震臺網(wǎng)中心對震源機(jī)制解進(jìn)行反演得出本次地震為走滑型事件，地震破裂持續(xù)45 s，震后30 s 內(nèi)釋放大量能量，最大滑移量4.5 m，地表破裂超過160 km（鄧文澤等，2021）。本次地震為汶川8.0 級地震以來，中國境內(nèi)發(fā)生的震級最高、地表破裂最長的地震事件。本文將對此次地震震前b值進(jìn)行研究，尋求該區(qū)域b值與應(yīng)力的變化關(guān)系，為該區(qū)域后續(xù)地震預(yù)報提供科學(xué)依據(jù)和參考。

Gutenberg 等（1944）發(fā)現(xiàn)震級M與地震頻次N之間存在統(tǒng)計學(xué)關(guān)系logN=a?bM，a表示研究區(qū)域地震活動水平，b表示研究區(qū)域不同震級的相對分布。Kun 等（2013）通過數(shù)值模擬研究驗(yàn)證了震級、頻次的數(shù)學(xué)關(guān)系。研究人員在進(jìn)行巖石實(shí)驗(yàn)時發(fā)現(xiàn)介質(zhì)的應(yīng)力水平不斷增加時，b值會呈現(xiàn)逐漸減小的趨勢（Mogi，1962；Scholz，1968；Rivière 等，2018）。Mousavi 等（2017）對整個板塊進(jìn)行研究，同樣發(fā)現(xiàn)b值的變化與介質(zhì)的應(yīng)力變化存在一定關(guān)系。因此，研究學(xué)者使用b值研究區(qū)域應(yīng)力的大小（Nuannin 等，2012；Scholz，2015）。易桂喜等（2006）對龍門山斷裂帶b值進(jìn)行分析，得到斷裂帶區(qū)域的強(qiáng)震危險性；劉艷輝等（2015）對青藏高原東南緣b值進(jìn)行時空掃描，發(fā)現(xiàn)在大地震發(fā)生前后b值會出現(xiàn)較大變化，地震發(fā)生前b值長時間處于低值；Tormann 等（2015）研究得到全球平均b值約為1，中強(qiáng)震多發(fā)地區(qū)b值小于1，小地震多發(fā)地區(qū)b 值大于1； Nanjo 等（2018）利用b值與應(yīng)力關(guān)系，得到海槽地震的初始破裂部位；謝卓娟等（2019）利用b值與應(yīng)力關(guān)系研究首都圈地震活動性；曾憲偉等（2021）通過b值分析得到寧夏吳忠?靈武地區(qū)的強(qiáng)震危險性。前人已開展大量區(qū)域應(yīng)力變化與b值之間的關(guān)系研究，在此基礎(chǔ)上Nievas 等（2020）提出可以利用b值變化研究區(qū)域應(yīng)力變化，尋求b值變化與地震發(fā)生的關(guān)系。本文對瑪多MS7.4 地震震前地震序列進(jìn)行分析，研究地震前b值的變化特征，為研究區(qū)域地震預(yù)測和地震發(fā)生風(fēng)險提供參考。

1 區(qū)域地質(zhì)背景

新生代以來，青藏高原不斷受到歐亞板塊和印度板塊的擠壓，使青藏高原呈現(xiàn)擴(kuò)張趨勢，表現(xiàn)出板塊運(yùn)動活躍、地震多發(fā)的特點(diǎn)，21 世紀(jì)以來發(fā)生7 級以上地震10 次（Molnar 等，2009）。巴顏喀拉塊體是青藏高原的次一級構(gòu)造單元，該塊體被東昆侖、甘孜-玉樹-鮮水河、龍門山3 條大型斷裂帶包圍，成為中國目前地震活動最為強(qiáng)烈的地區(qū)之一（張培震等，2003；Rui 等，2016；聞學(xué)澤，2018）。其中，巴顏喀拉塊體東緣最為活躍，十年內(nèi)相繼發(fā)生了2008 年汶川8.0 級地震、2013 年蘆山7.0 級地震、2017 年九寨溝7.0 級地震（圖1）。

圖1 青藏高原主要活動斷裂及歷史地震分布Fig. 1 Simplified distribution of main active faults and historic earthquakes in the Qinghai-Tibet plateau

東昆侖斷裂位于巴顏喀拉塊體北邊界，EW 走向，長約2 000 km，是一條巨型左旋走滑斷裂帶，傾角55°～85°。該斷裂自西向東走滑速率不斷降低，瑪沁段為12.5 mm/a，至瑪曲段降為5 mm/a，至塔藏段已不足3 mm/a（Kirby 等，2007；Ren 等，2013）。有記錄以來，該斷裂發(fā)生7 級以上地震6 次，最大震級為2001 年MS8.1 地震，地表破裂達(dá)426 km（徐錫偉等，2002；鄧起東等，2003；Xu 等，2006）。

甘孜-玉樹-鮮水河斷裂帶位于巴顏喀拉塊體南邊界，為左旋走滑斷裂帶。學(xué)者對該斷裂帶滑動速率有不同見解，聞學(xué)澤、石峰等認(rèn)為甘孜-玉樹段滑動速率為8～14 mm/a，但 Chevalier 等認(rèn)為甘孜-玉樹段滑動速率僅有5～8 mm/a（聞學(xué)澤等，2003；石峰等，2013；Chevalier 等，2018）。學(xué)者們普遍認(rèn)為鮮水河段滑動速率為10～13 mm/a（Zhang，2013；Chen 等，2016）。有記錄以來，該斷裂發(fā)生7 級以上地震17 次，最大震級為1997 年MS7.9 地震。

龍門山斷裂位于巴顏喀拉塊體東邊界，為逆沖右旋走滑斷裂帶，由多個近平行斷裂組成。21 世紀(jì)以來，該斷裂活動加劇，2008 年汶川MS8.0 地震、2013 年蘆山MS7.0 地震都發(fā)生在該斷裂上。

巴顏喀拉塊體除四周存在3 條大型斷裂帶外，其內(nèi)部也發(fā)育著一系列活動斷裂?，敹?甘德斷裂帶位于東昆侖斷裂帶以南70 km，為全新世活動的左旋走滑斷裂，整體呈NW-SE 走向，長約650 km，水平滑動速率小于6～8 mm/a，活動速率明顯小于位于其南北兩邊的甘孜-玉樹-鮮水河斷裂帶和東昆侖斷裂帶（熊仁偉等，2010）。

2 數(shù)據(jù)及方法

2.1 地震目錄數(shù)據(jù)

本次研究地震目錄由中國地震臺網(wǎng)提供，目錄選取時間為2009 年1 月1 日至2021 年4 月30 日。青?，敹郙S7.4 地震研究范圍選取震中所在的青藏高原東北部巴顏喀拉塊體內(nèi)部，南北范圍選取東昆侖斷裂帶、甘孜-玉樹-鮮水河斷裂帶之間區(qū)域（32°N～37°N），東西范圍覆蓋瑪多-甘德斷裂帶且包含本次地震所有余震（96°E～100°E）。

由于巴顏喀拉塊體內(nèi)部發(fā)育著眾多活動斷裂，導(dǎo)致區(qū)域內(nèi)5 級以上地震較為分散。其中瑪多MS7.4 地震及其余震沿瑪多-甘德斷裂帶分布，呈NW-SE 走向，空間分布長約150 km（詹艷等，2021；張喆等，2021；潘家偉等，2021）（圖2）。

圖2 研究區(qū)域中強(qiáng)震（MS≥5.0）及瑪多地震余震（ML≥3.0）分布圖Fig. 2 Simplified distribution of major earthquake（MS≥5.0）and Maduo earthquake aftershocks（ML≥3.0）in the study area

在b值計算時，MS震級適用于大范圍、長時間、高震級（一般指一級板塊以上范圍、100 年以上時間、6 級以上震級）范圍內(nèi)b值計算（陳培善等，2003）。本次研究所選區(qū)域較小，時間尺度為12 年，地震目錄99.8%震級小于5 級，因此本次研究采取ML震級進(jìn)行計算。

地震活動性研究中，余震是否應(yīng)該剔除尚未有統(tǒng)一答案。眾多學(xué)者認(rèn)為中強(qiáng)地震引發(fā)的余震會造成研究區(qū)域地震頻次變化，因此在b值變化研究時，應(yīng)對地震目錄內(nèi)的余震進(jìn)行剔除（Keilis-Borok 等，1980；Reasenberg，1985；陳凌等，1998；韓曉明等，2016）。也有學(xué)者認(rèn)為余震去除會影響地震目錄完整性，在研究時應(yīng)保持地震目錄的原始性，選取最小完整性震級以上地震進(jìn)行計算即可（史海霞等，2018）。為研究余震序列對b值計算結(jié)果的影響，本文采用完整地震目錄和去余震目錄 2 種地震目錄研究b值的時空變化，為準(zhǔn)確對比2 種地震目錄得到的b值差異，除地震目錄外其它參數(shù)均相同，其中余震剔除時采用G-C 方法，即對5 級以上地震余震時間窗915 天、余震空間窗半徑R0內(nèi)的余震進(jìn)行剔除（陳凌等，1998）：

2.2 地震b 值計算方法

目前b值計算方法主要有2 種，最小二乘法和最大似然法。最小二乘法是直接通過震級-頻次的關(guān)系確定b值的大小，該方法簡單易行，但其穩(wěn)定性受地震頻次影響，當(dāng)?shù)卣疠^少時結(jié)果容易出現(xiàn)偏差（Weichert，1980；Sandri 等，2007）。最大似然法作為概率方法具有更好的穩(wěn)定性，可降低少數(shù)地震對區(qū)域整體b值結(jié)果的影響。由于研究區(qū)域臺站間距較大，導(dǎo)致區(qū)域內(nèi)最小完備震級較大，地震頻次較少，因此，采用最大似然法計算本區(qū)域b值，該方法最早由Aki（1965）提出：

式中，M為震級；Mc為最小完備震級； ?M為震級分檔，本次地震目錄震級分檔為0.1。b標(biāo)準(zhǔn)誤差可由以下公式進(jìn)行估算（Aki，1965）：

式中，N為樣本量，即地震目錄中大于等于Mc的 所有地震的個數(shù)?？梢钥闯鯪與δb成反比，因此提高樣本數(shù)量可降低b值誤差。該方法先將研究區(qū)域網(wǎng)格化，然后計算每一個網(wǎng)格的b值。計算b值時需要震級為連續(xù)隨機(jī)變量，而實(shí)際震級計算時，會對震級進(jìn)行歸檔，使震級只保留一位小數(shù)。因此震級分檔?M越小，b值的計算結(jié)果越準(zhǔn)確，本次地震目錄?M為0.1。

2.3 最小完備震級

最小完備震級通過最大曲率法確定，該方法認(rèn)為震級-頻次曲線的一階導(dǎo)數(shù)最大值對應(yīng)的震級即最小完備震級（Wiemer 等，2000），研究區(qū)域最小完備震級為ML2.1（圖3）。最小完備震級反映該地區(qū)的地震監(jiān)測水平，青海地區(qū)雖地震頻發(fā)，但人口稀少、臺站密度較小，因此最小完備震級較大。為避免因地震監(jiān)測能力差異導(dǎo)致地震目錄不完整，影響b值的時空變化研究，本次研究將選取ML2.1 以上地震對瑪多7.4 級地震震源區(qū)進(jìn)行b值研究。

圖3 研究區(qū)域2009 年1 月1 日?2021 年4 月30 日震級-頻次圖Fig. 3 The frequency-magnitude distribution of earthquakes in the study area during January, 2009?April, 2021

3 地震b 值時空掃描特征

3.1 b 值時間掃描特征

為研究時間段區(qū)域內(nèi)b值隨時間的變化規(guī)律，選用ML2.1 以上地震目錄進(jìn)行計算，為避免地震數(shù)量較少對結(jié)果準(zhǔn)確性造成影響，以4 個月為時間段進(jìn)行b值時間掃描，分別計算余震剔除前后b值隨時間的變化（圖4）。由圖4 可知，2 種地震目錄得到的b值整體變化趨勢相同，在5 級地震發(fā)生前兩者吻合的較好，但在5 級以上地震發(fā)生后的一段時間內(nèi)存在明顯差異。圖4 中2 條水平直線分別表示2 種地震目錄b值均值，其中完整目錄b值均值為0.86、去余震目錄b值均值為0.85，均值相差0.01。

圖4 b 值時間曲線Fig. 4 The temporal variation of b value

從b值時間掃描特征來看，是否剔除余震對b值時間變化影響不明顯，僅在中強(qiáng)震發(fā)生后一段時間內(nèi)存在差異。這種現(xiàn)象是由于研究區(qū)域中強(qiáng)震多表現(xiàn)為主震-余震型地震序列類型，余震主要集中在主震發(fā)生后，因此是否剔除余震對中強(qiáng)震發(fā)生后一段時間較為明顯，對其他時間段影響較小。同時，雖然中強(qiáng)震余震持續(xù)時間較長，但隨時間的推移余震不斷減小，而該區(qū)域最小完備震級較大，不管是否剔除余震，震級小于ML2.1 的地震都被刪除，后續(xù)較小余震對2 種地震目錄未產(chǎn)生影響，中強(qiáng)震發(fā)生后一段時間后2 種地震目錄的b值變化曲線又互相重合。因此，是否剔除余震對中強(qiáng)震余震窗內(nèi)b值時間變化存在一定影響。

瑪多MS7.4 地震前1 年b值持續(xù)低于均值，并處于下降趨勢，b值有上升趨勢時地震發(fā)生，擴(kuò)大到區(qū)域內(nèi)其它5 級以上地震，也發(fā)現(xiàn)此規(guī)律。同時，2 種地震目錄的b值曲線均存在3 段較長的低b值時間段，且起止時間差異不大，分別為2009 年5 月?2011 年12 月、2014 年5 月?2016 年4 月、2019 年9 月至今，基本為中強(qiáng)震分割點(diǎn)，地震發(fā)生后b值明顯上升，短時間內(nèi)又下降到較低位置，并一直處于較低位置直至下次地震發(fā)生。這是由于地震發(fā)生前，該區(qū)域處于應(yīng)力積累階段，b值處于較低位置，地震發(fā)生后應(yīng)力得以釋放，區(qū)域內(nèi)b值恢復(fù)到正常水平。

3.2 地震b 值空間掃描特征

b值時間掃描結(jié)果表明，b值變化與中強(qiáng)震發(fā)生存在一定關(guān)系，可以作為地震預(yù)測時間上的參考。但研究區(qū)域較大，時間掃描特征僅能反映研究區(qū)域時間尺度上發(fā)生地震的可能，無法準(zhǔn)確給出地震發(fā)生的空間區(qū)域。要實(shí)現(xiàn)通過b值分析給出研究區(qū)內(nèi)中強(qiáng)震發(fā)生的空間預(yù)測范圍，就要對b值進(jìn)行空間掃描，分析研究區(qū)域內(nèi)b值處于異常狀態(tài)的地區(qū)。由b值時間掃描特征知，是否剔除余震對中強(qiáng)震余震窗內(nèi)b值存在一定影響。為驗(yàn)證余震剔除是否對b值空間掃描特征同樣具有影響，在b值空間掃描特征研究中同樣計算ML2.1 及以上震級去除余震地震目錄及完整地震目錄，分別計算2009 年1 月?2021 年4 月、2009 年1 月?2012 年12 月、2013 年1 月?2016 年12 月、2017 年1 月?2021 年4 月的b值空間特征，研究瑪多地震前該區(qū)域b值空間變化特征及b值空間分布與中強(qiáng)地震的關(guān)系?，敹郙S7.4 地震前去除余震目錄及完整目錄b值空間掃描結(jié)果分別如圖5、圖6 所示。

根據(jù)我國西部地震多發(fā)的特點(diǎn)，本文將研究區(qū)域經(jīng)緯度掃描步長設(shè)為0.1°，并以此進(jìn)行網(wǎng)格化（易桂喜等，2007）。由于研究區(qū)域地震較少，網(wǎng)格內(nèi)地震掃描半徑設(shè)為70 km，為保證結(jié)果的可信度，僅選取地震數(shù)量超過30 個的網(wǎng)格進(jìn)行計算，以每個網(wǎng)格點(diǎn)為單位計算b值，網(wǎng)格內(nèi)的b值誤差最大為0.18b，小于地震b值的變化范圍，可信度較高。再根據(jù)每個網(wǎng)格點(diǎn)的b值做出研究區(qū)域b值空間分布，若網(wǎng)格點(diǎn)地震數(shù)量達(dá)不到要求則進(jìn)行空白處理。

圖5（a）、圖6（a）分別為2009 年1 月1 日?2021 年4 月30 日去除余震目錄及完整目錄b值空間掃描結(jié)果，可以看出兩者在瑪多震中東南部低b值區(qū)域存在差別，圖5（a）中該區(qū)域b值約為0.7，圖6（a）中該區(qū)域b值約為0.8。區(qū)域內(nèi)整體b值范圍為0.5～2.0，瑪多MS7.4 震中處于低b值區(qū)域內(nèi)部，低b值區(qū)域以瑪多MS7.4 震中所在的瑪多-甘德斷裂為中心向兩側(cè)擴(kuò)展，該區(qū)域b值整體低于0.6，瑪多MS7.4 震中所在區(qū)域?yàn)閎值最低區(qū)域，b值約為0.5，表明該斷裂震前構(gòu)造應(yīng)力處于較高狀態(tài)。由圖5（a）、圖6（a）可知，b值的空間分布與中強(qiáng)地震的發(fā)生有較好的對應(yīng)，為研究瑪多MS7.4 地震前震中區(qū)域b值變化過程，本文將地震目錄時間進(jìn)行分割。

圖5（b）、圖6（b）分別為2009 年1 月1 日?2012 年12 月31 日去除余震目錄及完整目錄b值空間掃描結(jié)果，兩者在稱多西南部低b值區(qū)域存在差別，圖5（b）中該區(qū)域b值約為0.6，圖6（b）中該區(qū)域b值約為0.7。2010 年玉樹MS7.4 地震和2011 年囊謙MS5.2 地震均發(fā)生于該區(qū)域，兩者差異可能由是否剔除這2 次地震余震引起。圖5（b）、圖6（b）中瑪多MS7.4 地震震中位置空白，震中東部存在低b值區(qū)域，區(qū)域面積較小，b值在0.6 附近。

圖5（c）、圖6（c）為2013 年1 月1 日?2016 年12 月31 日去除余震目錄及完整目錄b值空間掃描結(jié)果，兩者沒有差異。

圖5 去余震目錄b 值空間分布圖像Fig. 5 Spatial distribution of b value with deleting aftershock

圖5（d）、圖6（d）為2017 年1 月1 日?2021 年4 月30 日去除余震目錄及完整目錄b值空間掃描結(jié)果，兩者在瑪多震中東南部存在差別，差異區(qū)域與圖5（a）、圖6（a）一致。造成該差異的主要原因?yàn)?020 年石渠MS5.9 地震發(fā)生于此區(qū)域，此次地震余震較多，去除余震會引起差異區(qū)域b值發(fā)生變化，因此該區(qū)域是否去除余震b值空間掃描結(jié)果呈現(xiàn)較為明顯的差異性。

圖6 完整目錄b 值空間分布圖像Fig. 6 Spatial distribution of b value without deleting aftershock

圖5（b）、（c）、（d）顯示隨著時間推移，瑪多震中附近低b值區(qū)域范圍逐漸擴(kuò)大，由正東不斷向震中擴(kuò)展。圖5（b）中低b值區(qū)域位于東昆侖斷裂帶與瑪多-甘德斷裂帶之間，圖5（d）中低b值區(qū)域遷移至瑪多-甘德斷裂帶，呈現(xiàn)南移現(xiàn)象。圖5（d）中低b值區(qū)域貫通震中東南區(qū)域，該區(qū)域b值逐年降低，因此存在發(fā)生強(qiáng)震的危險性。

對比2 種地震目錄在不同時間段的b值空間掃描結(jié)果，整體差異不大，但在低b值區(qū)域存在差異。該區(qū)域去余震目錄的b值相對較低，而完整地震目錄b值較高，這將掩蓋中強(qiáng)震區(qū)域的低b值特性。對比同一種地震目錄在不同時間段的b值空間掃描結(jié)果，可發(fā)現(xiàn)瑪多MS7.4 地震發(fā)生前低b值區(qū)域向震中的遷移現(xiàn)象，表明地震發(fā)生前震中附近應(yīng)力不斷集中。

4 結(jié)論及討論

本文根據(jù)中國地震臺網(wǎng)提供的2009 年1 月1 日?2021 年4 月30 日地震目錄，對瑪多MS7.4 地震震源區(qū)震前地震活動性進(jìn)行研究，分析區(qū)域內(nèi)中強(qiáng)震發(fā)生前b值的變化特征。通過對完整地震目錄和去余震地震目錄分別進(jìn)行b值時空掃描，發(fā)現(xiàn)是否剔除余震將影響中強(qiáng)震余震窗時間段內(nèi)b值結(jié)果，空間掃描結(jié)果表明，在低b值區(qū)域完整地震目錄b值高于剔除余震地震目錄b值，從而掩蓋了中強(qiáng)震區(qū)域低b值特性，而去余震地震目錄可以較好凸顯低b值區(qū)域，對中強(qiáng)震預(yù)測的指示考意義更強(qiáng)。因此，在進(jìn)行b值時空變化研究時，尤其是通過空間掃描結(jié)果研究中強(qiáng)震的位置時，應(yīng)剔除余震。對于文中結(jié)果整體差異較小的現(xiàn)象，主要是由于該地區(qū)最小完備震級較大，大量余震都處于最小完備震級以下，導(dǎo)致2 種地震目錄差異較小，因此所得結(jié)果差異也較小。雖然2 種地震目錄差異較小，但在較為關(guān)心的低b值區(qū)域仍存在一定差異，這將掩蓋中強(qiáng)震區(qū)域的低b值特性。因此，后續(xù)研究可選則地震監(jiān)測能力強(qiáng)、最小完備震級小的地區(qū)進(jìn)行去余震地震目錄與完整地震目錄對比的研究，驗(yàn)證是否隨著余震的增多，2 種地震目錄所得結(jié)果差異也隨之增大。具體結(jié)論如下：

（1）瑪多MS7.4 地震前1 年b值開始低于均值，且不斷下降，至b值有上升趨勢時地震發(fā)生，擴(kuò)大至區(qū)域內(nèi)其它5 級以上地震亦符合此規(guī)律，地震發(fā)生后b值明顯上升，短時間內(nèi)又下降至較低位置，并一直處于較低位置直至下次地震發(fā)生。這是由于地震發(fā)生前，該區(qū)域處于應(yīng)力積累狀態(tài)，b值處于較低位置，地震發(fā)生后應(yīng)力被釋放，區(qū)域內(nèi)的b值即恢復(fù)到正常水平。

（2）從b值空間掃描結(jié)果看，瑪多MS7.4 地震前震中位于低b值區(qū)域中心，該位置為b值最低處。

（3）不同時間段的b值空間掃描結(jié)果表明，瑪多MS7.4 地震發(fā)生前低b值區(qū)域向震中逐漸遷移，表明地震發(fā)生前應(yīng)力向震中附近不斷集中，因此低b值區(qū)域的遷移現(xiàn)象可為地震預(yù)報提供參考。

（4）余震剔除可以有效避免中強(qiáng)地震余震序列對b值結(jié)果的干擾，b值空間掃描時完整地震目錄掩蓋了中強(qiáng)震區(qū)域的低b值特性，而去余震地震目錄可以較好凸顯低b值區(qū)域，對中強(qiáng)震預(yù)測的指示意義更強(qiáng)。因此，在b值時空掃描研究時，應(yīng)剔除余震，力求更準(zhǔn)確的表現(xiàn)b值時空變化特征。