尚園程，史保平

中國科學(xué)院大學(xué)地球與行星科學(xué)學(xué)院,北京 100049

0 引言

震群(earthquake swarm)是一類特殊的地震，通常指的是某一特定小區(qū)域在較短時間內(nèi)發(fā)生的一系列沒有明顯主震的中小地震(Scholz,2002).震群的持續(xù)時間可以是一些天，幾個月，甚至幾年.震群中會存在多個震級較大且震級相近的地震，但整個震群的能量釋放過程同常見的主震-余震序列由一個主震或幾個大震占主導(dǎo)的方式明顯不同.主震后發(fā)生的一系列余震不應(yīng)該被稱為震群，而震級特別大的多次破壞性地震，一般也不稱為震群，比如1976年松潘—平武地區(qū)發(fā)生的2次7.2級強(qiáng)震和最大震級達(dá)6.7級的余震(梁曉峰，2020).理解震群的成因機(jī)制目前仍是地震科學(xué)的重要挑戰(zhàn).

北京時間2013年10月31日11時3分，在我國吉林省松原市前郭爾羅斯蒙古族自治縣發(fā)生了MS5.5(ML5.6)地震，震中位置為44.67°N，124.12°E，11時10分，震中再次發(fā)生了MS5.0(ML5.3)地震，截至2013年11月23日，吉林前郭—乾安地區(qū)共發(fā)生MS≥5.0地震5次(吳微微等，2014)，最大震級為MS5.8(ML5.6).此次地震是一起典型的板內(nèi)震群事件(以下簡稱為前郭震群)，其發(fā)震構(gòu)造是一條與長嶺—大山子斷裂幾乎垂直的隱伏斷層(李君和王勤彩，2018).前人的研究表明前郭震群的發(fā)生與流體的運(yùn)移有關(guān)(Zhang et al.,2019；劉俊清等，2017)，但前郭震群的地震發(fā)生率相較于2008年9月—2013年9月提高了至少100倍，并且b值(古登堡-里克特定律)出現(xiàn)了明顯的降低，現(xiàn)有的研究仍不能對震群出現(xiàn)的這兩種現(xiàn)象進(jìn)行統(tǒng)一的力學(xué)機(jī)制解釋.因此研究前郭震群的成因機(jī)制對理解震群的形成機(jī)制具有重要的科學(xué)意義.

慢滑移是指以幾個小時到幾個月的時間尺度緩慢發(fā)生了有限的斷層位錯，但沒有高頻地震波輻射的斷層慢滑動事件(吳忠良和許忠淮，2013).近期的研究總結(jié)了震群的驅(qū)動機(jī)制，表明板塊內(nèi)部震群的物理驅(qū)動機(jī)制主要是慢滑移或流體運(yùn)移(王鵬，2019).而前人的觀測結(jié)果還表明慢滑移現(xiàn)象的出現(xiàn)時常伴隨著流體運(yùn)移的發(fā)生(Szeliga et al.,2004;Shelly et al.,2006;Schwartz and Rokosky,2007;Audet and Bürgmann,2014)，即慢滑移與流體運(yùn)移有關(guān).但流體運(yùn)移和慢滑移如何影響斷層的演化進(jìn)而觸發(fā)震群以及流體運(yùn)移與慢滑移相關(guān)的物理機(jī)制仍需進(jìn)一步研究.

隨著巖石物理學(xué)實(shí)驗(yàn)的發(fā)展，前人已對斷層演化的物理機(jī)制取得了深入的認(rèn)識，表明地震的演化過程就是斷層在不同的時空狀態(tài)下克服摩擦的過程.基于實(shí)驗(yàn)室內(nèi)巖石實(shí)驗(yàn)結(jié)果得出的速率-狀態(tài)摩擦定律(Rate and State Dependent Friction Law，以下簡稱RSF)已成為定量化研究斷層之間的復(fù)雜行為的一個基本的框架(Dieterich,1972,1978,1979a,b,1992,1994;Ruina,1980,1983;Rice,1983).

本文將在RSF框架下分析慢滑移和流體運(yùn)移對斷層演化的影響及其觸發(fā)震群的機(jī)制，探討流體運(yùn)移與慢滑移之間的相關(guān)性，及探討前郭震群成因的力學(xué)機(jī)制與其地震發(fā)生頻率升高與b值減小的可能的原因.

