賈若 蔣海昆 康建紅 唐春呈 張羽 宋金

1)中國地震局地球物理研究所，北京 100081

2)吉林省地震局，長春 130022

3)中國地震臺網(wǎng)中心，北京 100045

0 引言

已有研究顯示，無論是基于統(tǒng)計還是從物理角度分析，R-S模型中各參數(shù)存在著一定的相關(guān)性，即某一參數(shù)的變化會引起其他參數(shù)的變化，對模型預(yù)測結(jié)果也會產(chǎn)生一定的影響(Belardinelli et al，1999；Toda et al，2003；Gomberg et al，2005；Catalli et al，2008；Hainzl et al，2009；Cocco et al，2010；Woessner et al，2011)。但到目前為止，在詳細考慮這些相關(guān)性后應(yīng)如何針對特定序列來確定最優(yōu)擬合方案的研究尚不多見。為探討這一問題，本文以2013年吉林前郭5.8級強震群為例進行研究。根據(jù)某些特定參數(shù)的相關(guān)性，提出基于對數(shù)最大似然法的2種模型參數(shù)擬合方案：一是假定序列余震持續(xù)時間未知，即余震持續(xù)時間ta非固定條件下的擬合；二是余震持續(xù)時間ta固定條件下的擬合?；趯η肮鹑河嗾鹉夸浀恼砑敖y(tǒng)計分析，分別討論了2種方案對該震群震后早期擬合的影響。

1 基本理論

Dieterich(1994)基于巖石力學(xué)實驗結(jié)果及速率-狀態(tài)依從摩擦定律，在地震成核理論的基礎(chǔ)上，針對1個地震序列，建立了地震活動率與孕震區(qū)應(yīng)力狀態(tài)之間的定量關(guān)系

(1)

(2)

(3)

假設(shè)應(yīng)力階躍前、后的應(yīng)力加載速率恒定，根據(jù)式(1)、(2)、(3)可得一次剪應(yīng)力階躍后的地震活動率R對時間的函數(shù)

(4)

進一步研究發(fā)現(xiàn)，式(2)中等號右邊的應(yīng)力變化部分，即dτ-(τ/σ-α)dσ與庫倫應(yīng)力變化給出的形式相一致(Dieterich，2000；Toda et al，2003、2005；Catalli et al，2008；仲秋等，2013)。換言之，式(2)可進一步表示為

(5)

通過式(5)的引入，以包含有正應(yīng)力信息的庫倫應(yīng)力階躍代替了原有的剪切應(yīng)力階躍，以此作為模型的初始應(yīng)力擾動。重復(fù)上述推導(dǎo)步驟，可以獲得一次庫倫應(yīng)力階躍ΔCFS后的地震活動率R(Catalli et al，2008)

(6)

(7)

利用R-S模型，已開展了多個地區(qū)的地震活動率的模擬研究(Toda et al，1998、2003、2005；Gross，2001；Parsons et al，2000；Catalli et al，2008；仲秋等，2013；賈若等，2014)。

式(6)中的庫倫應(yīng)力變化ΔCFS可由下式計算(King et al，1994；Harris，1998；Toda et al，1998、2003、2005)

ΔCFS=Δτrake+μ′Δσn

(8)

式中，Δτrake和Δσn分別為接收斷層面滑動方向上的靜態(tài)剪切應(yīng)力變化和靜態(tài)正應(yīng)力變化(拉張為正)；μ′為視摩擦系數(shù)，包括孔隙流體和斷層面上介質(zhì)特性影響，一般取μ′為0.4～0.6(Okada，1992；King et al，1994；Toda et al，1998、2003、2005、2008)。由于基于靜態(tài)庫倫應(yīng)力觸發(fā)理論的余震活動研究較多，本文不再贅述。

