蘆山兩次強震序列活動特征及余震預(yù)測效能對比分析

畢金孟 宋程 馬永

摘要：利用可充分考慮小震信息的Omi-R-J模型，對2013年蘆山MS7.0地震和2022年蘆山MS6.1地震序列活動特征進行對比分析，發(fā)現(xiàn)蘆山MS7.0地震序列較蘆山MS6.1地震序列完備性隨時間變化更加明顯，震后早期檢測能力較低，兩次地震序列穩(wěn)定時段檢測得到的完整性震級均為1.8級左右。蘆山兩次地震序列p值差異較小，展示了相對正常的衰減過程；蘆山MS7.0地震序列k值明顯大于蘆山MS6.1地震序列，這或許與蘆山MS7.0地震較為發(fā)育的余震或強余震有關(guān)；蘆山MS7.0地震序列b值小于蘆山MS6.1地震序列，表明MS7.0地震之后蘆山地區(qū)仍處于較高的應(yīng)力狀態(tài)，而MS6.1地震之后應(yīng)力處于相對較低的水平。利用N-test方法開展效能評估，結(jié)果顯示，兩次地震序列初期階段余震發(fā)生率展示了較好的預(yù)測效能，基于震后初期1天內(nèi)的數(shù)據(jù)開展的未來1天的余震預(yù)測中，僅有1次預(yù)測失效，而在對未來3 d的余震預(yù)測中，不存在預(yù)測失效的情況。

關(guān)鍵詞：蘆山地震；Omi-R-J模型；檢測率；序列參數(shù)；余震預(yù)測效能評估

中圖分類號：P315.725?? 文獻標(biāo)識碼：A?? 文章編號：1000-0666（2023）02-0204-12

doi：10.20015/j.cnki.ISSN1000-0666.2023.0021

0 引言

余震是一個地震序列中的重要組成部分，余震發(fā)生率在震后初期呈現(xiàn)指數(shù)級的增長。開展震后初期的余震序列活動分析和強余震預(yù)測不僅能夠降低震后風(fēng)險，服務(wù)應(yīng)急決策，還可檢驗地震的可預(yù)測性（Reasenberg，Jones，1989；Gerstenberger et al，2005；Nanjo et al，2012；Ogata et al，2013；Omi et al，2019）。震后早期序列參數(shù)是震源區(qū)應(yīng)力調(diào)整、余震序列衰減、主震破裂、斷層愈合等地震活動特征在一定時空范圍內(nèi)的集中反映，如b值直接關(guān)系到地殼中應(yīng)力的分布狀態(tài)（Utsu，1961），p值與地殼內(nèi)部的結(jié)構(gòu)異質(zhì)性、溫度以及斷裂帶的性質(zhì)有關(guān)（Utsu et al，1995；Enescu et al，2011）。地震序列參數(shù)是開展余震預(yù)測的關(guān)鍵基礎(chǔ)，如Reasenberg和Jones（1989）根據(jù)模型參數(shù)的中值構(gòu)建了“普通加州模型”的概率預(yù)測模型，用于震后早期余震的概率預(yù)測。

目前國際上正在開展“地震可預(yù)測性合作研究”計劃（Collaboratory for the Study of Earthquake Predictability，CSEP）和“可操作的余震預(yù)測”（Operational Aftershock Forecasting，OAF），短期余震概率預(yù)測模型由此得到快速發(fā)展（Schorlemmer et al，2018）?；赗-J模型（Reasenberg，Jones，1989）和ETAS模型（Ogata，1989）可以利用震后幾天的數(shù)據(jù)開展較為有效的預(yù)測工作，但由于震后初期大量余震的缺失，很難在震后極短時間內(nèi)進行有效預(yù)測（Ogata，1983；Utsu et al，1995；Iwata，2008；Peng，Zhao，2009；Lippiello et al，2019），有效的余震概率預(yù)測只能在主震后24小時之后開始（Omi et al，2019），嚴(yán)重影響了震后初期預(yù)測的實效性。針對此類問題，Omi等（2013）充分考慮小震信息，在R-J模型的基礎(chǔ)上開發(fā)了一種Omi-R-J方法，并用于其檢驗2011年日本“3·11”大地震以及2017年Karaburun-Lesvos 6.3級地震，利用震后一天包括小震在內(nèi)的全部數(shù)據(jù)成功預(yù)測了未來一周內(nèi)的余震（Utkucua et al，2021）。以此為基礎(chǔ)，Omi 等（2016，2019）在日本建立了實時的可操作余震預(yù)測系統(tǒng)，并逐步發(fā)展和完善。對于實時預(yù)測的探索也一直在進行中（Ogata，Katsura，2006；Omi et al，2013），地震學(xué)家已在新西蘭（Gerstenberger et al，2016）、意大利（Marzocchi et al，2014）、日本（ERC，2016）和美國（Michael et al，2020）針對余震預(yù)測（Reverso et al，2018；Omi et al，2019）開展了前瞻性研究。

近年來龍門山斷裂發(fā)生多次強震，造成了巨大的人員傷亡和財產(chǎn)損失。為了探究龍門山斷裂南段地震序列的活動特征和可預(yù)測屬性，本文聚焦震后初期的地震序列活動特點，利用可充分考慮小震震級記錄技術(shù)優(yōu)勢的Omi-R-J模型，對2013年蘆山MS7.0和2022年蘆山MS6.1兩次強震進行地震序列活動分析和初期強余震預(yù)測等對比研究，以期為該地區(qū)“可操作性”的余震預(yù)測研究提供參考，為防震減災(zāi)工作服務(wù)。

