姜叢, 蔣長勝*, 尹鳳玲*, 張延保, 畢金孟, 龍鋒, 司政亞, 尹欣欣,5

1 中國地震局地球物理研究所, 北京 100081 2 天津市地震局, 天津 300201 3 四川省地震局, 成都 610041 4 北京市地震局, 北京 100080 5 甘肅省地震局, 蘭州 730000

0 引言

描述地震震級-頻度分布(FMD)的古登堡-里克特(G-R)關(guān)系中的b值，被證實(shí)與地殼中的差應(yīng)力負(fù)相關(guān)(Wyss, 1973; Mori and Abercrombie, 1997; Toda et al., 1998)，使其成為間接反映地震孕育發(fā)生時(shí)空物理過程的重要參量，因此在地震危險(xiǎn)性分析、地震工程和地球動力學(xué)研究中得到廣泛應(yīng)用.除了在地表二維平面、考慮深度的三維空間、時(shí)-空間研究b值的非均質(zhì)性分布和演化規(guī)律外(Wiemer and Wyss, 1997; Schorlemmer and Wiemer, 2005; Murru et al., 2007; Wang et al., 2021)，b值在地震時(shí)間序列中也同樣具有重要的應(yīng)用價(jià)值.這是由于前震/余震被限制在與主震破裂尺度相關(guān)的震區(qū)內(nèi)，干熱巖和頁巖氣開發(fā)等新型工業(yè)活動以及深部礦山開采等引起的誘發(fā)地震活動也被局限在注采區(qū)較小空間范圍內(nèi)集中發(fā)生，因此可簡化為只需要考慮時(shí)間序列b值變化，這極大降低了研究的復(fù)雜性.目前時(shí)間序列b值已在誘發(fā)地震的最大震級Mmax估計(jì)(Hallo et al., 2014; Van Der Elst et al., 2016)、前瞻型預(yù)測的風(fēng)險(xiǎn)管控“交通燈系統(tǒng)”的構(gòu)建(Ader et al., 2020)，以及天然地震的前震與主震識別(Van Der Elst, 2021)、強(qiáng)震前的前兆性變化(史海霞等, 2018; Xie et al., 2019)、后續(xù)強(qiáng)余震的風(fēng)險(xiǎn)分級(Gulia and Wiemer, 2019)等多場景應(yīng)用中發(fā)揮關(guān)鍵作用.

但時(shí)間序列b值的計(jì)算始終面臨多方面的挑戰(zhàn)，包括傳統(tǒng)采用的固定地震數(shù)目和時(shí)間長度的固定窗口法(Hutton et al., 2010; Ogata, 2011)在計(jì)算規(guī)則設(shè)置上過強(qiáng)的人為主觀性，以及由于地震序列在時(shí)間上的叢集性和稀疏性造成對b值的突變難于準(zhǔn)確識別，都使得時(shí)間序列b值的計(jì)算以及對震例的解釋存在較多爭議(Kamer and Hiemer, 2013).針對上述問題，F(xiàn)iedler 等(2018)發(fā)展了一種利用貝葉斯因子來選擇模型并識別時(shí)間序列b值突跳點(diǎn)的方法，但在震級-頻度關(guān)系上仍選用了G-R關(guān)系.針對主震后大量余震缺失造成的完備震級以上可使用地震數(shù)目過少的問題，Van Der Elst(2021)發(fā)展了稱為“b-positive”的時(shí)間序列b值計(jì)算方法，用連續(xù)地震之間震級差的正子集來估算b值并規(guī)避了此問題.然而這些計(jì)算方法尚且不能將兩類問題同時(shí)解決，想要得到廣泛的認(rèn)可仍有待更多的震例驗(yàn)證和方法改進(jìn).