1 理論模型

1.1 速率-狀態(tài)摩擦定律

速率-狀態(tài)摩擦定律(簡稱RSF)是通過對大量巖石摩擦實(shí)驗(yàn)結(jié)果的總結(jié)而得出的，該定律如今已成為研究地震機(jī)制和斷層演化的有力工具.在RSF中，摩擦系數(shù)主要由滑移速率和狀態(tài)變量θ決定.考慮斷層上有效正應(yīng)力對狀態(tài)變量的影響，巖石摩擦實(shí)驗(yàn)給出(Linker and Dieterich，1992；Dieterich，1992)RSF可表示為：

(1)

(1)式中，τf為斷層上的摩擦強(qiáng)度；σ為斷層上的正應(yīng)力；μ為摩擦系數(shù)；A與B為摩擦參數(shù)；V為斷層的滑動速度，V*為參考速度，μ*為V=V*時的參考摩擦系數(shù)；Dc為臨界滑動位移；θ為考慮正應(yīng)力影響下的狀態(tài)變量；α為實(shí)驗(yàn)測得的無量綱常數(shù)，其取值范圍一般為0.2～0.5(Linker and Dieterich,1992).

1.2 1D彈簧-滑塊模型

本文采用1D彈簧-滑塊模型來類比真實(shí)斷層，該模型由彈性的彈簧和剛性的滑塊組成.其中，彈簧與滑塊的接觸面代表斷層面，彈簧本身則同斷層周圍的彈性介質(zhì)及斷層或裂紋自身尺度相關(guān)，而剛度系數(shù)k則與斷層面的幾何形態(tài)有關(guān).根據(jù)彈性理論給出的應(yīng)力-應(yīng)變關(guān)系，斷層內(nèi)部的剪切應(yīng)力τd的方程可表示為：

τd=τ0+τ-kδ,

(2)

(2)式中，τ0為斷層內(nèi)部的初始應(yīng)力或剩余應(yīng)力；τ為外界加載于斷層面的剪切應(yīng)力，δ為滑塊的位移；k為單位接觸面上的彈簧剛度k≈G/l(Dieterich，1992,1994)；其中，G為巖石的剪切模量；l為斷層的特征尺度.影響外部加載應(yīng)力τ的因素可以有很多，比如，遠(yuǎn)場板塊運(yùn)動、地震波擾動、地下流體擾動、人類活動等等，這些復(fù)雜因素綜合后共同決定著斷層面上的剪切應(yīng)力狀態(tài).

1.3 Dieterich 模型近似解

設(shè)定斷層在演化過程中始終保持準(zhǔn)靜態(tài)，即斷層面上的摩擦力始終等于斷層內(nèi)部的剪切應(yīng)力.當(dāng)斷層演化處于自加速階段時，(1)式中的演化方程可近似為(Dieterich,1994)：

(3)

將公式(1)中的摩擦定律和(3)式中的狀態(tài)定律代入(2)式中，我們可以得到關(guān)于滑移速率V的演化方程(Dieterich，1994；Dieterich et al.,2015)：

(4)

(4)式中，H0=B/Dc-k/σ0>0；S(t)=τ(t)-(μ0-α)σ(t)為修正的庫侖應(yīng)力方程(Dieterich et al.,2015)；τ(t)和σ(t)分別為加載于斷層上的剪切和有效正應(yīng)力；摩擦系數(shù)μ0≈0.6(Byerlee，1978)；S0為斷層內(nèi)部的初始應(yīng)力，利用公式(4)可以定量研究斷層成核過程時滑移速率V同庫侖應(yīng)力變化之間的關(guān)聯(lián)性，其中庫侖應(yīng)力變化可包含靜態(tài)或動態(tài)的應(yīng)力擾動.與公式(4)相對應(yīng)的地震發(fā)生率R可改寫為(Dieterich，1994；Segall and Lu，2015)：

(5)

1.4 流體作用下孔隙壓力對斷層滑移速率的影響

地下流體運(yùn)移是指巖漿、水、氣、油等流體在地下孔隙介質(zhì)中的流動.地下流體運(yùn)移可使斷層上的孔隙壓力升高，進(jìn)而導(dǎo)致斷層面上有效正應(yīng)力降低.為使問題不失去一般性，本文假設(shè)遠(yuǎn)場加載剪切和有效正應(yīng)力均隨時間呈線性變化，即遠(yuǎn)場應(yīng)力加載函數(shù)可表示為：