1 研究背景和數(shù)據(jù)

龍門山斷裂帶地處巴顏喀拉地塊東邊界，受板塊推擠作用形成，構(gòu)造活動復(fù)雜多樣，強震頻發(fā)，其北段和南段以及塔藏斷裂等處接連發(fā)生了2008年5月12日汶川MS8.0、2013年4月20日蘆山MS7.0和2017年8月8日九寨溝MS7.0等強震。對于2013年蘆山MS7.0地震，由于反演的多解性或不確定性以及地表調(diào)查未發(fā)現(xiàn)明顯沿區(qū)域斷裂發(fā)育的同震破裂，目前包括基底滑脫帶（李傳友等，2013；張岳橋等，2013）、未知盲逆斷層（徐錫偉等，2013a，b）、山前隱伏斷裂（大邑斷裂）（李勇等，2013；蘇金蓉等，2013；周榮軍等，2013）等在內(nèi)的多種斷裂或構(gòu)造都被認為是該次地震可能的發(fā)震構(gòu)造。

根據(jù)房立華等（2013）的雙差重定位結(jié)果，2013年蘆山MS7.0地震的發(fā)震斷層面向NW傾斜，淺部傾角較陡，深部略緩，表現(xiàn)為“鏟形”逆沖斷層的特征，且主破裂面逆沖前段受阻反向逆沖出現(xiàn)一條向SE傾斜的反沖斷層，與發(fā)震斷層相交構(gòu)成“Y”字形。2022年蘆山MS6.1地震震中位于2013年蘆山MS6.1地震震中西北約10 km處。從余震震源深度分布來看，2022年蘆山地震的發(fā)震斷層和2013年蘆山地震明顯不同、向SE傾斜，可能是因為2013年蘆山地震“Y”型斷層系統(tǒng)中的反沖斷層發(fā)生了活動（中國地震局地球物理研究所，2022）。盡管兩次地震震中相距較近，但發(fā)生在不同的斷層上，顯示了該區(qū)域較為復(fù)雜的構(gòu)造特征。

為對蘆山兩次強震的序列活動特征開展研究，本文使用了全國統(tǒng)一正式編目地震目錄系統(tǒng)（中國地震臺網(wǎng)中心，2022），并選用兩次地震震后10 d內(nèi)的數(shù)據(jù)。2013年蘆山MS7.0地震后10 d內(nèi)共記錄到4 855次地震，其中5級以上余震4次，4.0～4.9級地震45次，3.0～3.9級地震244次；2022年蘆山MS6.1地震后10 d內(nèi)共記錄到916次余震，其中4級以上地震1次，3.0～3.9級地震7次。兩次地震序列的震級-時間（M-t）分布如圖1a-1、b-1所示，震級-序號圖如圖1a-2、b-2所示。由此可知，兩次蘆山地震，特別是2013年MS7.0地震震后早期余震事件缺失嚴(yán)重。

2 地震序列檢測率的變化特征

Reasenberg和Jones（1989）以余震衰減關(guān)系Omori-Utsu公式（Omori，1894；Utsu，1961）作為頻次約束，以震級-頻度分布關(guān)系式（Gutenberg，Richter，1944）作為強度約束，提出了余震短期概率預(yù)測模型。地震序列中，在t時刻震級不低于M的余震強度函數(shù)可表示為：

λ（t，M）=k（t+c）p×10-bM（1）

式中：t為余震的發(fā)震時刻；k控制余震及強余震活躍程度，與余震數(shù)目成正比；p的大小與序列的衰減快慢呈正相關(guān)關(guān)系；c用于調(diào)節(jié)主震發(fā)生后極短時間內(nèi)余震記錄的不完整性，與震源深度成負相關(guān)關(guān)系（Shebalin，Narteau，2017）；b與應(yīng)力累積水平呈反比關(guān)系（Wiemer，Katsumata，1999；Enescu et al，2011）。

為將地震序列早期階段完整性震級以下的余震考慮到模型參數(shù)擬合及后續(xù)的余震發(fā)生率預(yù)測中，Omi等（2013）在R-J模型（Reasenberg，Jones，1989）的基礎(chǔ)上發(fā)展了“Omi-R-J”模型，即利用Ogata和Katsura（1993）給出的檢測率函數(shù)q（M），來描述地震事件記錄不完整部分的檢測程度。實際記錄到的地震概率密度函數(shù)可表示為：

P（Mβ，μ，σ）=e-βMq（Mμ，σ）∫+∞-∞e-βMq（Mμ，σ）dM=e-βMq（Mμ，σ）e-βμ+β2 σ2/2/β=βe-β（M-μ）+β2 σ2/2q（Mμ，σ）（2）

式中：μ表示在檢測率為50%時對應(yīng)的震級；σ（σ>0）表示部分地震事件被檢測到的震級范圍；e-βM為震級-頻度分布（FMD）中的冪指數(shù)函數(shù)，其中M為震級，β與G-R關(guān)系中的b值成線性關(guān)系，即β=bln10。實際記錄到的余震的檢測率隨時間呈現(xiàn)一種動態(tài)變化，由于震后初期余震的缺失，使得μ值在地震序列的早期階段是一個相對高值，隨著序列持續(xù)時間的增加，逐漸恢復(fù)到正常水平，μ為時間t的函數(shù)，即μ（t）=μi（ti-1