近年來數(shù)據(jù)驅(qū)動研究快速興起，其核心思想是利用大規(guī)模的隨機(jī)建模和模型優(yōu)選，獲得無參數(shù)的客觀估算結(jié)果.將數(shù)據(jù)驅(qū)動技術(shù)引入地震活動參數(shù)計(jì)算的研究也已逐步開始.例如Kamer和Hiemer(2015)利用數(shù)據(jù)驅(qū)動技術(shù)和G-R關(guān)系估算了美國加州地區(qū)的b值空間分布，Sambridge等(2013)利用數(shù)據(jù)屬性推斷模型參數(shù)的空間分布信息，Nandan等(2017)發(fā)展了基于數(shù)據(jù)驅(qū)動的時(shí)空傳染型余震序列(ETAS)模型參數(shù)估算方法，Si和Jiang(2019)、Jiang等(2021)發(fā)展了基于數(shù)據(jù)驅(qū)動和OK1993模型的b值空間和時(shí)空分布計(jì)算方法.對比上述方法研究，時(shí)間序列b值的計(jì)算能否也引入數(shù)據(jù)驅(qū)動技術(shù)，并同時(shí)解決前述兩個(gè)瓶頸問題呢？為此我們做了方法研究，在本文中發(fā)展了“基于數(shù)據(jù)驅(qū)動的時(shí)間序列b值計(jì)算”(文中簡稱TbDD)方法，并進(jìn)行了合成地震目錄的理論測試，以及應(yīng)用于剛剛發(fā)生的2021年5月21日云南漾濞MS6.4地震序列中.本文的TbDD方法研究以及實(shí)際應(yīng)用，試圖對科學(xué)認(rèn)識2021年云南漾濞MS6.4地震序列的活動規(guī)律、震區(qū)應(yīng)力狀態(tài)，以及在其他地震序列的危險(xiǎn)性分析、工業(yè)開采誘發(fā)地震的地震災(zāi)害風(fēng)險(xiǎn)管控提供科學(xué)分析方法.

1 基于數(shù)據(jù)驅(qū)動的時(shí)間序列b值計(jì)算

本文提出的TbDD方法的構(gòu)建，包括震級-頻度分布函數(shù)的選擇、時(shí)間軸的隨機(jī)段落劃分、最佳模型選擇三個(gè)主要步驟.

步驟一：震級-頻度分布函數(shù)的選擇.為科學(xué)描述震級-頻度分布，尤其是為充分利用記錄不完整的、數(shù)量較大的小震級地震事件(Iwata, 2013; Omi et al., 2013)，目前國際上廣泛采用Ogata和Katsura(1993)給出的連續(xù)分布函數(shù)(簡稱OK1993模型).OK1993模型是在通常意義上的冪指數(shù)函數(shù)的震級-頻度分布λ0(m)基礎(chǔ)上，疊加了地震事件的檢測率函數(shù)q(m)，從而給出包括未被完整檢測的小震級事件的震級-頻度分布模型的一般表達(dá)式：

λ(m|β)=exp(-βm)q(m)，

(1)

式中的exp(-βm) 即為冪指數(shù)函數(shù)的FMD，m為地震事件的震級，冪指數(shù)β與b值之間為線性關(guān)系β=bln10，因此獲得β的估計(jì)值后即可直接折算b值.檢測率函數(shù)q(m)是取值在[0,1]范圍的正態(tài)累積分布函數(shù)，其表達(dá)式為：

(2)

式中的M是計(jì)算所用地震目錄的最大震級.參數(shù)μ是指檢測率為50%時(shí)的震級，σ表示相應(yīng)的震級范圍，兩個(gè)參量一般可結(jié)合起來共同描述地震臺網(wǎng)記錄地震的完備程度.Iwata(2013)研究認(rèn)為，由于OK1993模型是連續(xù)、漸變的非線性函數(shù)，可分別用98%置信水平下完整記錄的μ+2σ或者在99.9%置信水平下完整記錄的μ+3σ來近似描述最小完整性震級MC.

根據(jù)公式(1)和(2)可得到實(shí)際觀測地震的概率密度函數(shù)，表示為：

=e-β mq(m|μ,σ)/e(-β μ+β2σ2/2)/β

=βe-β(m-μ)-β2σ2/2q(m|μ,σ).

(3)

對于一組實(shí)際觀測地震事件{m1,m2,…,mn}，Ogata和Katsura(1993)給出了公式(3)對應(yīng)的OK1993模型的對數(shù)似然函數(shù)：

(4)

可利用最大似然法對(4)式的對數(shù)似然函數(shù)進(jìn)行擬合，獲得OK1993模型的三個(gè)參數(shù)[β,μ,σ].

步驟二：時(shí)間軸的隨機(jī)段落劃分.數(shù)據(jù)驅(qū)動技術(shù)思路實(shí)現(xiàn)的關(guān)鍵之一是構(gòu)建大量的隨機(jī)“模型”.所謂的模型是指對研究的時(shí)空間進(jìn)行隨機(jī)剖分，每一種劃分方案被稱為一個(gè)模型，這些模型將為后續(xù)的計(jì)算和最優(yōu)模型的選擇創(chuàng)造條件.在二維空間/時(shí)-空間、三維空間的隨機(jī)剖分上，常采用Voronoi網(wǎng)格剖分方法自動生成隨機(jī)網(wǎng)格，例如政治投票分區(qū)和預(yù)測(Svec et al., 2007)、物流匹配和最優(yōu)運(yùn)輸(Rubner et al., 2000)等數(shù)據(jù)驅(qū)動技術(shù)思路的實(shí)現(xiàn).由于本文的研究對象是時(shí)間序列，因此在模型構(gòu)建上需要對時(shí)間軸上的起始時(shí)間t0與終止時(shí)間t1之間的[t0,t1]時(shí)段進(jìn)行隨機(jī)的段落切分.