(6a)

(6b)

進(jìn)一步，我們設(shè)定斷層上的剪切應(yīng)力的增加僅由遠(yuǎn)場加載函數(shù)(6a)式所主導(dǎo)，即剪切應(yīng)力函數(shù)形式保持不變，而流體運(yùn)移作用后斷層上正應(yīng)力的擾動量為Δσ(孔隙壓力)可以是隨時間變化的函數(shù).在本文的研究中，我們僅考慮Δσ/σ0?1情形時，Δσ的變化量對斷層滑移速率的影響.孔隙壓力Δσ與流體應(yīng)力Δp之間的關(guān)系為Δp(t)=-B1Δσ(t)，其中B1為Skempton系數(shù)(Cocco and Rice，2002).因此，當(dāng)Δp>0，斷層上的有效正應(yīng)力則會出現(xiàn)明顯的下降.若孔隙應(yīng)力變化Δσ(t)可由擾動函數(shù)f(t)=Δσ0sin(ωt)表示，因而擾動量的變化由幅值Δτ0和角頻率ω所決定，而角頻率與擾動周期Tp相關(guān)(Tp=2π/ω).若流體壓力變化發(fā)生在t0

(7)

當(dāng)t0≤t≤t1時，且有Tp?Δt=t1-t0及V0≈Vp l，公式(7)也可近似為：

(8)

相應(yīng)的滑動位移為：

(9)

1.5 斷層相互作用模型和慢滑移/蠕變滑移觸發(fā)震群的機(jī)制

圖1 斷層相互作用概念模型圖中用藍(lán)色線圍成的方塊代表一個地質(zhì)塊體，地質(zhì)塊體中的黑色線段代表斷層(或裂紋)，數(shù)字i=1,2,3,…,N代表斷層的編號，Vpl為遠(yuǎn)場對塊體的加載速度，δR代表遠(yuǎn)場對塊體的加載位移.Fig.1 Conceptual model of fault interactionIn the figure,the block surrounded by blue lines represents a geological block,the black line segment in the geological block represents the fault (or crack),the number i=1,2,3,…,N represents the number of the fault,Vpl represents the loading speed and δR represents the displacement.

(10)

(11)

(12)

圖2 由慢滑或非地震滑移所觸發(fā)的地震發(fā)生率隨時間演化特征示意圖(a)，(b)和(c)由本文公式(11)與(12)所給出.(a)的取值對地震發(fā)生率的影響；的取值對地震發(fā)生率的影響；(c)t1的取值對地震發(fā)生率的影響.Fig.2 Schematic diagram of the seismicity rate of earthquake swarms triggered by slow slip or aseismic slip(a),(b)and (c)are given by formulas (11)and (12)in this paper.(a)The influence of value on seismicity rate;(b)The influence of value on seismicity rate;(c)The influence of t1 value on seismicity rate.

1.6 基于地震復(fù)發(fā)周期的近似解分析慢滑或非地震滑移對震群的影響

完整的地震復(fù)發(fā)周期一般可分為4個主要階段(Stein and Wysession,2003)：①閉鎖階段；②自加速/成核階段；③同震相；④震后松弛/滑移階段.當(dāng)震后松弛階段結(jié)束后，斷層又進(jìn)入到閉鎖階段，表明新的地震循環(huán)的開始.根據(jù)Dieterich模型(1994)，復(fù)發(fā)周期Tr可近似為：

(13)

(14)

2 2013年前郭震群序列分析

松遼盆地位于中國的東北部，地跨黑龍江、吉林、遼寧和內(nèi)蒙古，第四紀(jì)以來在盆地周邊斷裂控制下，總體活動趨勢趨于沉降.控制松遼盆地演化發(fā)展的主要地球動力是地殼深部地幔物質(zhì)的驅(qū)動和西太平洋板塊向東亞大陸邊緣的俯沖(高立新，2008；李恩澤等，2012).2013年前郭震群位于松遼盆地中部北東向扶余—肇東斷裂和西北向查干花—道字井?dāng)嗔训慕粎R部位附近(劉俊清等，2017)，地處松遼盆地中央凹陷地區(qū)，是盆地內(nèi)沉降速度最快、沉降幅度最大、坳陷最深的區(qū)域(吳微微等，2014).