Pβ，σ（Mμ）=∏ni=1βe-β（Mi-μi）+β2 σ2/2q（Miμi，σ）（3）

對于其平滑程度，需要利用參數(shù)V來控制μ（t），其通過懲罰的二階差分來進行平滑，可表示為：

PV（μ）=P（μ1，μ2）∏n-2i=112πVe-（μi+2-2μi+2+μi）22V（4）

同時令μ1和μ2均符合均勻先驗分布，即P（μ1，μ2）=const。根據(jù)貝葉斯準(zhǔn)則，后驗函數(shù)可表示為：

Pβ，σ，V（μM）=Pβ，σ（Mμ）PV（μ）∫Pβ，σ（Mμ′）PV（μ′）dμ′（5）

Omi等（2013）利用牛頓迭代法（Newtons iteration）對超參數(shù)β、σ和V進行估計。在地震參數(shù)估計中，如果采取多次重復(fù)的計算方式，很難計算出參數(shù)估計中的偏差，獲得較為穩(wěn)定、可靠的參數(shù)估計值。針對數(shù)據(jù)不完整性問題。Veen和Schoenberg（2008）引入期望最大化（expectation-maximization，EM）算法來獲得參數(shù)估計結(jié)果，本文也借鑒最大期望算法對參數(shù)進行優(yōu)化。

Omi-R-J方法通過對每個地震事件逐步迭代的方式，可獲得研究時段內(nèi)的地震檢測率隨時間的變化。由于基于OK1993模型（Ogata，Katsura，1993）的Omi-R-J模型的μ是連續(xù)估計的參數(shù)，可分別用95.45%置信水平下完整記錄的μ+2σ或者在99.74%置信水平下完整記錄的μ+3σ，來近似描述最小完整性震級MC（Mignan，Woessner，2012；Iwata，2013）。圖2a-1、b-1分別給出了2013年蘆山MS7.0地震和2022年蘆山MS6.1地震序列按照Ogata和Katsura（1993）模型的50%地震檢測率μ（t）（50%）以及μ（t）+2σ（95.45%），μ（t）+3σ（99.74%）的分布情況。由圖2a-1中的3個范圍的檢測率可知，2013年蘆山地震序列的完備性隨時間變化較為明顯，震后缺失事件較多，余震記錄事件有限，檢測能力較低，在震后2 d逐漸恢復(fù)至穩(wěn)定階段，如果以μ（t）+3σ作為最小完整性震級，則被檢測到余震的震級約為1.8級。2022年蘆山地震序列的完備性隨時間變化比2013年蘆山地震序列平穩(wěn)，如圖2b-1所示，初期階段μ（t）+3σ檢測到余震的震級為2.5級，然后逐漸過渡至平穩(wěn)時段的1.8級。造成蘆山兩次地震早期余震序列檢測結(jié)果差異的原因為：一是2013年蘆山地震震級較大，面波發(fā)育，震后早期更多的地震事件被“淹沒”在主震中，較多的小地震無法被及時識別出來；二是近幾年理記技術(shù)的發(fā)展和地震臺網(wǎng)布局的優(yōu)化使得更多的小震、微震被檢測到，震后早期地震序列的檢測率逐步提升。為更加定量化地分析檢測率，本文基于0～10.00 d的觀測數(shù)據(jù)，利用震級頻度變化（FMD）進行了OK1993模型參數(shù)擬合，圖2a-2、b-2分別給出了兩次地震的擬合結(jié)果，菱形為地震目錄的震級-頻度分布，實線為OK1993模型擬合結(jié)果。擬合所得兩組參數(shù)分別為β=1.250 9、μ=0.727 2、σ=0.278 5和β=2.180 9、μ=1.099 1、σ=0.320 7，擬合結(jié)果較好，驗證了檢測率變化的可靠性。

3 地震序列參數(shù)及不確定性分析

可靠穩(wěn)定的地震序列參數(shù)是開展可操作的余震預(yù)測和地震危險性評估等研究的基礎(chǔ)。Omi-R-J方法利用記錄的全部余震事件，在震后短期內(nèi)快速開展參數(shù)擬合，獲得較為可靠的地震序列參數(shù)。通過利用狀態(tài)—空間模型獲得參數(shù)β、σ和V以及檢測率的動態(tài)變化μ（t）后，考慮不完整地震記錄的地震發(fā)生率函數(shù)可以表示為ν（t，M）=λ（t，M）q（Mμ（t），σ），與k、p、c相關(guān)的對數(shù)似然函數(shù)表示為：

lnL（k，c，p）=∑Mi≥Mclnν（ti，Mi）-∫∞McdM∫0Tdtν（t，M）（6）

式中：ti和Mi分別為在模型擬合的“學(xué)習(xí)時段” ［0，T］內(nèi)發(fā)生的第i個余震的發(fā)震時刻與震級。