利用上述模型對實際地震活動進行擬合的關(guān)鍵步驟及相關(guān)注意事項有：①確定研究區(qū)范圍。所選研究區(qū)范圍應(yīng)盡可能地排除其他叢集活動的干擾，這就要求研究區(qū)范圍不能過大，但過小的范圍又會導(dǎo)致應(yīng)力計算的不完全。通常情況下，對于同震庫倫應(yīng)力擾動而言，合適的研究區(qū)范圍應(yīng)在盡可能小的原則下，包含有全部0.01MPa以上的應(yīng)力擾動區(qū)域及全部目標地震活動。②網(wǎng)格劃研究區(qū)，計算同震庫倫應(yīng)力擾動。網(wǎng)格步長會影響計算效率、結(jié)果的精度。網(wǎng)格并非越小越好，應(yīng)保證大部分網(wǎng)格內(nèi)能夠包含一定數(shù)量的地震活動。③利用每個網(wǎng)格內(nèi)得到的應(yīng)力擾動計算地震活動率，從而得到全區(qū)域地震活動率R隨時間、空間的變化。④將理論值與實際值進行數(shù)值擬合，確定模型參數(shù)的最優(yōu)解。

2 數(shù)據(jù)準備

2.1 震例簡介及資料完備性分析

2013年10～11月，吉林松原前郭爾羅斯自治縣發(fā)生5.8級強震群(以下簡稱前郭震群)，這是松遼盆地內(nèi)部罕見的強震群型地震活動。根據(jù)防災(zāi)科技學(xué)院及松原市地震局項目研究結(jié)果(1)沈軍，2012，松原市活斷層探測與地震危險性評價。，穿過震區(qū)的主要活動斷裂有NE向的松原-肇東斷裂、NW向的查干泡-道子井斷裂以及平行于松原-肇東斷裂并位于其北部30km附近的NE向克山-大安斷裂，此次震群發(fā)生在NE向的松原-肇東斷裂與NW向的查干泡-道子井斷裂交匯地帶。根據(jù)吉林地震臺網(wǎng)記錄，該震群在1個月的時間內(nèi)先后發(fā)生5次5級以上地震，分別是10月31日MS5.5和MS5.0地震，11月22日MS5.3地震以及11月23日MS5.8和MS5.0地震。截至2016年10月24日，共記錄ML>0地震1353次，此后震區(qū)持續(xù)無震狀態(tài)，直到2018年4月12日發(fā)生ML3.6地震。據(jù)此可認為2013年前郭5.8級強震群余震活動在2016年10月24日前后已趨于結(jié)束，共持續(xù)約1089天。

圖1為前郭震群2013年10月31日以來的余震活動及附近區(qū)域(44.3°～45.0°E，123.6°～124.8°N)1968年以來的歷史地震分布，可見此次震群余震活動在空間上非常集中，基本分布在主震周圍。對比附近區(qū)域1968年以來的歷史地震活動(圖2(a))，前郭震群強震多且時間集中，伴隨后續(xù)余震頻次起伏衰減，多次主震具有多次應(yīng)力觸發(fā)特點(賈若等，2016)。

圖1 2013年吉林前郭MS5.8震群震中附近歷史地震及主要斷層分布紅色圓圈為前郭震群；藍色圓圈為1968年以來該地區(qū)歷史地震；f1：松原-肇東斷裂，f2：查干泡-道子井斷裂，f3：克山-大安斷裂

圖2 前郭震群序列及歷史地震M-t分布圖

利用最大曲率法MAXC確定前郭序列最小完整性震級MC(Wiemer et al，2000)。MAXC方法獲得的MC非常接近序列非累積頻次最高點對應(yīng)的震級值，也對應(yīng)序列開始偏離G-R關(guān)系的震級。圖3 為前郭震群ML>0地震最小完整性震級隨時間的變化(頻次統(tǒng)計步長取100)，可見平均最小完整性震級在ML1.0左右，因此下文中選擇ML1.0作為擬合震級下限，盡可能消除資料不完備對結(jié)果的影響。

圖3 基于最大曲率法的2013年前郭震群最小完備震級分析

2.2 震源機制解及同震庫倫應(yīng)力變化