具體操作步驟為：(1)利用隨機(jī)數(shù)發(fā)生器在[t0,t1]時(shí)段生成隨機(jī)時(shí)間節(jié)點(diǎn)T={T1,T2,…,Tn}，n為時(shí)間節(jié)點(diǎn)的個(gè)數(shù)，由此將[t0,t1]時(shí)段切分為時(shí)段S={S1,S2,…,Sn+1}并構(gòu)成了一個(gè)模型；(2)大量重復(fù)上述隨機(jī)過程w次，獲得時(shí)間軸上的w組劃分方案{Pi,i=1,2,…,w}，也即w個(gè)模型；(3)對每一個(gè)模型(每一組劃分方案)的時(shí)段S={S1,S2,…,Sn+1}，可對落入每個(gè)時(shí)段內(nèi)的地震事件求取OK1993模型參數(shù)[β,μ,σ]，涉及的計(jì)算量為w×(n+1)次.

步驟三：最佳模型選擇.最終的參數(shù)計(jì)算是通過模型選擇來實(shí)現(xiàn)的.這是按照某種懲罰函數(shù)來約束和選擇一些最優(yōu)模型.在地震活動參數(shù)的模型選擇中，Kamer和Hiemer(2015)、Nandan等(2017)、Si和Jiang(2019)均采用貝葉斯信息準(zhǔn)則(Bayesian information criteria, 簡稱BIC)(Schwarz, 1978)來實(shí)現(xiàn)，本文也同樣采用BIC進(jìn)行最優(yōu)化模型選擇.BIC值由公式(4)的對數(shù)似然函數(shù)lnL(θ)，以及模型的自由度k、為計(jì)算使用的地震事件數(shù)目N計(jì)算獲得：

(5)

由于每個(gè)模型的每個(gè)時(shí)段的劃分是由起點(diǎn)時(shí)刻和終點(diǎn)時(shí)刻唯一決定，因此每個(gè)時(shí)段的參數(shù)由起點(diǎn)時(shí)刻、終點(diǎn)時(shí)刻、時(shí)段內(nèi)的OK1993模型參數(shù)[β,μ,σ]組成，公式(5)涉及到的每個(gè)模型的整體自由度為k=5×(n+1).此外，每一個(gè)模型的BIC值是對所有時(shí)段內(nèi)計(jì)算結(jié)果的求和Sum(BIC1,BIC2,…,BICn+1).

最佳模型是在w個(gè)模型中選擇一定比例的、BIC值最低的部分模型，例如通常選擇模型數(shù)量的5%或10%，一些研究也認(rèn)為具體的數(shù)量選擇對最終結(jié)果影響不大(Si and Jiang, 2019).最終的參數(shù)值是對這些最佳模型的參數(shù)值在時(shí)間軸上計(jì)算的集合中值，用相應(yīng)的絕對中位差(MAD)，即四分位距的一半來表示結(jié)果的不確定度.

2 方法的理論測試

2.1 合成地震目錄的生成

為驗(yàn)證TbDD方法對時(shí)間序列b值計(jì)算的有效性，需要進(jìn)行方法的理論測試，即構(gòu)建合成地震目錄并測試計(jì)算結(jié)果對輸入?yún)?shù)的還原程度.為接近真實(shí)地模擬地震活動，在合成地震目錄時(shí)設(shè)定了與后續(xù)的實(shí)際地震觀測數(shù)據(jù)計(jì)算的2021年云南漾濞MS6.4地震序列一致的時(shí)間范圍、接近的生成地震數(shù)量，也對震級-頻度分布采用OK1993模型的概率密度函數(shù).