利用后驗分布對Omi-R-J模型參數(shù)的不確定性進行檢驗、評估，參數(shù)的邊際后驗分布范圍越寬（越窄），表示不確定性越大（越小）。為直觀顯示后驗分布的離散圖，采用馬爾可夫鏈蒙特卡羅方法（Markov Chain Monte Carlo，MCMC）從后驗分布中抽取Omi-R-J參數(shù)集，其優(yōu)點是能夠以幾乎精準(zhǔn)的方式評估不確定性，而無需對后驗函數(shù)進行近似。采用Metropolis算法（Metropolis et al，1953）模擬一個多變量馬爾可夫鏈，在高維后驗分布中，可隨機選擇一個分量來進行轉(zhuǎn)換。由于Omi-R-J模型參數(shù)的維數(shù)較少，同時創(chuàng)建所有參數(shù)分量來轉(zhuǎn)換。為保證模型參數(shù)的獨立性，需要重復(fù)模擬過程，直到平穩(wěn)馬爾可夫鏈的樣本足夠長，然后在時間間隔內(nèi)對參數(shù)進行采樣。本文基于Metropolis算法模擬10 000個參數(shù)集，然后選取1 000套（t=10，20，30，…，10 000）估算Omi-R-J模型參數(shù)。圖3給出了2022年蘆山MS6.1地震的參數(shù)估計值（黑十字）、采樣參數(shù)集（灰點）以及采樣參數(shù)的統(tǒng)計分布情況（柱狀圖），揭示了后驗分布的復(fù)雜形式，以及參數(shù)之間的統(tǒng)計相關(guān)性。各數(shù)據(jù)采樣點參數(shù)呈現(xiàn)出較好的優(yōu)勢分布，邊際后驗分布范圍較窄，即存在較小的不確定性，揭示了參數(shù)擬合結(jié)果的可靠性。此外，各參數(shù)間存在一定的相關(guān)性，如反映應(yīng)力水平的b（β=bln10）值與余震數(shù)目的k值具有較強的負相關(guān)，從另一側(cè)面說明了應(yīng)力水平較高的區(qū)域更易觸發(fā)較多的強余震。

為更好地描述參數(shù)的變化特征，本文對蘆山兩次強震序列分別采用了連續(xù)、滑動的參數(shù)擬合，擬合起始時間為主震后0.05 d，以0.05 d為步長，持續(xù)增加至序列的10.00 d，進行了100個時段的滑動擬合。圖4給出了擬合所得p值、c值、k值和b值隨序列持續(xù)時間的變化圖。由圖可知，各參數(shù)值在序列的早期階段變化相對劇烈，后期逐漸平穩(wěn)。與基于ETAS模型所得的2013年蘆山地震序列參數(shù)相比（蔣長勝等，2013），Omi-R-J模型在震后早期能夠更早更快地獲得較為穩(wěn)定的序列參數(shù)。此外Omi-R-J方法曾應(yīng)用于2008年四川汶川MS8.0地震（Bi，Jiang，2020a）、2017年四川九寨溝MS7.0地震（蔣長勝等，2018）、華北地區(qū)（Bi，Jiang，2020b）以及中國大陸其它板內(nèi)地震（畢金孟等，2022），并與ETAS模型和R-J模型做了對比分析，發(fā)現(xiàn)Omi-R-J模型在震后早期能更快獲得穩(wěn)定參數(shù)。穩(wěn)定時段，2013年四川蘆山MS7.0地震的模型參數(shù)為p=1.017±0.026，c=0.014±0.005，k=0.029±0.006，b=0.784±0.015；2022年四川蘆山MS6.1地震的模型參數(shù)為p=1.079±0.032，c=0.022±0.007，k=0.002±0.001，b=0.942±0.049，其中b值根據(jù)β=ln10×b所得。

由圖4a-1、b-1可見，兩個序列的p值相差不大，2013年蘆山MS7.0地震序列在震后5 d內(nèi)呈現(xiàn)一定的變化幅度，之后逐漸穩(wěn)定在1.017左右，略低于基于ETAS模型的計算結(jié)果1.22（蔣長勝等，2013）；2022年蘆山MS6.1地震序列p值在震后3 d內(nèi)變化較為劇烈，之后逐漸穩(wěn)定在1.079左右，較2013年蘆山MS7.0地震序列變化更為劇烈；兩個地震序列的p值與以中國大陸6.0級以上地震序列計算的平均值（0.984 2）和中值（1.020 2）基本接近（畢金孟等，2022），展示了相對正常的序列衰減過程。

由圖4a-2、b-2可知，兩個序列的c值相差不大，2013年蘆山MS7.0地震序列穩(wěn)定值趨向于0.014附近，2022年蘆山MS6.1地震序列穩(wěn)定值趨向于0.022附近，略大于2013年蘆山MS7.0地震序列穩(wěn)定值。由于壓力和溫度通常都隨著深度的增加而增加，c隨深度的增加而減?。⊿hebalin， Narteau，2017）?；趨^(qū)域地震臺網(wǎng)記錄，利用“裁剪—粘貼”（CAP）方法對兩次蘆山地震進行矩心矩張量反演得到2013年蘆山MS7.0地震的震源矩心深度約為14 km（呂堅等，2013），2022年蘆山MS6.1地震的最佳矩心震源深度約16 km（梁皓等，2022）。同樣，基于雙差定位的方法測得兩次地震事件的震源深度分別為17.6 km和20 km（房立華等，2013；中國地震局地球物理研究所，2022），兩種方法測得兩次地震震源深度略有不同，c值的差異也證明了兩次地震的震源深度存在一定的差異。

