章 陽(yáng) ,陳 婷 ,王利兵 ,羅 娜 ,張朋杰

（1.河北省地震局石家莊中心臺(tái),石家莊 050021;2.河北省地震局預(yù)測(cè)研究中心,石家莊 050021;3.河北省地震局紅山基準(zhǔn)臺(tái),河北 邢臺(tái) 054000）

0 引言

中國(guó)地震臺(tái)網(wǎng)中心正式測(cè)定,2019年6月17日22時(shí)55分在四川宜賓市長(zhǎng)寧縣（28.34°N,104.90°E）發(fā)生6.0級(jí)地震,震源深度8 km,震源機(jī)制解為逆沖兼走滑型地震[1]。本次地震序列較為活躍,截至8月20日,共記錄到3級(jí)以上余震64次,其中3.0～3.9級(jí)地震52次,4.0～4.9級(jí)地震8次,5.0～5.9級(jí)地震4次,分別為6月17日23時(shí)36分5.1級(jí)地震,6月18日7時(shí)34分5.3級(jí)地震,6月22日22時(shí)29分5.4級(jí)地震,7月4日10時(shí)17分5.6級(jí)地震。1900年以來(lái),震中50 km內(nèi)未發(fā)生6級(jí)以上地震;震中100 km內(nèi)發(fā)生過(guò)2次6級(jí)以上地震,震中100～200 km內(nèi)發(fā)生6.0～6.9級(jí)地震5次,這些6級(jí)以上地震均發(fā)生在川滇塊體東邊界的馬邊—大關(guān)地震帶。此次地震位于川東南地區(qū)北東向華鎣山斷裂的東南方向,發(fā)生在四川盆地邊緣的長(zhǎng)寧背斜構(gòu)造附近規(guī)模不大的次級(jí)斷層上,根據(jù)主震和余震分布方向,推測(cè)發(fā)震斷層為長(zhǎng)寧—雙河大背斜[1],這與上述6級(jí)以上地震以及發(fā)生在龍門(mén)山斷裂帶的汶川8.0級(jí)地震和蘆山7.0級(jí)地震處于不同的構(gòu)造位置。由于本次地震余震比較豐富,對(duì)于余震的危險(xiǎn)性評(píng)估很有必要,因此,本文綜合考慮同震位錯(cuò)和震后粘滯松弛效應(yīng),采用粘彈性模型計(jì)算主震的庫(kù)侖應(yīng)力變化,分析主震對(duì)余震的觸發(fā)作用,根據(jù)長(zhǎng)時(shí)間尺度庫(kù)侖應(yīng)力動(dòng)態(tài)變化特征為預(yù)測(cè)后續(xù)余震提供地點(diǎn)依據(jù)。長(zhǎng)寧地震區(qū)域地質(zhì)構(gòu)造如圖1。

圖1 長(zhǎng)寧地震區(qū)域地質(zhì)構(gòu)造圖Fig.1 Geological structure of Changning earthquake region

靜態(tài)庫(kù)侖應(yīng)力變化被應(yīng)用在分析大地震對(duì)余震分布的影響中,研究認(rèn)為靜態(tài)庫(kù)侖應(yīng)力改變量超過(guò)0.01MPa（0.1bar）可以影響余震的發(fā)生位置[2]。庫(kù)侖應(yīng)力增強(qiáng)觸發(fā)后續(xù)地震,反之則推遲后續(xù)地震發(fā)生[3],因此,庫(kù)侖應(yīng)力在余震危險(xiǎn)性評(píng)估和地震中長(zhǎng)期預(yù)報(bào)中應(yīng)用廣泛。國(guó)內(nèi)眾多學(xué)者開(kāi)展利用庫(kù)侖應(yīng)力動(dòng)態(tài)演化分析斷裂帶上地震活動(dòng)時(shí)空變化過(guò)程,進(jìn)一步分析區(qū)域地震危險(xiǎn)性。萬(wàn)永革等[4]計(jì)算1900年以來(lái)青藏高原東北部發(fā)生的7級(jí)以上地震庫(kù)侖應(yīng)力積累和演化的過(guò)程,并且討論大地震發(fā)生之間的影響,發(fā)現(xiàn)庫(kù)侖應(yīng)力對(duì)大震的觸發(fā)率達(dá)85%;邵志剛等[5]計(jì)算汶川8.0級(jí)地震引起的庫(kù)侖應(yīng)力動(dòng)態(tài)演化,結(jié)合地震學(xué)和地質(zhì)學(xué)資料分析活動(dòng)斷裂危險(xiǎn)性,并指出在分析強(qiáng)震庫(kù)侖應(yīng)力時(shí)不能忽略震后效應(yīng);徐晶等[6]研究川滇菱形塊體東邊界斷裂帶內(nèi)庫(kù)侖應(yīng)力演化及危險(xiǎn)性,認(rèn)為強(qiáng)震之間可能存在觸發(fā)作用,并根據(jù)庫(kù)侖應(yīng)力的顯著增加,給出斷裂帶上的強(qiáng)震危險(xiǎn)性區(qū);靳志同等[7]根據(jù)前續(xù)斷層在后續(xù)斷層的庫(kù)侖應(yīng)力變化給出汶川地震的破裂傳播方向;靳志同等[8]基于均勻彈性半空間模型分析九寨溝7.0地震的同震庫(kù)侖應(yīng)力變化對(duì)余震的觸發(fā);