基于速率-狀態(tài)依從摩擦定律的余震活動率模擬，需要計算主震對震區(qū)產(chǎn)生的應(yīng)力階躍。如前所述，在一定假設(shè)條件下，這種應(yīng)力階躍可由同震庫倫應(yīng)力變化表示。本文使用Coulomb3軟件(Toda et al，2003、2005、2008)，在彈性半無限空間中計算了前郭震群5次5級以上主震的同震庫倫應(yīng)力變化，并假定接收斷層參數(shù)與主震破裂面參數(shù)一致。

不同機構(gòu)(2)http://www.cea-igp.ac.cn/(3)https://www.usgs.gov/和研究人員(吳微微等，2014；劉俊清等，2017)分別給出前郭震群主震的震源機制解，本文按以下原則確定用于庫倫破裂應(yīng)力計算的震源機制解結(jié)果：一是盡量使用近震臺網(wǎng)結(jié)果，二是結(jié)果與其他單位相差不大，三是5次主震參數(shù)盡可能來源于同一單位。最終確定2013年10月31日MS5.5、11月MS5.3、MS5.8和MS5.0地震震源機制采用CEAIGP結(jié)果，10月31日MS5.0地震震源參數(shù)選用劉俊清等(2017)結(jié)果，5次地震用于庫倫應(yīng)力計算的震源機制解列于表1。由于震群震源機制解一致性較高，因而接收斷層統(tǒng)一采用5次主震震源參數(shù)的平均值：走向314°、傾角60°、滑動角30°。

表1 前郭震群5次5級以上地震震源機制解

依據(jù)余震精確定位(薛艷等，2015；張洪艷等，2015；劉俊清等，2017)以及野外地質(zhì)考察，認為前郭震群主發(fā)震構(gòu)造為NW向系列斷裂，且多次主震可能存在類似的破裂機制。依據(jù)斷層破裂尺度與震級的統(tǒng)計關(guān)系，計算同震破裂尺度(Wells et al，1994)

M=a1+b1·lgSRL

lgRW=a2+b2·M

(9)

其中，SRL(surface rupture length)為沿斷層面水平走向的最短破裂長度，RW(downdip rupture width)為沿斷層面深度方向的最短破裂寬度。參數(shù)ai、bi數(shù)值對不同類型斷層有不同的優(yōu)勢分布范圍(Wells et al，1994)，由于5次5級地震多以逆沖及走滑為主，而逆沖斷層的a1、b1統(tǒng)計值為4.78～5.22、1.06～1.38，a2、b2統(tǒng)計值約為-1.41～-1.81、0.38～0.44，本文選擇多組ai、bi值進行破裂尺度試算，參考余震重定位結(jié)果，認為最接近余震展布尺度的破裂尺度最為合理，最終確定a1=4.78、b1=1.08；a2=-1.61、b2=0.41。結(jié)果顯示，幾次5級地震的平均破裂尺度在10km左右。

5次5級地震的同震庫倫應(yīng)力變化計算結(jié)果如圖4 所示。計算中的視摩擦系數(shù)取為0.4，深度8km?？梢?1月23日MS5.8地震的同震庫倫應(yīng)力擾動影響范圍最大。根據(jù)所選定的接收斷層參數(shù)，在破裂面附近的一些區(qū)域內(nèi)，5次主震均產(chǎn)生了不同程度的應(yīng)力卸載，即“應(yīng)力影區(qū)”。在下文對地震活動率的擬合計算中，并未忽略這些“應(yīng)力影區(qū)”的影響，因為負的應(yīng)力階躍會引起地震活動相對背景地震活動的降低(Dieterich，1994)。

圖4 前郭震群5次主震同震庫倫應(yīng)力變化結(jié)果

3 模型參數(shù)擬合與余震活動率計算

根據(jù)式(6)可見，R-S模型中包含多個未知參數(shù)，若僅根據(jù)先驗信息確定這些參數(shù)，在實際工作中會存在很大困難；而如果假定各參數(shù)完全相互獨立，直接采用窮舉法等進行數(shù)值擬合又會導(dǎo)致大量的循環(huán)計算，效率不高且還可能會因為多參數(shù)帶來的多解性，出現(xiàn)收斂不穩(wěn)定的情況。