由圖4c-1、c-2可見，兩個序列的k值存在很大的差異，2013年蘆山MS7.0地震的地震序列k值在震后早期階段（4 d內(nèi)）變化較為劇烈，隨后緩慢減小至0.029附近；2022年蘆山MS6.1地震序列R值在震后早期階段（3 d內(nèi)）基本穩(wěn)定在0.002附近，較2013年蘆山MS7.0地震地震序列值明顯偏小。k值控制余震的活躍程度，在中國大陸板內(nèi)地震中與強余震數(shù)量呈正相關(guān)關(guān)系（畢金孟等，2022），這很大程度上取決于單個余震序列，且與主震震級關(guān)系不大，這也解釋了2013年蘆山MS7.0地震較為發(fā)育的余震序列和2022年蘆山MS6.1地震序列余震尤其是強余震數(shù)量相對較少的原因。

由圖4d-1、d-2可見，兩個序列的b值存在一定的差異，在 震后初期階段（2 d內(nèi)）呈現(xiàn)出一定的變化，之后2013年蘆山MS7.0地震序列b值逐漸趨于0.784附近，與基于最大似然法獲得的b值結(jié)果接近（蔣長勝等，2013），2022年蘆山MS6.1地震序列b值逐漸趨于0.942附近，明顯高于2013年蘆山MS7.0地震序列。2013年蘆山MS7.0地震序列的b值低于中國大陸的均值線（0.849 4）和中值線（0.830 6）（畢金孟等，2022），顯示出了余震區(qū)較高的應(yīng)力水平，而2022年蘆山MS6.1地震序列b值高于中國大陸的b值均值線和中值線，表明該地區(qū)的應(yīng)力水平已得到很大程度的釋放，也表明了該地區(qū)發(fā)生強余震的風(fēng)險大大降低。2013年蘆山MS7.0地震序列b值呈現(xiàn)出先上升后下降的過程，反映了主震后應(yīng)力得到釋放，然后繼續(xù)積累，繼而觸發(fā)余震的過程，而2022年蘆山MS6.1地震序列b值呈現(xiàn)出先下降后逐漸上升最終趨于穩(wěn)定的過程，反映了震后短暫的應(yīng)力調(diào)整后就持續(xù)釋放的過程。

4 余震預(yù)測結(jié)果及效能評估

基于上述的序列參數(shù)擬合結(jié)果，通過式（1）計算在M>MP和任意時間間隔［t2，te］內(nèi)的余震預(yù)測數(shù)目：

N=∫tet2∫∞Mpλ（t，M）dtdM（7）

為檢驗蘆山兩次強震在震后初期（1 d內(nèi)）的“M-3”強余震發(fā)生率的預(yù)測效能，即距主震震級3個震級檔以內(nèi)的余震事件視為強余震，選用CSEP計劃中針對地震數(shù)檢驗的N-test方法（Kagan，Jackson，1995；Schorlemmer et al，2007；Zechar，2010）來驗證Omi-R-J模型預(yù)測的地震數(shù)目與實際數(shù)目的偏離程度，并利用有效顯著性水平αeff=0.025對N-test檢驗結(jié)果進行單邊檢驗，當(dāng)δ1<αeff時，表示預(yù)測強余震數(shù)目過少；當(dāng)δ2<αeff時則表示預(yù)測強余震數(shù)目過多。

以蘆山兩次地震后0.20～1.20 d、0.20～3.20 d的強余震預(yù)測為例，圖5a分別給出了基于Omi-R-J模型對蘆山兩次地震未來1 d（圖5a-1、a-2）、3 d（圖5a-3、a-4）預(yù)測結(jié)果與實際觀測結(jié)果在震級-累積頻次上的比較，圖中黑色圓點為實際觀測結(jié)果，黑色線段為預(yù)測結(jié)果，灰色陰影為預(yù)測結(jié)果的95%置信區(qū)間。圖5b分別給出了蘆山兩次地震未來1 d（圖5b-1、b-2）、3 d（圖5b-3、b-4）距主震震級3個震級檔的預(yù)測結(jié)果與實際觀測結(jié)果在時間-累積頻次上的比較，圖中黑色曲線為實際觀測結(jié)果，灰色曲線為預(yù)測結(jié)果，灰色虛線標(biāo)出了預(yù)測結(jié)果的95%置信區(qū)間，黑色垂直虛線標(biāo)出了預(yù)測起始時刻。實際發(fā)生地震數(shù)均在預(yù)測數(shù)的95%置信區(qū)間，其評分值依次為［0.821 6，0.256 1］、［0.425 1，0.758 0］、［0.830 8，0.227 7］、［0.748 4，0.423 4］，顯示了較好的預(yù)測效能。由于大量的強余震發(fā)生在主震后的數(shù)天內(nèi)，尤其是震后1～3 d內(nèi)，是減輕地震災(zāi)害風(fēng)險的“黃金時間”，因此，此次僅針對震后1 d內(nèi)數(shù)據(jù)分別在［0～0.2，0～0.4，0～0.6，0～0.8，0～1.0］時段開展震后1 d和3 d的強余震（M>3）發(fā)生率預(yù)測及“量化”的效能評估，其結(jié)果見表1，在不考慮預(yù)測時段未發(fā)生地震的情況下，僅有一次地震事件未通過N-Test檢驗，出現(xiàn)預(yù)測過少的情況，其他時段均表現(xiàn)出了較好的預(yù)測效能。

5 結(jié)論