然而,基于彈性變形的模擬[9]只能估計(jì)應(yīng)力變化的彈性響應(yīng),能模擬地震后較短時(shí)間內(nèi)的應(yīng)力應(yīng)變場(chǎng)變化特征。但對(duì)于較長(zhǎng)的時(shí)間間隔（數(shù)十年）,便不能忽略粘彈性松弛效應(yīng)。在地震應(yīng)力轉(zhuǎn)移的應(yīng)力觸發(fā)研究及地震危險(xiǎn)性估計(jì)過(guò)程中,考慮粘彈性應(yīng)力松弛十分重要。程佳等[10]利用粘彈性地殼模型計(jì)算歷史地震對(duì)九寨溝地震的同震和震后庫(kù)侖應(yīng)力作用,認(rèn)為歷史強(qiáng)震促進(jìn)九寨溝地震的發(fā)生,并計(jì)算九寨溝地震對(duì)周邊斷層的地震危險(xiǎn)性影響;尹鳳玲等[11]基于分層半無(wú)限空間粘彈性模型,計(jì)算紅河斷裂帶上25次強(qiáng)震由于同震應(yīng)力階變、震后粘滯松弛和震間構(gòu)造應(yīng)力加載綜合作用的庫(kù)侖應(yīng)力變化的演化過(guò)程,進(jìn)而分析地震危險(xiǎn)區(qū);劉琦等[12]基于位錯(cuò)理論和粘彈本構(gòu)關(guān)系,分析汶川地震和蘆山地震的同震位錯(cuò)和震后粘彈松弛效應(yīng)對(duì)周邊區(qū)域的地殼形變場(chǎng)、應(yīng)力場(chǎng)、斷層活動(dòng)狀態(tài)及地震危險(xiǎn)性等的影響;黃祿淵等[13]綜合考慮同震靜態(tài)應(yīng)力和震后粘彈性應(yīng)力分析汶川地震對(duì)九寨溝地震的影響。Pollitz等[14]認(rèn)為,用Burgers體模型能更好模擬震后的粘滯松弛效應(yīng),計(jì)算震后庫(kù)侖應(yīng)力變化。Wang等[15]編寫(xiě)粘彈性應(yīng)力應(yīng)變程序PSGRN/PSCMP程序,該程序基于Burgers粘彈性模型[16],考慮地球重力對(duì)形變場(chǎng)的影響,可計(jì)算同震和震后長(zhǎng)時(shí)間尺度的應(yīng)力變化。

因此,本文綜合考慮同震位錯(cuò)和震后粘滯松弛效應(yīng),基于Burgers粘彈性模型,利用PSGRN/PSCMP程序計(jì)算地震破裂面上產(chǎn)生的庫(kù)侖破裂應(yīng)力變化,分析主震對(duì)3級(jí)以上余震的觸發(fā)作用,進(jìn)而為后續(xù)的余震預(yù)測(cè)提供地點(diǎn)依據(jù)。同時(shí),計(jì)算粘彈性模型的長(zhǎng)時(shí)間尺度庫(kù)侖應(yīng)力動(dòng)態(tài)變化,根據(jù)粘彈性模型的震后長(zhǎng)時(shí)間尺度庫(kù)侖應(yīng)力動(dòng)態(tài)演化特征進(jìn)一步分析區(qū)域地震危險(xiǎn)性。

1 庫(kù)侖應(yīng)力變化

根據(jù)庫(kù)侖破裂準(zhǔn)則,當(dāng)巖石發(fā)生破裂時(shí),破裂面上的庫(kù)侖應(yīng)力（CFS）可定義為