本文在起始時(shí)間t0=2021-05-18 08∶00∶00，終止時(shí)間t1= 2021-05-26 15∶30∶00的時(shí)段內(nèi)生成合成地震目錄.按照時(shí)長t1-t0的5/12、3/12、4/12比例順次分成三個(gè)時(shí)段、設(shè)定不同的OK1993模型合成地震目錄輸入?yún)?shù)，其中令三個(gè)時(shí)段均使用μ=0.8000和σ=0.2000，而β值則分別為0.60、0.85、0.50.分別利用公式(3)的概率密度函數(shù)在每個(gè)時(shí)段生成1000個(gè)地震事件，震級范圍為0～6.4.由于生成的是有限數(shù)量的隨機(jī)地震事件，對其再次進(jìn)行OK1993模型擬合時(shí)獲得的參數(shù)往往與合成地震目錄輸入?yún)?shù)存在偏差，例如本研究中的三組參數(shù)分別為β=1.4125、μ=0.8186、σ=0.2081，β=2.0404、μ=0.7924、σ=0.2076以及β=1.1095、μ=0.8138、σ=0.2059，我們將這些新生成的合成地震目錄參數(shù)作為后續(xù)比較和評判的參考值.上述獲得的隨機(jī)地震事件的FMD和OK1993模型參數(shù)擬合值如圖1a—c所示，震級-時(shí)間分布(M-t)如圖1d所示，相應(yīng)三個(gè)時(shí)段的b0值分別為0.6134、0.8861和0.4818.

圖1 對TbDD方法測試的理論地震目錄 (a)—(c)理論地震目錄的FMD以及OK1993模型擬合結(jié)果.三個(gè)子圖分別對應(yīng)假設(shè)時(shí)段2021-05-18 08∶00∶00—05-26 15∶30∶00的前5/12、中間3/12、后4/12，均為1000個(gè)隨機(jī)生成事件.菱形為理論地震目錄的FMD，實(shí)曲線為OK1993模型擬合結(jié)果，理論地震目錄生成后的實(shí)際參數(shù)β、μ和σ分別被標(biāo)注在各子圖上；(d)理論地震目錄的震級-時(shí)間分布(M -t).圓點(diǎn)為生成的地震事件，顏色與圖(a)—(c)中相應(yīng)時(shí)段 的OK1993模型擬合曲線的顏色相一致，藍(lán)色虛線為三個(gè)時(shí)段對應(yīng)的b值、與圖(a)—(c)的β /log10相一致.Fig.1 Synthetic earthquake catalogs for testing the TbDD method (a)—(c) The magnitude-frequency distribution (FMD) of the synthetic catalogs and the fitting results of OK1993 model. The three subgraphs respectively correspond to the first 5/12, the middle 3/12, and the last 4/12 of the hypothetical period 2021-05-18 08∶00∶00—05-26 15∶30∶00. The number of earthquakes in each time segment is 1000. The diamonds show the FMD of the synthetic catalogs. The solid curves represent the fitting results of OK1993 model. The actual parameters of synthetic catalogs are marked in each subgraph in the order of β, μ, and σ; (d) The magnitude-time distribution (M -t) of the synthetic catalogs. The dots indicate generated seismic events, whose colors correspond to that of the OK1993 model′s fitting curves of the same segment in figure (a)—(c). The blue dashed line denotes the b-values of the three time segments, which are consistent with β /log10 in figure (a)—(c).

2.2 理論測試結(jié)果

對生成的t0=2021-05-18 08∶00∶00至t1=2021-

05-26 15∶30∶00時(shí)段的合成地震目錄，我們根據(jù)時(shí)間序列的復(fù)雜程度，設(shè)定TbDD方法的時(shí)間節(jié)點(diǎn)的個(gè)數(shù)n=4，也即將計(jì)算所用的時(shí)間跨度隨機(jī)切分成5段.生成m=10000個(gè)模型，其中BIC值最小的1000個(gè)被選作最佳模型.獲得的計(jì)算結(jié)果如圖2所示，圖2a和2b中的紅色曲線為b值的集合中值，灰色區(qū)域?yàn)橄鄳?yīng)b值的絕對中位差.通過與合成地震目錄的輸入?yún)?shù)b0值比較可見，由TbDD方法獲得的b值基本可以恢復(fù)b0值，尤其是當(dāng)b0值在差異較大的0.8861與0.4818之間突跳和在b0=0.6134與b0=0.4818時(shí)段內(nèi)部時(shí).這一結(jié)果展示了TbDD方法在恢復(fù)初始參數(shù)上的有效性.

為了與傳統(tǒng)的時(shí)間序列b值計(jì)算方法進(jìn)行比較，我們采用了兩種方式進(jìn)行測試：(1)采用固定地震數(shù)目的窗長和步長的固定窗口法，分別設(shè)置兩組窗長或步長為330個(gè)事件和180個(gè)事件，具體計(jì)算仍采用OK1993模型和最大似然法計(jì)算參數(shù)β并折算為b值；(2)采用180個(gè)事件的步長，但窗長不斷累積擴(kuò)大，分別沿時(shí)間軸和逆時(shí)間軸選擇地震，采用OK1993模型和最大似然法計(jì)算參數(shù)β并折算為b值.上述兩種測試結(jié)果如圖2所示.