本研究采用中國地震臺網(wǎng)中心提供的2008年9月—2017年7月前郭區(qū)域(44.3°N—45°N，123.6°E—124.7°E)內(nèi)的地震目錄，其M-T圖如圖4a所示.為了便于定量研究，我們把2013年前郭震群的發(fā)生時間定為2013年10月31日—2013年12月13日，把2008年9月至2013年10月30日的地震定為背景地震，把2013年12月14日—2017年7月的地震定為前郭震群結(jié)束后的“余震”.前郭震群的空間分布和時間分布如圖3和圖4b所示，ML≥5.0的地震的空間位置比較集中.前郭震群序列起始于2013年10月31日11時3分的ML5.6(MS5.6)地震，結(jié)束于2013年12月13日ML1.8地震，此后研究區(qū)10日內(nèi)沒有探測到地震發(fā)生.前郭震群序列中地方震級ML≥5.0的地震共有4個(若采用面波震級，則MS≥5.0的地震共有5個，因?yàn)榍肮鹑褐?013年11月23日6時32分發(fā)生的地震的震級為ML4.9或MS5.0)，其中最大震級地震為2013年11月23日6時4分發(fā)生的ML5.6(MS5.8)地震.前郭區(qū)域背景地震(ML≥2.0)發(fā)生頻率平均約為0.01個/d，前郭震群地震發(fā)生頻率(ML≥2.0)平均約為1.56個/d，所以前郭震群地震發(fā)生頻率相較于前郭區(qū)域背景地震發(fā)生頻率平均提高了約156倍.圖5給出了前郭震群地震發(fā)生頻率隨時間的變化關(guān)系，表明前郭震群序列之中有兩個地震發(fā)生頻率較高的時段(如圖5中箭頭所指)，一段為2013年11月4日左右達(dá)到峰值，另一段在2013年11月23日左右達(dá)到峰值，這兩個峰值點(diǎn)所在的時刻接近ML≥5.0地震發(fā)生的時刻.

圖3 前郭震群序列空間分布圖Fig.3 The spatial distribution of the Qianguo earthquake swarm

圖4 不同時間尺度下的前郭區(qū)域內(nèi)地震的M-T圖(a)前郭區(qū)域2008年9月—2017年7月地震序列M-T圖；(b)前郭震群序列M-T圖.Fig.4 M-T diagram of Qianguo earthquakes with different time spans(a)M-T diagram of earthquake sequence in Qianguo area from September 2008 to July 2017;(b)M-T diagram of Qianguo earthquake swarm.

圖5 前郭震群地震發(fā)生頻率圖由時窗推移法求得，時窗分別取5 d，7 d和9 d，每間隔2d計算一次，完備震級MC取2.0.Fig.5 Schematic of Qianguo earthquakes swarm seismicity rateThe seismicity rate is obtained by moving time window methods,the time window takes 5 days,7 days,9 days and calculated every 2 days.Minimum magnitude takes 2.0

Ross等(2020)利用美國南加州密集的臺站記錄到的連續(xù)地震波形數(shù)據(jù)和深度神經(jīng)網(wǎng)絡(luò)獲取精確地震目錄的方法，對南加州圣哈辛托(San Jacinto)斷層和埃爾西諾(Elsinore)斷裂之間的卡維拉谷(Cahuilla Valley)巖體附近2016—2019年發(fā)生的震群進(jìn)行了分析，表明該地區(qū)地震活動的遷移是因?yàn)榱黧w持續(xù)注入和遷移導(dǎo)致的，其中慢慢擴(kuò)散的前鋒面暗示這個震群主要是由流體運(yùn)移來驅(qū)動的，而不是如典型的地震序列那樣由靜態(tài)或者動態(tài)應(yīng)力觸發(fā)來驅(qū)動.這為震群是由流體運(yùn)移驅(qū)動而形成的機(jī)制提供了一個重要的依據(jù).