式（1）中,τ為破裂面上的剪應(yīng)力;μ為內(nèi)摩擦系數(shù);σn為正應(yīng)力,張開(kāi)為正,壓縮為負(fù);p為地殼內(nèi)部的孔隙產(chǎn)生的張性應(yīng)力,S為剪切強(qiáng)度。當(dāng)剪應(yīng)力τ 越接近于S-μ(σn+p),庫(kù)侖破裂應(yīng)力CFS越接近0,材料越容易破裂。倘若μ 和S不隨時(shí)的變化而變化,那么,庫(kù)侖應(yīng)力變化可以表示為

當(dāng)ΔCFS＞0時(shí)斷層面處于加載狀態(tài),觸發(fā)后續(xù)地震,反之則推遲后續(xù)地震發(fā)生。研究中常用的庫(kù)侖應(yīng)力變化的近似表達(dá)形式

其中,μ’=μ(1-β’)為視摩擦系數(shù),β’是斯肯普頓系數(shù),取值范圍為0.7～1.0[17],研究表明,視摩擦系數(shù)不是材料的特質(zhì),主要依靠介質(zhì)的應(yīng)力變化率,不同取值只影響應(yīng)力變化的大小,不影響庫(kù)侖破裂應(yīng)力變化的空間分布[18]。在計(jì)算庫(kù)侖應(yīng)力變化過(guò)程中依據(jù)研究區(qū)域的背景資料確定視摩擦系數(shù)的大小,一般取0.2～0.8。本文采用汪榮江研究員開(kāi)發(fā)的PSGRN/PSCMP程序計(jì)算強(qiáng)震的震后庫(kù)侖應(yīng)力變化,模擬同震響應(yīng)和震后松弛效應(yīng)。

2 計(jì)算模型

2.1 地殼介質(zhì)模型

地殼介質(zhì)模型采用分層半無(wú)限空間粘彈性模型,分為完全彈性的上地殼和中地殼、Burgers粘彈性的下地殼和上地幔。本研究中的P波速度參考趙珠等[19]的研究成果,S波速度根據(jù)求得;殼幔分層結(jié)構(gòu)、地層介質(zhì)密度參考王芃等[20]的研究結(jié)果;粘滯系數(shù)參考石耀霖等[21]、李艷娥等[22]的結(jié)果,下地殼的穩(wěn)態(tài)粘滯系數(shù)取為1.0×1021Pa?s,上地幔的粘滯系數(shù)取為1.0×1022Pa?s,下地殼瞬時(shí)與穩(wěn)態(tài)粘滯系數(shù)一般相差1～1.5個(gè)數(shù)量級(jí)[23,11]。因此,下地殼的瞬時(shí)粘滯系數(shù)取為1.0×1020Pa?s,。計(jì)算庫(kù)侖應(yīng)力變化時(shí)的具體地殼介質(zhì)模型見(jiàn)表1。

表1 地殼介質(zhì)模型Table 1 Crustal medium model

圖2 余震序列的震源深度分布圖Fig.2 Focal depth distribution map of aftershock sequences

2.2 長(zhǎng)寧6.0級(jí)地震的震源參數(shù)及同震位錯(cuò)模型

長(zhǎng)寧6.0級(jí)地震斷層參數(shù)引用臺(tái)網(wǎng)中心的震源機(jī)制解：走向323°,傾角57°,滑動(dòng)角65°。同震位錯(cuò)模型參數(shù)中,破裂長(zhǎng)度和破裂寬度參考程佳等[24]經(jīng)驗(yàn)公式,計(jì)算結(jié)果約為長(zhǎng)28 km和寬9 km。走滑同震位錯(cuò)參考朱航等[25]經(jīng)驗(yàn)公式,然后根據(jù)斷層面滑動(dòng)角估算出傾滑位錯(cuò),計(jì)算得到走滑位錯(cuò)為0.638 m,傾滑位錯(cuò)為1.368 m。前人研究認(rèn)為改變摩擦系數(shù)對(duì)庫(kù)侖應(yīng)力變化影響不大,本文借鑒經(jīng)驗(yàn)取最優(yōu)值為0.4[26]。

3 計(jì)算結(jié)果

3.1 同震庫(kù)侖應(yīng)力變化及對(duì)余震的觸發(fā)作用

根據(jù)前文給出的地殼介質(zhì)模型和震源參數(shù)及同震位錯(cuò)模型計(jì)算庫(kù)侖應(yīng)力變化,接收斷層參數(shù)以該斷層周邊歷史大震的破裂參數(shù)為依據(jù),選取為走向145°,傾角90°,滑動(dòng)角0°[24]。計(jì)算庫(kù)侖應(yīng)力變化時(shí),計(jì)算深度是一個(gè)很重要的控制參數(shù),余震震源深度分布如圖2。采用基于粘彈性模型的PSGRN/PSCMP程序計(jì)算不同典型深度處的庫(kù)侖應(yīng)力變化,進(jìn)而分析6.0級(jí)主震對(duì)3級(jí)以上余震的觸發(fā)作用。所得結(jié)果見(jiàn)圖3。