為分析2013年蘆山MS7.0和2022年蘆山MS6.1兩次強震的序列活動特征和震后初期的強余震預(yù)測效能，本文基于Omi-R-J模型，利用余震序列早期階段大量不完整的小震級事件，對序列的檢測率變化和序列參數(shù)的變化進行了對比研究，并在震后初期（1 d內(nèi)）開展強余震發(fā)生率預(yù)測以及利用針對地震數(shù)的N-test方法進行效能評估，得到以下結(jié)果：

（1）2013年蘆山MS7.0地震序列的完備性隨時間變化明顯，震后早期階段余震記錄事件有限，檢測能力較低，此后兩個地震序列的完備性逐漸恢復(fù)，至穩(wěn)定階段μ（t）+3σ（99.73%）檢測到的余震的震級約為1.8。

（2）從總體上看，由于地震震級大小、發(fā)震構(gòu)造等方面的不同，蘆山兩次地震的參數(shù)變化以及穩(wěn)定時段的參數(shù)值均存在一定的差異。在震后早期階段，2013年蘆山MS7.0地震序列參數(shù)變化較為劇烈且達到穩(wěn)定時段所需時間較長，這或許是與其序列較復(fù)雜有關(guān)。2013年蘆山MS7.0地震穩(wěn)定時段的序列參數(shù)為p=1.017±0.026，c=0.014±0.005，k=0.029±0.006，b=0.784±0.015；2022年蘆山MS6.1地震序列參數(shù)為p=1.079±0.032，c=0.022±0.007，k=0.002±0.001，b=0.942±0.049。

（3）對比兩次蘆山地震序列參數(shù)發(fā)現(xiàn)，p值相差不大，且與中國大陸6.0級以上地震序列的平均值和中值接近，展示了相對正常的序列衰減過程；兩次地震的震源深度的不同使得c值存在略微差異；兩次地震序列的k值存在很大差異——2013年蘆山MS7.0地震序列的k值明顯大于2022年蘆山MS6.1地震序列，這與2013年蘆山地震較為發(fā)育的余震或強余震有關(guān)；兩次地震序列的b值存在一定的差異，2013年蘆山MS7.0地震序列的b值小于2022年的蘆山MS6.1地震序列，表明2013年之后蘆山地區(qū)仍是高應(yīng)力地區(qū)，依然存在發(fā)生中強震或強余震的風(fēng)險，而2022年蘆山MS6.1地震之后，應(yīng)力得到釋放，應(yīng)力水平處于一個相對較低的狀態(tài)。

（4）利用N-test方法對蘆山兩次地震震后初期強余震預(yù)測結(jié)果的效能開展評估，結(jié)果顯示，Omi-R-J模型對兩個地震序列震后初期階段展示了較好的預(yù)測效能，在基于震后1 d內(nèi)的數(shù)據(jù)開展的分段和分尺度預(yù)測中，兩個序列對未來1 d預(yù)測，僅存在一次預(yù)測失效的情況，而對未來3 d的預(yù)測，不存在預(yù)測失效的情況，這證明了Omi-R-J模型可以在震后初期階段就能開展較為有效的余震預(yù)測。

本文研究中使用了中國地震臺網(wǎng)中心“全國地震編目系統(tǒng)”提供的“統(tǒng)一正式目錄”，中國地震局地球物理研究所蔣長勝研究員和日本東京大學(xué)生產(chǎn)技術(shù)研究所Takahiro Omi博士為本研究提供了程序和技術(shù)支持，評審專家提出了諸多建設(shè)性修改建議，對稿件質(zhì)量提升幫助很大，在此一并表示感謝。

參考文獻：

畢金孟，蔣長勝，來貴娟，等.2022.中國大陸強震的早期余震概率預(yù)測效能評估與制約因素［J］.地球物理學(xué)報，65（7）：2532-2545.

房立華，吳建平，王未來，等.2013.四川蘆山MS7.0級地震及其余震序列重定位［J］.科學(xué)通報，58（20）：1901-1909.

蔣長勝，畢金孟，王福昌，等.2018.利用早期余震預(yù)測的Omi-R-J方法對2017年四川九寨溝MS7.0地震的應(yīng)用研究［J］.地球物理學(xué)報，61（5）：2099-2110.

蔣長勝，莊建倉，龍鋒，等.2013.2013年蘆山MS7.0地震序列參數(shù)的早期特征：傳染型余震序列模型計算結(jié)果［J］.地震學(xué)報，35（5）：661-669.

李傳友，徐錫偉，甘衛(wèi)軍，等.2013.四川省蘆山MS7.0地震發(fā)震構(gòu)造分析［J］.地震地質(zhì)，35（3）：671-683.

李勇，周榮軍，趙國華，等.2013.龍門山前緣的蘆山地震與逆沖—滑脫褶皺作用［J］.成都理工大學(xué)學(xué)報（自然科學(xué)版），40（4）：353-363.

梁皓，戴丹青，楊志高，等.2022.2022年6月1日四川蘆山6.1級地震應(yīng)急產(chǎn)品及震源參數(shù)初步分析［J］.中國地震，38（2）：360-369.

呂堅，王曉山，蘇金蓉，等.2013.蘆山7.0級地震序列的震源位置與震源機制解特征［J］.地球物理學(xué)報，56（5）：1753-1763.