由靜態(tài)應(yīng)力或動態(tài)應(yīng)力觸發(fā)的地震序列類似于由主震觸發(fā)的余震序列，通常滿足大森(Omori)定律和巴特(B?th)定律.前郭震群序列的地震發(fā)生率首次開始下降后，在2013年11月17日左右又開始回升(圖4)，因此前郭震群不滿足大森定律.序列中的最大震級為MS5.8(ML5.6)，第二大震級為MS5.6(ML5.6)，面波震級差僅為0.2，且MS≥5.0的地震共有5個，因此前郭震群也不滿足巴特定律.Dieterich(1994)地震發(fā)生率模型賦予了大森定律物理意義，描述了由應(yīng)力擾動引起的地震序列的地震發(fā)生頻率逐漸衰減的過程，然而同前郭震群不滿足大森定律的原因相同，前郭震群也無法滿足該物理模型.因此，前郭震群有可能同美國南加州卡維拉谷巖體附近的震群成因機(jī)制類似，簡單的靜態(tài)或動態(tài)應(yīng)力觸發(fā)無法解釋地震的叢集現(xiàn)象.

古登堡-里克特定律(Gutenberg and Richter,1944)可表示為lgN(≥M)=a-bM，其中N(≥M)表示區(qū)域內(nèi)震級大于等于M的地震總個數(shù)，a值代表區(qū)域內(nèi)地震活動總體水平，b值的大小同區(qū)域介質(zhì)特征、應(yīng)力狀態(tài)和不均勻性等有關(guān)，反映了區(qū)域內(nèi)地震的震源特性(Rundle,1989;Marzocchi and Sandri,2003).大量資料統(tǒng)計的b值接近于1.0(陳培善等，2003；El-Isa and Eaton,2014)，針對中國東北地區(qū)的研究表明東北地區(qū)地震的b值為0.73(吳兆營等，2005).利用古登堡-里克特定律，采用最大似然法(Marzocchi and Sandri,2003)，我們分別計算了前郭區(qū)域內(nèi)背景地震、前郭震群和前郭震群的余震的a和b值，得到的結(jié)果如圖6所示，其中MC為完備震級.計算結(jié)果表明2008年9月—2017年7月期間前郭區(qū)域內(nèi)地震的b值(<0.7)普遍小于東北地區(qū)的平均b值，前郭震群的b值小于前郭區(qū)域背景地震和前郭震群的余震的b值，說明該區(qū)域內(nèi)的b值在前郭震群期間降低，而在前郭震群結(jié)束后又有上升.通過形變率與地震矩的關(guān)系，近期的研究定量地給出了b值的變化同斷層上應(yīng)力加載速率變化的近似關(guān)系(Hainzl et al.,2010；尚園程和史保平，2020)，其關(guān)系式可表示為：

圖6 不同時間尺度下的前郭區(qū)域內(nèi)地震序列震級-頻度(M-F)圖，完備震級MC取值2.0Fig.6 The Magnitude-Frequency plot for Qianguo earthquake sequences with different time spans,Minimum magnitude takes 2.0

(15)

前人的研究表明2013年前郭震群是由流體運(yùn)移的驅(qū)動而產(chǎn)生，這種流體可能與地幔上升所導(dǎo)致的深部流體向淺部的運(yùn)移相關(guān)(Zhang et al.,2019)，也可能與附近的油田開采相關(guān)(劉俊清等，2017).然而受限于觀測條件和資料，目前無法對前郭震群序列中的全部地震進(jìn)行高精度定位，從而也無法認(rèn)識前郭震群發(fā)震區(qū)域內(nèi)斷層帶的精細(xì)結(jié)構(gòu)與震群的精細(xì)動態(tài)演化過程.由于該震群的震源深度平均分布在地下10 km處(吳微微等，2014)，因此無法直接確定前郭震群震源區(qū)域及其附近是否存在斷層的緩慢滑動現(xiàn)象.本文從理論上闡述了流體的運(yùn)移可以觸發(fā)慢滑移或非地震滑移，而慢滑移或非地震滑移可導(dǎo)致臨近斷層上剪切應(yīng)力加載速率(庫侖應(yīng)力加載速率)上升，進(jìn)而導(dǎo)致b值降低.因此地下流體活動所導(dǎo)致的慢滑或非地震滑移可能為前郭震群成因的力學(xué)機(jī)制.