根據(jù)以往研究，無論是基于統(tǒng)計實驗還是物理推導(dǎo)，余震活動率模型參數(shù)間都存在一定的相關(guān)性(Hainzl et al，2009；Cocco et al，2010；Woessner et al，2011)。本文結(jié)合前人研究給出一些統(tǒng)計關(guān)系，首先通過先驗信息確定一些基本參數(shù)，進而簡化模型的可變獨立參數(shù)為(r，Aσ)，然后進一步根據(jù)r和Aσ的相關(guān)性，提出2種基于對數(shù)最大似然法的數(shù)值擬合方式，一是考慮相關(guān)性的“余震持時ta固定”條件下的擬合，二是不考慮相關(guān)性的“余震持時ta不固定”擬合。下文以前郭5.8級震群為例，簡要介紹這2種擬合方法，并在擬合效率、計算精度、結(jié)果可靠性、震例適用性等方面對比探討二者的優(yōu)劣勢。

3.1 基本參數(shù)確定及最大似然擬合

(10)

(11)

(12)

式中，ti為地震事件的發(fā)生時刻，θ表示擬合參數(shù)向量，最優(yōu)擬合的θ分布與最大lnL相對應(yīng)。在本文中，λ為地震活動率R(或余震活動率Ra)，θ則為(r，Aσ)。

需要指出的是，通常情況下，背景地震活動率r和b一樣，可通過歷史地震活動資料確定，但在吉林前郭地區(qū)，由于歷史地震活動弱、資料少，因而本文將背景地震活動率r和固有狀態(tài)參數(shù)Aσ同作為需要擬合來確定的獨立變量。同時，假定長時間尺度下的背景地震活動與序列在震級分布上具有近似特征，直接通過對余震活動的G-R關(guān)系擬合確定式(10)中b值，約為0.5。下文將進一步考慮r與Aσ的相關(guān)性，分2種情況進行擬合計算。選取用于擬合的實際觀測資料為前郭震群10月31日～12月30日之間ML>1.0的全部余震活動。

3.2 余震持續(xù)時間非固定與固定條件下的擬合結(jié)果對比

(13)

(14)

針對式(14)有2種擬合思路，一是不考慮這種相關(guān)性，直接與觀測資料進行擬合，確定Aσ與r，進而利用式(13)計算ta，稱之為“余震持時非固定”條件下的擬合；二是考慮相關(guān)性，通過實際資料預(yù)先確定ta，將Aσ用r表示，再進行擬合，稱之為“余震持時固定”條件下的擬合。2種思路擬合結(jié)果的比較如下：

(1)“余震持時ta非固定”條件下的擬合結(jié)果

此時有Aσ和r2個可變參數(shù)，基于式(12)采用窮舉法確定最優(yōu)Aσ和r。首先要確定的是2個參數(shù)的擬合范圍和搜索步長。

參考以往研究(Harris et al，1998；Toda et al，1998、2003；Guatteri et al，2001；Belardinelli et al，1999；Console et al，2006)，統(tǒng)計結(jié)果顯示，Aσ一般分布在[0.0012，0.6]，本文取搜索范圍為[0.001，0.6]，搜索步長0.001。

根據(jù)上文所述，前郭地區(qū)1978年以來的歷史地震統(tǒng)計顯示，所選研究區(qū)1978～2013年前郭震群發(fā)生前，共發(fā)生ML>1.0地震117次，則對應(yīng)的該區(qū)平均歷史地震活動率約為1.0×10-6d-1km-2。進一步考慮到歷史地震資料的不完備，實際背景地震活動率r可能大于該量級。然而，當步長較小時，跨越多個數(shù)量級的搜索范圍會使計算變得非常耗時。為了盡快找到r的精確擬合解，這里采用一種分步窮舉的方案：首先分別在4個數(shù)量級下對r進行大步長的最大似然擬合搜索，分別為1.0×10-6～9.0×10-6，1.0×10-5～9.0×10-5，1.0×10-4～9.0×10-4，1.0×10-3～9.0×10-3，每個數(shù)量級下對應(yīng)的搜索步長為分別為1.0×10-n(n=6，5，4，3)，找到對數(shù)最大似然函數(shù)值的所在范圍為1.0×10-4～9.0×10-4；然后，在該范圍內(nèi)，再取步長為0.1×10-4進行擬合。