圖2 對理論地震目錄的TbDD方法計(jì)算結(jié)果 及其與固定窗口法的比較 (a) TbDD計(jì)算結(jié)果及其與330個(gè)事件的窗長和步長的固定窗口法的比較.黑色折線為固定窗口法計(jì)算獲得的b值； (b) TbDD計(jì)算結(jié)果及其與180個(gè)事件的窗長和步長的固定窗口法的比較.綠色折線和紫色折線為采用180個(gè)事件的步長進(jìn)行累積窗長、分別沿時(shí)間軸和逆時(shí)間軸的b值.(a)和(b)中的紅色曲線為利用TbDD方法計(jì)算獲得的b值中值，計(jì)算采用5段劃分、10000次隨機(jī)剖分、1000個(gè)最佳模型，灰色區(qū)域?yàn)?相應(yīng)的b值的MAD.藍(lán)色虛折線為理論地震目錄的b0值.Fig.2 The calculation results from the TbDD method and its comparison with the fixed window method by using the synthetic earthquake catalogs (a) The results of the TbDD are compared with that of the fixed window method of the window length and step length of 330 events. The black polyline indicates the b-values calculated by the fixed window method; (b) The results of the TbDD are compared with that of the fixed window method of the window length and step length of 180 events. The green polyline and the purple polyline indicate the b-values calculated along the time axis and against the time axis respectively, and both use a step of 180 events and a continuous accumulation window to select data. The red curves in (a) and (b) show the median values of the b-values calculated by the TbDD method. Five-segment division, 10000 random partitions and 1000 best models are used in the calculation. The gray area show the MAD of the corresponding b-values. The blue dashed polylines show the b0-values of the synthetic catalogs.

對比TbDD方法和傳統(tǒng)的固定窗口法獲得的b值結(jié)果可見，(1)當(dāng)固定地震數(shù)目的窗長或步長較大時(shí)，在合成輸入?yún)?shù)穩(wěn)定的時(shí)段也可獲得相對穩(wěn)定的結(jié)果，例如圖2a的330個(gè)事件的窗長或步長.當(dāng)固定地震數(shù)目的窗長或步長較小時(shí)，b0則會由于擬合OK1993模型不確定性的提高漲落較大.這對傳統(tǒng)的固定窗口法提出了較大的約束限制；(2)在b0值出現(xiàn)突跳的位置，傳統(tǒng)的固定窗口法同樣受限于窗長的大小，無法準(zhǔn)確地恢復(fù)突跳變化.例如在差異較大的b0=0.8861與b0=0.4818之間突跳時(shí)，窗長分別為330個(gè)事件和180個(gè)事件的兩種設(shè)置均難以恢復(fù)這種突跳.甚至在采用累積窗長的順時(shí)間軸與逆時(shí)間軸的計(jì)算中，由于經(jīng)過突升的b0=0.8861時(shí)段時(shí)，累積窗長已顯著大于該時(shí)段的長度，計(jì)算獲得的b值幾乎無法反映該時(shí)段的突升現(xiàn)象.這些均揭示了傳統(tǒng)固定窗口法的明顯弊端，同時(shí)也對比展示了TbDD方法的相應(yīng)優(yōu)勢.

3 2021年云南漾濞MS6.4地震序列的實(shí)際應(yīng)用

我們將本文發(fā)展的TbDD方法用于實(shí)際的地震序列b值計(jì)算.據(jù)中國地震臺網(wǎng)測定，2021年5月21日21時(shí)48分，云南大理州漾濞縣附近發(fā)生MS6.4地震，震中位于25.67°N、99.87°E，震源深度8 km.此次地震前的5月18日、19日，震區(qū)發(fā)生兩組最大超過4級的前震活動，甚至在主震前半小時(shí)內(nèi)的21時(shí)21分震區(qū)還發(fā)生了震級高達(dá)MS5.6的前震.本文使用的數(shù)據(jù)來自中國地震臺網(wǎng)中心的地震編目系統(tǒng)《統(tǒng)一快報(bào)目錄》，研究時(shí)段為2021-05-18 08∶00∶00至2021-05-26 15∶30∶00，研究區(qū)域設(shè)定為包含此次云南漾濞MS6.4地震序列的空間范圍25.45°N—25.83°N、99.70°E—100.15°E，震級范圍M≥ML0.0.此次地震序列較為發(fā)育，在上述研究時(shí)段和研究區(qū)域內(nèi)發(fā)生MS3.0～3.9地震75次、MS4.0～4.9地震20次、MS5.0～6.9地震2次(MS5.6和MS5.2)，M≥ML0.0地震總次數(shù)達(dá)到2939次.圖3a給出了上述地震序列的震級-時(shí)間分布.