圖7 擬合2013年前郭震群地震發(fā)生頻率(日期：2013-08-30—2014-03-03)圓圈代表實(shí)際數(shù)據(jù)，曲線代表擬合結(jié)果.實(shí)際數(shù)據(jù)由時窗推移法求得，完備震級MC取2.0，背景地震發(fā)生頻率取0.011個/d.Fig.7 Fitting the seismicity rate of Qianguo earthquakes swarm in 2013 (date is 2013-08-30—2014-03-03)The circles represent the actual data,the curve represents the fitting result.The actual data is obtained by moving time window methods.Minimum magnitude takes 2.0.Background seismicity rate takes 0.011/d.

3 討論

前震泛指主震發(fā)生之前的較短時間內(nèi)，在主震附近發(fā)生的中小地震，加速發(fā)生的小震與古登堡-里克特定律中的b值變小是前震序列的典型特征.本文給出由流體運(yùn)移而引發(fā)的慢滑移或非地震滑移所觸發(fā)的震群序列的地震發(fā)生頻率公式(11)在早期階段也呈現(xiàn)出與前震相同的逐漸加速的特征.這是因?yàn)榧铀侔l(fā)生的前震是受主震成核所驅(qū)動的，前震斷層上的剪應(yīng)力加載速率在主震成核的作用下升高，與由慢滑或非地震滑移所觸發(fā)的震群的力學(xué)機(jī)制相同，但前震與震群的物理過程和驅(qū)動機(jī)制是完全不同的.前震序列會伴隨著破壞力較大的主震，而震群序列中大多是中小地震，若能在地震序列的初始階段判斷這段序列是前震還是震群將會十分有利于地震災(zāi)害的預(yù)測及防范，而我們對震群的研究可表明在現(xiàn)有的臺網(wǎng)精度下利用地震是否加速發(fā)生的現(xiàn)象來分辨前震序列和受慢滑移或流體運(yùn)移驅(qū)動的震群序列是不可行的.Gulia和Wiemer(2019)的研究表明，矩震級MW>6.0地震的余震的b值一般比該區(qū)域內(nèi)背景地震的b值大20%左右，而若MW>6.0地震發(fā)生后所伴隨的“余震序列”的b值相對背景地震的b值不變或變小，則該序列可能為后續(xù)更大震級地震的前震序列，并依此提出了FTLS(Foreshock Traffic Light System)以用于評估地震風(fēng)險.流體運(yùn)移所引發(fā)的慢滑或非地震滑移可以導(dǎo)致斷層上剪應(yīng)力加載速率的升高進(jìn)而觸發(fā)震群，斷層上剪應(yīng)力加載速率升高同樣會導(dǎo)致震群序列的b值降低，因此，b值降低不是前震序列所獨(dú)有的特征，利用b值變小判別前震序列是不完全準(zhǔn)確的.FTLS成立的條件是矩震級MW>6.0，而前郭震群并不滿足這一條件，未來還需針對含有MW>6.0地震的震群序列進(jìn)行研究.

前郭震群序列中MS≥5.0的地震共有5個，均呈逆沖斷層兼少量走滑性質(zhì)(吳微微等，2014)，且震源機(jī)制解整體上具有良好的一致性(李君和王勤彩，2018).針對這5個地震，表1給出了吳微微等(2014)及李君和王勤彩(2018)給出的精確定位結(jié)果，表明這5個地震的空間位置十分相近.本文的擬合結(jié)果表明前郭震群期間2級以上地震的平均復(fù)發(fā)周期降為原來的1/130，所以區(qū)域內(nèi)的一些斷層在前郭震群期間重復(fù)失穩(wěn)多次的可能性是存在的.再考慮到這5個地震具有一定的相似性，因此本文認(rèn)為目前不能排除前郭震群序列中MS≥5.0的5個地震是一個或幾個(小于等于4)5級斷層的重復(fù)失穩(wěn)的可能性.但由于數(shù)據(jù)精度的限制，目前只能從理論上提出這種可能性，其驗(yàn)證方式與方法仍是一個尚待解決的問題.

表1 前郭震群序列中MS≥5.0地震的精確定位結(jié)果(吳微微等，2014；李君和王勤彩，2018)Table 1 Accurate location results of MS≥5.0 earthquakes in Qianguo earthquakes swarm (Wu et al.,2014;Li and Wang,2018)

4 結(jié)論