最終，得到最大對數(shù)似然函數(shù)值為1099.65，對應(yīng)的r=1.1×10-4d-1km-2，Aσ=0.006MPa。再根據(jù)式(10)和式(13)，計算得到余震持續(xù)時間ta=435天，但這一結(jié)果與實際偏差較大。

(2)“余震持時ta固定”條件下的擬合結(jié)果

如前所述，從實際資料來看，前郭震群自2013年10月31日開始至2016年10月24日基本結(jié)束，因而可粗略判定實際余震持續(xù)時間ta=1089天，利用式(10)、(13)，此時未知參數(shù)只有r。同上文所述的對r的分步窮舉搜索方案，得到最優(yōu)結(jié)果r=1.0×10-4d-1km-2，計算得到Aσ=0.0136MPa。對應(yīng)的最大對數(shù)似然函數(shù)值1091.10。

對比上述2種條件下的擬合結(jié)果可見：

①背景地震活動r差別不大，表明模型背景地震活動率受擬合資料影響不大；

②Aσ結(jié)果差異較大，余震持時ta不固定條件下Aσ為0.006MPa，余震持時ta固定條件下Aσ為0.0136MPa，二者近乎相差1個數(shù)量級。余震持時ta不固定條件下計算得到的理論余震持續(xù)時間約為435天，與實際差異較大(實際余震活動持續(xù)時間約1089天)，因此當可以明確認定實際余震活動持續(xù)時間時，利用實際余震持續(xù)時間對模型擬合進行約束可能是有益做法；

③余震持時ta不固定及固定2種條件下得到的lnL分別為1099.65、1091.10，進一步采用AIC方法(Akaike信息準則)對2種方法進行了綜合判定：根據(jù)Akaike(1974)的研究，模型評價指標AIC=-2[lnL最大值]+2[參數(shù)個數(shù)]，AIC值越小，模型評價越好。計算得到ta不固定AIC=-2195.3，ta固定AIC=-2180.2。由此說明，對于震后早期(實際余震持時尚不能確定)而言，余震持續(xù)時間作為未知變量的擬合方式效果更好。

圖5給出2種擬合條件下理論地震活動率和實際地震活動率的對比。由圖5 可見，盡管模型評價略有差別，但2者對前郭震群實際資料的擬合差別并不大，均能夠較好地反映出實際余震頻次的總體衰減特征。另外，實際余震活動的衰減特征更為復(fù)雜，具有不規(guī)則起伏，而理論計算結(jié)果則平滑漸變。由此可知，理論模擬的頻次衰減實際上反映了擬合時段下頻次衰減的總體特征，若要反映衰減特征隨時間的變化，則可能需要更為復(fù)雜的模型構(gòu)建，且模型參數(shù)隨時間變化而變化。另一方面，從2種方法得到的前郭震區(qū)背景地震活動率r值均在約10-4d-1km-2量級分布，說明該地區(qū)具有較低的背景地震活動特征，這與實際觀測也較為一致。

圖5 2種擬合條件下對前郭震群早期余震活動率的擬合結(jié)果

3.3 余震活動率回溯性預(yù)測檢驗及討論