圖3 不同深度處長(zhǎng)寧6.0級(jí)地震同震庫(kù)侖應(yīng)力變化及對(duì)3級(jí)以上余震的觸發(fā)作用Fig.3 The co-seismic Coulomb stress change of Changning M6.0 earthquake and its triggering effect on aftershocks with M≥3 at different depths

圖3（a）、圖3（c）和圖3（e）分別為5 km深度處、8 km深度處和10 km深度處長(zhǎng)寧地震的同震庫(kù)侖應(yīng)力變化,主震破裂面上的庫(kù)侖應(yīng)力變化分別為0.24 bar、0.9 bar和2.1 bar,均超過(guò)應(yīng)力觸發(fā)閾值0.1 bar,有利于觸發(fā)后續(xù)地震,而且隨著計(jì)算深度的增加,觸發(fā)作用有一定的增強(qiáng)。同震庫(kù)侖應(yīng)力變化均呈現(xiàn)“蝴蝶狀分布特征”,主震區(qū)的破裂方向?yàn)楸蔽飨?。?yīng)力觸發(fā)區(qū)主要集中在震區(qū)的北北西、南南東、北東東和南西西方向。圖3（b）、圖3（d）和圖3（f）分別為5 km深度處、8 km深度處和10 km深度處長(zhǎng)寧地震對(duì)3級(jí)以上余震的觸發(fā)作用,不同深度計(jì)算結(jié)果均表現(xiàn)為64個(gè)3級(jí)以上余震中,有62個(gè)發(fā)生在庫(kù)侖應(yīng)力增加的應(yīng)力觸發(fā)區(qū),2個(gè)發(fā)生在應(yīng)力影區(qū),綜合觸發(fā)率為96.9%。根據(jù)主震庫(kù)侖應(yīng)力變化特征可以較好地預(yù)測(cè)3級(jí)以上余震的發(fā)震地點(diǎn),后續(xù)余震震中位置主要集中在震區(qū)的北北西、南南東、北東東和南西西方向。

3.2 庫(kù)侖應(yīng)力動(dòng)態(tài)演化

研究長(zhǎng)寧地震長(zhǎng)時(shí)間尺度的庫(kù)侖應(yīng)力動(dòng)態(tài)演化特征,計(jì)算不同深度處Burgers粘彈性模型同震、震后50年、震后100年以及震后200年的庫(kù)侖應(yīng)力變化,接收斷層選取為走向145°,傾角90°,滑動(dòng)角0°。計(jì)算5 km深度、8 km深度和10 km深度主震破裂面上震中位置處的不同時(shí)間尺度庫(kù)侖應(yīng)力變化與同震庫(kù)侖應(yīng)力變化的差值,結(jié)果如表2所示。以10 km計(jì)算深度為例,分析不同粘滯系數(shù)的選取對(duì)計(jì)算結(jié)果的影響,結(jié)果如圖4、圖5。

表2 主震破裂面上震中位置處的不同深度不同時(shí)間尺度庫(kù)侖應(yīng)力變化與同震庫(kù)侖應(yīng)力變化差值Table 2 Difference between Coulomb stress changes and co-seismic Coulomb stress changes at different depths and different time scales on the main shock rupture surface of the epicenter location

圖4 粘彈性模型各個(gè)時(shí)間點(diǎn)庫(kù)侖應(yīng)力變化（η1=1.0x1018Pa?s、η2=1.0x1019Pa?s）Fig.4 The Coulomb stress change of viscoelastic model at various time points (η1=1.0x1018Pa?s、η2=1.0x1019Pa?s)

圖5 粘彈性模型各個(gè)時(shí)間點(diǎn)庫(kù)侖應(yīng)力變化（η1=1.0x1020Pa?s、η2=1.0x1021Pa?s）Fig.5 The Coulomb stress change of viscoelastic model at various time points (η1=1.0x1020Pa?s、η2=1.0x1021Pa?s)