在將TbDD方法應(yīng)用于云南漾濞MS6.4地震序列中，我們采用與理論測試相同的建模策略和參數(shù)設(shè)置.即設(shè)定時(shí)間軸的隨機(jī)段落劃分為5段，隨機(jī)生成m=10000個(gè)模型，并選取其中BIC值最小的1000個(gè)最優(yōu)模型.獲得的b值的集合中值和MAD的結(jié)果如圖3b所示.為了展示云南漾濞MS6.4地震序列的最小完整性震級MC，圖3a還給出了計(jì)算獲得的μ值的集合中值，以及μ+2σ和μ+3σ的結(jié)果.由圖3b的TbDD方法計(jì)算的b值結(jié)果可見，在研究時(shí)段的云南漾濞MS6.4地震序列發(fā)生過程中b值總體低于1.0.在MS6.4主震發(fā)生前b值為0.7左右、主震前20小時(shí)還出現(xiàn)了0.1左右的下降.在主震發(fā)生后，b值起伏明顯、總體趨勢增加直至研究時(shí)段結(jié)束的0.8左右.圖中b值的起伏變化也伴隨著較大的MAD，參考圖2的合成地震目錄測試結(jié)果可以確認(rèn)，這些時(shí)段b值的時(shí)間非平穩(wěn)性，可能與震后余震區(qū)斷層上較為劇烈的應(yīng)力調(diào)整有關(guān)(Xie et al., 2019; Jiang et al., 2021).

圖3 對2021年5月21日云南漾濞MS6.4地震序列 的TbDD方法計(jì)算結(jié)果及其與固定窗口法的比較 (a) 云南漾濞MS6.4地震序列的震級-時(shí)間分布(M -t)，黑色虛線為利用TbDD方法計(jì)算獲得的參數(shù)μ值中值，藍(lán)色虛線和實(shí)線分別為最小完備震級參考的μ+2σ和μ+3σ； (b) 利用TbDD方法計(jì)算獲得的b值結(jié)果.圖中紅色曲線為b值中值，灰色區(qū)域?yàn)橄鄳?yīng)的b值的MAD，計(jì)算采用5段劃分、10000次隨機(jī)剖分、1000個(gè)最佳模型； (c) 固定窗口法計(jì)算獲得的b值結(jié)果.在采用了300個(gè)事件的步長的情況下，黑色折線給出了300個(gè)事件的滑動窗口的方式的計(jì)算結(jié)果，綠色折線和紫色折線給出了計(jì)算結(jié)果采用累積窗分別沿時(shí)間軸和逆時(shí)間軸的計(jì)算結(jié)果.圖(b)和(c)中的垂虛線為與圖(a)中M≥ML4.0事件相對應(yīng)的發(fā)震時(shí)刻.計(jì)算的時(shí)段為2021-05-18 08∶00∶00—05-26 15∶30∶00、所用事件的震級為M≥ ML0.0.Fig.3 The calculation results from the TbDD method and its comparison with the fixed window method for the Yangbi MS6.4 earthquake sequence of May 21, 2021 in Yunnan (a) The M -t of the Yangbi MS6.4 earthquake sequence in Yunnan. The black dashed lines show the median μ values calculated by the TbDD method, and the blue dashed lines and solid lines show the values of μ+2σ and μ+3σ for reference of the minimum completeness magnitude, respectively; (b) The b-values obtained by the TbDD method. The red curves and gray areas show the median and MAD of b-values, respectively. Five-segment division, 10000 random partitions and 1000 best models are used in the calculation; (c) The b-values calculated by the fixed window method. In the case of using a step length of 300 events, the black polyline shows the calculation result with sliding window of 300 events, the green polyline and the purple polyline show the calculation results using cumulative window along the time axis and against the time axis, respectively. The vertical dashed lines in figures (b) and (c) show the origin time corresponding to the events of M≥ML4.0 in figure (a). The calculated period is 2021-05-18 08∶00∶00—05-26 15∶30∶00, and the magnitude of the events is M≥ML0.0.