基于上述討論，利用模型對前郭5.8級震群全部余震活動率隨時間的變化進行了回溯性預(yù)測檢驗。由于前郭震群2016年10月24日前后就已基本結(jié)束，實際余震持續(xù)時間ta=1089天，因此采用余震持續(xù)時間ta固定的擬合方式確定參數(shù)，即Aσ=0.0136MPa、r=1×10-4d-1km-2。圖6 為截至2016年10月24日模型計算的理論余震活動率(日頻次)變化和實際觀測結(jié)果的對比。為了更好地顯示早期的余震活動變化，取時間軸為對數(shù)坐標?？梢婎A(yù)測結(jié)果與實際資料有較好的正相關(guān)，對于2組主震群前后的余震頻次起伏有較好的呈現(xiàn)。大約100天后的理論地震頻次略高于實際，說明模型計算的余震活動率衰減較慢。在11月23日的最大主震MS5.8同震期間，理論頻次遠低于實際。分析其原因，一是同震庫倫應(yīng)力變化在主破裂面附近的負值分布，這些負值區(qū)對應(yīng)于地震活動率相對背景活動率的降低區(qū)域，而在主震后早期時段，負的庫倫應(yīng)力帶來的“降低”影響最大，隨著時間推移，這種影響越來越??；二是介質(zhì)粘彈性質(zhì)(Dieterich，1994；Helmstetteret al，2009；蔣海昆等，2012)對應(yīng)力擾動的滯后影響，事實上，主震同震應(yīng)力擾動導(dǎo)致的直接余震活動并非余震活動的全部，余震發(fā)生還包括余滑、粘彈應(yīng)力松弛等觸發(fā)機制(Deng et al，1999；Perfettini et al，2004；曲均浩等，2012)；三是本文僅考慮了5次主震的余震觸發(fā)，忽略了其他較強余震所導(dǎo)致的高階余震活動(Ogata，1988；Marsan，2006)。許多強余震能夠觸發(fā)其自身的余震，特別是在震后早期強余震較為豐富階段，序列實際地震頻次包含了大量高階余震活動，而這一部分在模型中并無體現(xiàn)。另外，2014年以后，雖然也發(fā)生了幾次4級強余震，但根據(jù)實際觀測目錄顯示，相對于5.8級主震同震期間而言，這幾次4級余震后續(xù)頻次的起伏變化并不顯著，因此對模型預(yù)測結(jié)果的影響相對較弱。

圖6 基于R-S模型的前郭震群ML>1.0余震活動率預(yù)測結(jié)果

4 結(jié)論

(1)基于速率-狀態(tài)依從地震活動率模型(Dieterich，1994)，以2013年前郭5.8級地震序列為例，探討了模型參數(shù)相關(guān)性對實際擬合及預(yù)測結(jié)果的影響。根據(jù)余震持時ta與Aσ、r的相關(guān)性，可分為余震持時固定和非固定2種擬合方式。實際擬合中，如果能依據(jù)實際資料確定余震持續(xù)時間，則采用余震持時固定的擬合方式可大大提升擬合效率。但當ta不確定時，同時擬合Aσ和r反而能夠獲得更大的對數(shù)似然函數(shù)lnL值，適用于震后早期情況。由此也可以看出，由于模型參數(shù)相關(guān)性的影響，模型參數(shù)擬合存在多種可能的最優(yōu)結(jié)果，對于參數(shù)的物理解釋需倍加謹慎。

(2)模擬了前郭5.8震群余震活動率隨時間的變化。結(jié)果顯示理論預(yù)測值與實際觀測值在變化趨勢上吻合較好，特別是前郭震群的前后2組主震對余震活動的連續(xù)觸發(fā)現(xiàn)象，在5.8級地震同震期間，理論余震頻次也達到了最高值。2種條件下的擬合結(jié)果均顯示，前郭地區(qū)背景地震活動率為10-4d-1km-2量級，基本符合前郭地區(qū)歷史地震活動較弱的實際情況，也從側(cè)面證明了模型參數(shù)擬合的合理性。存在的問題是，模型以同震庫倫應(yīng)力變化作為初始應(yīng)力擾動，可能忽略了應(yīng)力松弛、余滑等實際物理過程以及高階余震活動的影響，這一假設(shè)條件可能導(dǎo)致了主震后早期理論余震頻次與實際值之間的較大偏差。

致謝：本研究使用了中國地震局地球物理研究所提供的震源機制解結(jié)果以及吉林省地震局產(chǎn)出的序列目錄資料，使用了美國Golden software公司及MathWorks公司發(fā)布的一系列繪圖及計算軟件，謹此致謝。感謝編輯部老師對文章最終修改校稿的詳細意見。