蘇金蓉，鄭鈺，楊建思，等.2013.2013年4月20日四川蘆山M7.0級地震與余震精確定位及發(fā)震構(gòu)造初探［J］.地球物理學(xué)報，56（8）：2636-2644.

徐錫偉，陳桂華，于貴華，等.2013a.蘆山地震發(fā)震構(gòu)造及其與汶川地震關(guān)系討論［J］.地學(xué)前緣，20（3）：11-20.

徐錫偉，聞學(xué)澤，韓竹軍，等.2013b.四川蘆山7.0級強震：一次典型的盲逆斷層型地震［J］.科學(xué)通報，58（20）：1887-1893.

張岳橋，董樹文，侯春堂，等.2013.四川蘆山2013年MS7.0地震發(fā)震構(gòu)造初步研究［J］.地質(zhì)學(xué)報，87（6）：747-758.

中國地震局地球物理研究所.2022.2022年6月1日四川蘆山6.1級地震的震源特征［EB/OL］.（2022-06-02）［2022-06-28］.https：//www.cea-igp.ac.cn/kydt/279024.html.

中國地震臺網(wǎng)中心.2022.全國地震編目系統(tǒng)地震目錄［DB/OL］.（2022-06-25）［2022-06-25］.http：//10.5.160.18/console/exit.action.

周榮軍，李勇，蘇金蓉，等.2013.四川蘆山MW6.6級地震發(fā)震構(gòu)造［J］.成都理工大學(xué)學(xué)報（自然科學(xué)版），40（4）：364-370.

Bi J M，Jiang C S.2020a.Comparison of early aftershock forecasting for the 2008 Wenchuan MS8.0 earthquake［J］.Pure and Applied Geophysics，177（1）：9-25.

Bi J M，Jiang C S.2020b.Research on the forecasting strategy of early aftershocks in North China［J］.Annals of Geophysics，63（4）：SE441.

Enescu B，Enescu D，Ito K.2011.Values of b and p：their variations and relation to physical processes for earthquakes in Japan and Romania［J］.Romanian Journal of Physics，56（3-4）：590-608.

ERC（Earthquake Research Committee report）.2016.Information on earthquake forecasting after a large earthquake［J/OL］.（2016-08-19）［2022-07-02］.https：//www.jishin.go.jp/main/yosoku_info/honpen.pdf.

Gerstenberger M C，Rhoades D A，McVerry G H.2016.A hybrid time-dependent probabilistic seismic-hazard model for Canterbury，New Zealand［J］.Seismological Research Letters，87（6）：1311-1318.

Gerstenberger M C，Wiemer S，Jones L M，et al.2005.Real-time forecasts of tomorrows earthquakes in California［J］.Nature，435（7040）：328-331.

Gutenberg R，Richter C F.1944.Frequency of earthquakes in California［J］.Bulletin of the Seismological Society of America，34（4）：185-188.

Iwata T.2008.Low detection capability of global earthquakes after the occurrence of large earthquakes：Investigation of the Harvard CMT catalogue［J］.Geophysical Journal International，174（3）：849-856.

Kagan Y Y，Jackson D D.1995.New seismic gap hypothesis：five years after［J］.Journal of Geophysical Research：Solid Earth，100（B3）：3943-3959.

Lippiello E，Cirillo A，Godano C，et al.2019.Post seismic catalog incompleteness and aftershock forecasting［J］.Geosciences，9（8）：355.

Marzocchi W，Lombardi A M，Casarotti E.2014.The establishment of an operational earthquake forecasting system in Italy［J］.Seismological Research Letters，85（5）：961-969.

Metropolis N，Rosenbluth A W，Rosenbluth M N，et al.1953.Equation of state calculations by fast computing machines［J］.The Journal of Chemical Physics，21（6）：1087-1092.

Michael A J，McBride S K，Hardebeck J L，et al.2020.Statistical seismology and communication of the USGS operational aftershock forecasts for the 30 November 2018 MW7.1 Anchorage，Alaska，earthquake［J］.Seismological Research Letters，91（1）：153-173.

Mignan A，Woessner J.2012.Estimating the magnitude of completeness for earthquake catalogs［J/OL］.（2012-05）［2022-07-02］Community Online Resource for Statistical Seismicity Analysis，1-45.http：//www.corssa.org/export/sites/corssal.galleries/articlespaf/Migan-Woessner-2012-CORSSA-Magitude-of-completeness.pdf-2063069299.pdf.

Nanjo K Z，Tsuruoka H，Yokoi S，et al.2012.Predictability study on the aftershock sequence following the 2011 Tohoku-Oki，Japan，earthquake：first results［J］.Geophysical Journal International，191（2）：653-658.

Ogata Y，Katsura K，F(xiàn)alcone G，et al 2013.Comprehensive and topical evaluations of earthquake forecasts in terms of number，time，space，and magnitude［J］.Bulletin of the Seismological Society of America，103（3）：1692-1708.

Ogata Y，Katsura K.1993.Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues［J］.Geophysical Journal International，113（3）：727-738.

Ogata Y，Katsura K.2006.Immediate and updated forecasting of aftershock hazard［J］.Geophysical Research Letters，33（10）：L10305.

Ogata Y.1983.Estimation of the parameters in the modified Omori Formula for aftershock frequencies by the maximum likelihood procedure［J］.Journal of Physics of the Earth，31（2）：115-124.

Ogata Y.1989.Statistical model for standard seismicity and detection of anomalies by residual analysis［J］.Tectonophysics，169（1-3）：159-174.