由表2可以發(fā)現(xiàn),同一深度處隨著時(shí)間尺度的推移,破裂面上的庫(kù)侖應(yīng)力變化不斷增大,但幅度逐漸減小,隨著計(jì)算深度的增加,庫(kù)侖應(yīng)力變化逐漸增加,而且隨著時(shí)間尺度的推移,幅度逐漸增大。圖4和圖5為不同粘滯系數(shù)的Burgers模型計(jì)算得到的10 km深度處長(zhǎng)寧地震同震響應(yīng)、震后50年、震后100年以及震后200年的庫(kù)侖應(yīng)力變化圖,均呈現(xiàn)“蝴蝶狀”分布特征。由圖4可以發(fā)現(xiàn),當(dāng)?shù)貙诱硿禂?shù)較小時(shí),應(yīng)力傳遞速度較快,小尺度時(shí)間范圍內(nèi)庫(kù)侖應(yīng)力變化幅度較大（圖4（a）和圖（b））,隨著時(shí)間的推移,同震庫(kù)侖應(yīng)力變化將逐漸弛豫耗散,應(yīng)力觸發(fā)區(qū)和應(yīng)力影區(qū)不斷向外圍傳播,震后50年北西向出現(xiàn)新的應(yīng)力影區(qū)（圖4（b）,圖4（c）,圖4（d））,到達(dá)穩(wěn)態(tài)所需要的時(shí)間短。圖5中地層粘滯系數(shù)較大,應(yīng)力傳遞速度慢,小尺度時(shí)間范圍內(nèi)庫(kù)侖應(yīng)力變化幅度較小（圖5（a）和圖5（b））,震后100年北西向出現(xiàn)新的應(yīng)力影區(qū)（圖5（c）,圖5（d））,到達(dá)穩(wěn)態(tài)所需要的時(shí)間更長(zhǎng),因而百年尺度庫(kù)侖應(yīng)力變化幅度較小。依據(jù)長(zhǎng)寧地震發(fā)生至震后200年的庫(kù)侖應(yīng)力動(dòng)態(tài)演化（圖4和圖5）,可以看出地震發(fā)生后應(yīng)力觸發(fā)區(qū)主要集中在震區(qū)的北北西、南南東、北東東和南西西方向,其他區(qū)域?yàn)閼?yīng)力影區(qū)。將庫(kù)侖應(yīng)力動(dòng)態(tài)演化特征應(yīng)用到本區(qū)域的地震危險(xiǎn)性分析,重點(diǎn)關(guān)注上述應(yīng)力觸發(fā)區(qū),可為后續(xù)余震發(fā)生提供地點(diǎn)依據(jù)。

4 結(jié)論和討論

本文綜合考慮同震位錯(cuò)和震后粘滯松弛效應(yīng),利用PSGRN/PSCMP程序計(jì)算長(zhǎng)寧6.0地震產(chǎn)生的庫(kù)侖應(yīng)力變化,分析長(zhǎng)寧6.0級(jí)地震的庫(kù)侖應(yīng)力演化特征及其對(duì)3級(jí)以上余震的觸發(fā)作用。得出結(jié)論,3級(jí)以上余震發(fā)生在應(yīng)力觸發(fā)區(qū)的概率為96.9%,說(shuō)明3級(jí)以上余震與庫(kù)侖應(yīng)力變化有較好應(yīng)力觸發(fā)關(guān)系,可以將庫(kù)侖應(yīng)力動(dòng)態(tài)演化應(yīng)用到本區(qū)地震中長(zhǎng)期危險(xiǎn)性預(yù)測(cè)中。粘彈性模型長(zhǎng)時(shí)間尺度的庫(kù)侖應(yīng)力演化呈現(xiàn)“蝴蝶狀”分布特征,震區(qū)北北西、南南東、北東東和南西西向均位于庫(kù)侖應(yīng)力增加的應(yīng)力觸發(fā)區(qū),后期重點(diǎn)關(guān)注這些區(qū)域發(fā)生余震的可能性。

在計(jì)算庫(kù)侖應(yīng)力變化時(shí),存在不確定因素和基本假設(shè)。首先,我們假定長(zhǎng)寧6.0地震之前的累積庫(kù)侖破裂應(yīng)力為0;其次,建立地殼層狀模型時(shí),由于缺少本地區(qū)精細(xì)的三維速度模型及粘滯系數(shù)等參數(shù),本文建立的模型參數(shù)均是根據(jù)前人研究成果選取的;同時(shí)對(duì)于初始應(yīng)力場(chǎng)和長(zhǎng)期加載場(chǎng)認(rèn)識(shí)有限,對(duì)于庫(kù)侖應(yīng)力估計(jì)比較粗糙;這些都會(huì)給計(jì)算帶來(lái)一定誤差。雖不能準(zhǔn)確地預(yù)測(cè)破裂時(shí)間,但根據(jù)庫(kù)侖應(yīng)力動(dòng)態(tài)演化仍然能夠?yàn)榈卣鹬虚L(zhǎng)期預(yù)測(cè)提供重要信息。