采用與合成地震目錄測試類似的方式，我們還給出了固定窗口法的b值計(jì)算結(jié)果，包括采用300個(gè)事件的步長和窗長，以及采用300個(gè)事件的步長、累積窗長、分別沿時(shí)間軸和逆時(shí)間軸的b值，相應(yīng)的結(jié)果如圖3c所示.將上述固定窗口法的b值計(jì)算結(jié)果與TbDD方法的計(jì)算結(jié)果比較可見，采用固定步長和窗長得到的b值的總體變化過程和數(shù)值大小，與TbDD方法的計(jì)算結(jié)果接近；而采用兩個(gè)不同方向的累積窗長計(jì)算結(jié)果，在TbDD方法得到的b值起伏變化時(shí)段存在顯著差異，這也同樣驗(yàn)證了此時(shí)段“真實(shí)”的b值應(yīng)有起伏變化.

為了進(jìn)一步驗(yàn)證TbDD方法在建模規(guī)則上是否影響b值的計(jì)算結(jié)果，我們采用多種模式計(jì)算b值：(1)隨機(jī)生成m=10000個(gè)模型、選擇1000個(gè)最佳模型，但對時(shí)間軸分別進(jìn)行3段、4段、6段的隨機(jī)劃分；(2)采用5段的時(shí)間劃分方式，但使用m=2000、m=1000和m=500個(gè)模型來計(jì)算，相應(yīng)的最佳模型數(shù)量則分別為各自的總模型數(shù)量的10%.上述多種模式的TbDD方法建模計(jì)算結(jié)果如圖4所示.將這些結(jié)果進(jìn)行比較，并與圖3b的結(jié)果比較發(fā)現(xiàn)，在隨機(jī)模型數(shù)量相同時(shí)，時(shí)間軸的隨機(jī)劃分段落越多，能夠揭示的b值隨時(shí)間的微觀起伏過程越多(圖4a)，但b值隨時(shí)間的總體趨勢一致.而在時(shí)間軸的隨機(jī)段落劃分相同時(shí)，隨機(jī)模型數(shù)量無論是10000個(gè)還是2000、1000或500個(gè)，對b值集合中值結(jié)果的影響不明顯，但MAD值將隨著隨機(jī)模型數(shù)量的增加而減小.這些多模式的計(jì)算結(jié)果間的比較，顯示了TbDD方法在計(jì)算時(shí)間序列b值時(shí)的高度穩(wěn)定性.

圖4 對2021年5月21日云南漾濞MS6.4地震序列采用 多種計(jì)算模式設(shè)置的TbDD方法計(jì)算結(jié)果的比較 (a) 采用10000次隨機(jī)剖分和不同時(shí)間劃分方式的計(jì)算結(jié)果比較； (b) 采用5段的時(shí)間劃分方式和不同隨機(jī)剖分次數(shù)的計(jì)算結(jié)果比較.圖中的實(shí)線或虛線是b值的中值，各種顏色填充的區(qū)域是與同種顏色實(shí)線或虛線相對應(yīng)的b值的MAD值.圖中實(shí)線或虛線為b值中值，各顏色區(qū)域?yàn)橄鄳?yīng)的b值的MAD.計(jì)算使用的序列時(shí)段為2021-05-18 08∶00∶00—05-26 15∶30∶00， 震級為M≥ML0.0.Fig.4 Comparison of the results of the TbDD method using multiple sets of calculation models for the MS6.4 Yangbi, Yunnan earthquake sequence (a) Comparison of calculation results by using 10000 random partitions and different time division methods; (b) Comparison of calculation results obtained from models of different random partition times using five-segments time division methods. The solid or dashed lines indicate the median of b-values, and the areas filled with various colors show the MAD of b-values corresponding to the lines with the same color. The calculated period is 2021-05-18 08∶00∶00—05-26 15∶30∶00, and the magnitude of the events is M≥ML0.0.

4 結(jié)論和討論

在廣泛用于余震區(qū)地震序列分析、工業(yè)開采區(qū)誘發(fā)地震危險(xiǎn)性分析時(shí)間序列b值的計(jì)算中，為克服傳統(tǒng)計(jì)算方法計(jì)算規(guī)則設(shè)置的人為主觀性問題，本文借鑒基于數(shù)據(jù)驅(qū)動技術(shù)思路，利用震級-頻度分布的OK1993模型、時(shí)間軸的隨機(jī)劃分和隨機(jī)模型建立、BIC的最佳模型選擇等，構(gòu)建了基于數(shù)據(jù)驅(qū)動的時(shí)間序列b值計(jì)算方法——TbDD.在經(jīng)過合成地震目錄的理論測試、2021年云南漾濞MS6.4地震序列的實(shí)際應(yīng)用分析后，獲得的主要認(rèn)識如下：