Omi T，Ogata Y，Hirata Y，et al.2013.Forecasting large aftershocks within one day after the main shock［J］.Scientific Reports，3（1）：2218.

Omi T，Ogata Y，Shiomi K，et al.2016.Automatic aftershock forecasting：a test using real-time seismicity data in Japan［J］.Bulletin of the Seismological Society of America，106（6）：2450-2458.

Omi T，Ogata Y，Shiomi K，et al.2019.Implementation of a real-time system for automatic aftershock forecasting in Japan［J］.Seismological Research Letters，90（1）：242-250.

Omori F.1894.On aftershocks of earthquakes.Journal of the College of Science［J］.Imperial University of Tokyo，7：111-200.

Peng Z G，Zhao P.2009.Migration of early aftershocks following the 2004 Parkfield earthquake［J］.Nature Geoscience，2（12）：877-881.

Reasenberg P A，Jones L M.1989.Earthquake hazard after a mainshock in California［J］.Science，243（4895）：1173-1176.

Reverso T，Steacy S，Marsan D.2018.A hybrid ETAS-Coulomb approach to forecast spatiotemporal aftershock rates［J］.Journal of Geophysical Research：Solid Earth，123（11）：9750-9763.

Schorlemmer D，Gerstenberger M C，Wiemer S，et al 2007.Earthquake likelihood model testing［J］.Seismological Research Letters，78（1）：17-29.

Schorlemmer D，Werner M J，Marzocchi W，et al 2018.The collaboratory for the study of earthquake predictability：Achievements and priorities［J］.Seismological Research Letters，89（4）：1305-1313.

Shebalin P，Narteau C.2017.Depth dependent stress revealed by aftershocks［J］.Nature Communications，8（1）：1317.

Utkucua M，Süleyman S，Nalbant B，et al.2021.The June 12，2017 M6.3 Karaburun-Lesvos earthquake of the Northern Aegean Sea：Aftershock forecasting and stress transfer［J］.Tectonophysics，814：228945.

Utsu T，Ogata Y，Matsuura R S.1995.The centenary of the Omori formula for a decay law of aftershock activity［J］.Journal of Physics of the Earth，43（1）：1-33.

Utsu T.1961.A statistical study of on the occurrence of aftershocks［J］.Geophysical Magnitude，30：521-605.

Veen A，Schoenberg F P.2008.Estimation of space-time branching process models in seismology using an em-type algorithm［J］.Journal of the American Statistical Association，103（482）：614-624.

Wiemer S，Katsumata K.1999.Spatial variability of seismicity parameters in aftershocks zones［J］.Journal of Geophysical Research，104（13）：13135-13151.

Zechar J D.2010.Evaluating earthquake predictions and earthquake forecasts：a guide for students and new researchers［J/OL］.（2010-09-01）［2022-06-29］.http：//www.corssa.org/export/sites/corssa/.galleries/articles-pdf/zechar.pdf.

Comparative Analysis of the Activity Characteristics and the AftershockForecasting Efficiency of Two Earthquake Sequences in Lushan

BI Jinmeng1，2，SONG Cheng1，MA Yong1

（1.Tianjin Earthquake Agency，Tianjin 300201，China）

（2.Institute of Geophysics，China Earthquake Administration，Beijing 100081，China）

Abstract

We used the Omi-R-J model，which can fully make use of the information of small earthquakes，to compare and analyze the characteristics of the activity of two strong-earthquake sequences in Lushan，Sichuan Province.we found that the completeness of the Lushan MS7.0 earthquake sequence in 2013 changed more obviously with time than that of the Lushan MS6.1 earthquake sequence in 2022.In the early period after the two mainshocks，a few aftershocks were detected by the Omi-R-J model，and the completeness magnitude detected in the stable period of the two earthquake sequences was both about 1.8.The difference between the p-values of the two Lushan earthquake sequences was small，showing a relatively normal attenuation process of the two sequences.k-value of the MS7.0 earthquake sequence was significantly larger than that of the MS6.1 sequence；this may be the result of the relatively more aftershocks or the strong aftershocks of the MS7.0 sequence.b-value of the MS7.0 earthquake sequence was smaller than that of the MS6.1 sequence；this indicated that Lushan area was still in a high-stress state after the 2013 MS7.0 earthquake，but after the 2022 MS6.1 earthquake the stress reached a relatively low level.The efficiency evaluation of the N-test method showed that in the initial stage of the two earthquake sequences the aftershock forecasting was efficient.Based on the next days data，the forecasting of the aftershocks on the second day after the mainshock only failed once，while the forecasting of the aftershocks of the successive 3 days after the second day of the mainshock did not fail.

Keywords：the Lushan earthquake；the Omi-R-J model；detection rate；sequence parameters；effectiveness evaluation of aftershork forecasting

收稿日期：2022-07-09.

基金項目：地震預(yù)測開放基金（XH23072D）；中國地震局震情跟蹤定向工作任務(wù)（2022010116，2020010104）.

第一作者簡介：畢金孟（1989-），工程師，主要從事地震活動性研究.E-mail：jinmengbi@126.com.

畢金孟，宋程，馬永.2023.蘆山兩次強震序列活動特征及余震預(yù)測效能對比分析［J］.地震研究，46（2）：204-215，doi：10.20015/j.cnki.ISSN1000-0666.2023.0021.