(1)利用合成地震目錄的理論測試結(jié)果表明，本文新發(fā)展的TbDD方法可較為準(zhǔn)確地恢復(fù)b值的突跳變化、穩(wěn)定地獲得各時(shí)段的參考值，能夠克服傳統(tǒng)的固定窗口法在計(jì)算規(guī)則設(shè)置上的人為主觀性等瓶頸問題.

(2)對2021年云南漾濞MS6.4地震序列的實(shí)際應(yīng)用表明，b值在MS6.4主震前為較低的0.7左右、震前20小時(shí)出現(xiàn)了0.1左右幅度的下降，表明震區(qū)在序列發(fā)生前處于較高的差應(yīng)力水平狀態(tài)；b值在MS6.4主震發(fā)生后起伏明顯、總體趨勢增加至0.8左右，可能與震區(qū)較為劇烈的應(yīng)力調(diào)整有關(guān).

(3)TbDD方法計(jì)算的時(shí)間序列b值受隨機(jī)模型數(shù)量/計(jì)算量的影響較小，具備較強(qiáng)的穩(wěn)定性，但大量的隨機(jī)模型可降低結(jié)果的不確定度；而對時(shí)間軸的隨機(jī)劃分段落設(shè)置可影響b值隨時(shí)間較為微觀的起伏變化.

本文在構(gòu)建TbDD方法中，對震級-頻度關(guān)系選用了OK1993模型，這種連續(xù)函數(shù)由于充分利用了不完整記錄的大量地震事件(Omi et al., 2013)，相比使用“一刀切”的最小完整性震級以上地震的傳統(tǒng)的基于震級-頻度G-R關(guān)系的固定窗口法，甚至一些基于數(shù)據(jù)驅(qū)動技術(shù)的空間b值計(jì)算方法(Kamer and Hiemer, 2015)，在參數(shù)估計(jì)結(jié)果的統(tǒng)計(jì)顯著性、可靠性上更具有天然優(yōu)勢.此外，最小完整性震級是隨著時(shí)間動態(tài)變化的，例如在主震發(fā)生后的短期內(nèi)，大量的余震密集發(fā)生顯著降低了地震波形的信噪比、人工檢測地震的遺漏增加(Kagan, 2004; Hainzl, 2016)，都可造成最小完整性震級的升高或檢測能力的下降.因此在基于震級-頻度G-R關(guān)系的固定窗口法使用時(shí)，可使用的“一刀切”的最小完整性震級以上地震的數(shù)量大大降低.而TbDD方法在大量的數(shù)據(jù)驅(qū)動計(jì)算中，同時(shí)估算了b值和涉及檢測能力的參數(shù)μ和σ值，使得在最小完整性震級動態(tài)演變過程中也可最大限度地利用全部數(shù)據(jù).

在應(yīng)用TbDD方法建模過程中，使用更多的隨機(jī)段落劃分可以揭示部分b值時(shí)間變化的微觀起伏.當(dāng)一個(gè)模型的隨機(jī)段落中最短的段落上樣本數(shù)過少、擬合OK1993模型誤差較大、BIC值升高，從而“拖累”這個(gè)模型難于被選為最優(yōu)模型時(shí)，客觀上也限制了無法過度細(xì)分隨機(jī)段落.我們在TbDD方法用于2021年云南漾濞MS6.4地震序列的實(shí)際案例分析中，采用的5段劃分仍具有經(jīng)驗(yàn)性，未來在此類數(shù)據(jù)驅(qū)動方法的改進(jìn)研究中，尋找客觀的真正智能化的模型構(gòu)建策略可能是需要進(jìn)一步努力的方向，包括利用Akaike貝葉斯準(zhǔn)則(ABIC)等評價(jià)方法對分段數(shù)目參數(shù)進(jìn)行優(yōu)化.此外，本文提出的TbDD方法對隨機(jī)模型次數(shù)的弱依賴性表明，大規(guī)模計(jì)算與有限次數(shù)的計(jì)算在集合中值上差別不明顯，即使是數(shù)百次的隨機(jī)模型也可得到相對可靠的結(jié)果.但大規(guī)模計(jì)算可有助于減少M(fèi)AD并確保集合中值的可靠性，這使得通過大量隨機(jī)模型進(jìn)一步取得b值計(jì)算在精度上的優(yōu)勢成為可能.

致謝本文屬于云南漾濞6.4級地震科學(xué)考察工作，地震序列研究工作組對本文給予了幫助.兩位評審專家提出了諸多建設(shè)性修改建議、對稿件質(zhì)量提升幫助很大，在此一并表示感謝.