張珊珊，萬永革，王曉山，崔華偉，王曰風(fēng)，張秀萍

（1.張家口地震監(jiān)測中心站,河北 張家口 075000；2.防災(zāi)科技學(xué)院,河北 三河 065201；3.河北省地震動力學(xué)重點實驗室,河北 三河 065201；4.河北省地震局,石家莊 050021；5.山東省地震局,濟(jì)南 250000）

0 引言

唐山及其鄰區(qū)地質(zhì)構(gòu)造復(fù)雜，張家口-蓬萊地震構(gòu)造帶與華北平原地震構(gòu)造帶在這里交匯，多條斷裂遍布周圍區(qū)域（圖1）。1976年唐山地區(qū)先后發(fā)生唐山7.8級地震、灤縣7.1級地震和寧河6.9級地震，大地震雖已經(jīng)過去45年，但在這期間唐山及鄰近地區(qū)仍為小震頻發(fā)的態(tài)勢，近期又先后發(fā)生了2020年7月12日古冶5.1級地震、2021年4月16日灤州4.3級地震等中強地震?，F(xiàn)今唐山及鄰近地區(qū)發(fā)生的地震是否屬于唐山地震余震？這些中小地震受到了唐山歷史大地震怎樣的影響？仍是各專家學(xué)者研究的熱點。

圖1 唐山及鄰區(qū)斷裂分布圖

前人從多個方面研究了唐山地震與余震的對應(yīng)關(guān)系。其中，萬永革等分別用簡單破裂模型及基于分段破裂模型結(jié)合粘彈性力學(xué)模型，研究了唐山主震對灤縣地震、寧河地震及3次大震對唐山地區(qū)后續(xù)地震的庫侖破裂應(yīng)力作用情況，得到了唐山地震對灤縣地震及寧河地震有觸發(fā)作用、3次大震對后續(xù)地震有觸發(fā)作用的結(jié)論[1-2]；劉桂萍等運用3段主震破裂斷層模型，研究了唐山主震對3個余震區(qū)的彈性觸發(fā)作用，發(fā)現(xiàn)主震有助于后續(xù)地震的發(fā)生[3]；蔣長勝等利用Etas模型對2010年以來發(fā)生的3次4級以上地震進(jìn)行分析研究，得到了唐山現(xiàn)今地震為背景地震的可能性較大的結(jié)論[4]；王輝等對唐山地區(qū)2012年發(fā)生的2次M4.0以上地震進(jìn)行震源機(jī)制討論，結(jié)合對30余年唐山地區(qū)小震活動和地殼水平應(yīng)變率計算，認(rèn)為唐山強震區(qū)的余震數(shù)十年后仍在持續(xù)，震后粘性松弛可能是華北地區(qū)長時間余震持續(xù)的主要原因[5]；朱琳等基于Burgers流變模型，模擬了1976年唐山強震群引起的震后形變場以及同震和震后庫侖應(yīng)力變化，分析認(rèn)為強震群應(yīng)力調(diào)整過程已基本穩(wěn)定，古冶5.1級地震是1976年唐山大地震的余震的可能性不大[6]。

前人通過不同方法理論對唐山地區(qū)現(xiàn)今地震是否為歷史大震余震進(jìn)行分析得到的結(jié)論不盡相同，靜態(tài)庫侖破裂應(yīng)力的計算分析是用于判斷不同的地震及其余震活動分布的一個重要可靠方法。本文將在萬永革等基于粘彈性力學(xué)模型計算唐山、灤縣及寧河地震在唐山地區(qū)的庫侖破裂應(yīng)力基礎(chǔ)上[2]，采用萬永革等用形變資料計算1976年唐山地震序列的破裂分布后得出的唐山地區(qū)非均勻破裂模型[7]，對1976年唐山3次大地震發(fā)生后的庫侖破裂應(yīng)力的演變過程進(jìn)行計算，并將庫侖破裂應(yīng)力分布同地震時空特征分布對比分析。同時，收集唐山及鄰近地區(qū)2002—2020年部分中強地震震源機(jī)制解及近期發(fā)生的古冶5.1級、灤州4.3級地震的震源機(jī)制解，通過計算其受3次大地震產(chǎn)生的庫侖破裂應(yīng)力情況，進(jìn)一步探討歷史大震對現(xiàn)今地震的觸發(fā)情況。

1 庫侖破裂應(yīng)力計算

1.1 庫侖破裂應(yīng)力計算方法

根據(jù)庫侖破裂準(zhǔn)則，當(dāng)斷層面剪切應(yīng)力達(dá)到摩擦強度時，斷層面將發(fā)生剪切破壞。庫侖破裂應(yīng)力變化Δσf定義為

當(dāng)大地震發(fā)生后，短期的庫侖破裂應(yīng)力體現(xiàn)在彈性形變上，而隨著時間的積累粘彈性特征會越發(fā)明顯?，F(xiàn)有的認(rèn)知地殼的介質(zhì)為孔隙流體介質(zhì)同彈性體介質(zhì)的混合體。比較符合真實的介質(zhì)模型為麥克斯威爾體[9]，即一個彈性體和一個粘性體的串聯(lián)介質(zhì)。在短時間內(nèi)呈現(xiàn)的是彈性體的特征，在長時間尺度呈現(xiàn)的是粘性介質(zhì)特征。

麥克斯威爾體的總應(yīng)變可表示為

引入剪切模量和粘性系數(shù)的表達(dá)式為：

式中：μm、ηm分別為彈性介質(zhì)的剪切模量和粘性介質(zhì)的粘性系數(shù)，該方程為麥克斯威爾體的本構(gòu)方程[10]。

考慮到麥克斯韋粘彈性介質(zhì)模型的力學(xué)特征，平衡方程的傅立葉變換解可以用剪切模量和體變模量給出的解來表示，即

1.2 模型及參數(shù)選取

震源模型：萬永革等用形變資料計算1976年唐山地震序列的破裂分布，在計算過程中將唐山地區(qū)地震斷層分為4個大段，共計108個小斷層，各個子斷層的長度和寬度大約為5 km×5 km，最終得到了由108個子斷層組成的非均勻的破裂分布結(jié)果，將該結(jié)果作為震源模型來進(jìn)行計算[7]。

地殼模型：選取萬永革等[2]基于粘彈性力學(xué)模型計算唐山地震序列庫侖破裂應(yīng)力觸發(fā)的地殼結(jié)構(gòu)模型（表1）。其地殼模型綜合分析了劉昌銓等[12]、曾融生等[13]、徐錫偉等[14]地震探測研究結(jié)果以及于湘?zhèn)サ萚15]的反演計算結(jié)果，結(jié)合張學(xué)民等得到的體波速度結(jié)構(gòu)模型[16]得到唐山地區(qū)的綜合地殼模型；地殼密度模型采用鄭天愉等[17]運用近震記錄確定唐山地震震源過程時所給出的唐山地區(qū)地層密度；粘度模型分成了三大部分，其中上地殼以上部分（表1中1～4層）選取了一個較大的粘滯系數(shù)1.0×1030Pa·s，下地殼部分（表 1中 5～7層）和上地幔及以下（表1第8層）采用孫荀英等[18]給出的華北板塊下方深部物質(zhì)的粘性結(jié)構(gòu)，粘滯系數(shù)為7.1×1018Pa·s和 2.1×1019Pa·s。

表1 模擬中所采用的地殼結(jié)構(gòu)參數(shù)

接收斷層參數(shù)：在計算3次大地震對唐山地區(qū)的庫侖破裂應(yīng)力分布過程中需將總庫侖破裂應(yīng)力投影到可能的最優(yōu)破裂面及滑動方向上。前人對唐山地區(qū)現(xiàn)今應(yīng)力場進(jìn)行了諸多研究，本文選取了王曉山等[19]將京津冀地區(qū)分成了1°×1°的網(wǎng)格后基于標(biāo)量斷層類型值分類方法，采用MSATSI軟件進(jìn)行反演得到的結(jié)果，同時將本文研究區(qū)域分成4個區(qū)域，選取了其中對應(yīng)的4個計算結(jié)果來進(jìn)行庫侖破裂應(yīng)力計算，具體范圍及接收斷層參數(shù)結(jié)果見表2。

表2 模擬中所采用的接收斷層參數(shù)

地震數(shù)據(jù)：自1976年來唐山地區(qū)發(fā)生的地震數(shù)量較大，地震重定位可以更準(zhǔn)確地反映出地震發(fā)震位置，而雙差定位方法是一種目前應(yīng)用較為廣泛的相對定位方法，本文選取了雙差重定位后的地震分布結(jié)果同庫侖破裂應(yīng)力分布進(jìn)行對比分析。但鑒于1979年前的地震觀測資料收集較困難，故將研究區(qū)域內(nèi)1979年以來雙差重定位后的地震數(shù)據(jù)同庫侖破裂應(yīng)力分布進(jìn)行對比分析。

計算時間段劃分：張素欣等通過逐月對比分析了唐山及鄰區(qū)地震分布圖像、分布頻次及強度等特征，通過對唐山老震區(qū)的余震時空分布特征的研究，得到了1976—2015年中5個時空分布特征較為一致的時間段[20]。本文選取了其中自1979年開始的后4個時間段，并增加了2016—2021年4月這一時間段，共計5個時段來進(jìn)行庫侖破裂應(yīng)力分布情況計算。

2 計算結(jié)果

2.1 重定位后地震分布同庫侖破裂應(yīng)力分布對比結(jié)果

通過計算，1976年唐山MS7.8地震發(fā)生后的地震深度平均值為10.64 km，故選取了地下11 km的深度進(jìn)行庫侖破裂應(yīng)力分布計算?；趲靵銎屏褢?yīng)力計算的理論方法及所選取的模型參數(shù)，計算得到3次大震在指定點的總應(yīng)力變化，然后在選取的唐山地區(qū)可能的最優(yōu)破裂面上進(jìn)行投影，從而得到大地震在唐山地區(qū)庫侖破裂應(yīng)力演變過程。

研究區(qū)域內(nèi)地震雙差定位結(jié)果來自2個部分：1979—2008年定位結(jié)果由王曉山提供，該時間段共計5 906個地震重定位結(jié)果；2009—2021年重定位結(jié)果通過雙差重定位程序進(jìn)行重定位后得到。Waldhauser等在提出雙差重定位方法后，進(jìn)一步開發(fā)了雙差定位程序[21]。本文選取了分布在河北、北京、天津等地的168個地震臺站經(jīng)緯度信息，以及經(jīng)河北遙測臺網(wǎng)地震數(shù)據(jù)庫下載到的5 723個地震震相到時資料，通過雙差定位程序進(jìn)行雙差重定位后得到4 865個地震定位結(jié)果；定位均方根殘差為0.21 s，數(shù)值較小；在水平和垂直向的平均相對誤差分別為0.82 km、0.92 km、1.15 km，重定位前后震中分布及誤差分布見圖2～3。得到的重定位震中分布結(jié)果顯示，震中分布沿斷裂分布更加收斂。

圖2 重定位前后震中分布圖

圖3 水平及垂直方向重定位平均相對誤差分布

將研究區(qū)域內(nèi)對應(yīng)時間段雙差重定位后的地震分布同庫侖破裂應(yīng)力分布對比分析得到圖4。

由圖4可以看出，在區(qū)域I、II中地震分布同庫侖破裂應(yīng)力分布相關(guān)性更佳，而區(qū)域III、IV中地震分布同庫侖破裂應(yīng)力正區(qū)域整體分布在東北和西南方向有一相關(guān)性，但在細(xì)節(jié)部分對應(yīng)度不夠，尤其是大震破裂跡線周圍的地震分布十分密集但庫侖破裂正應(yīng)力分布卻無法與之完全吻合，出現(xiàn)此現(xiàn)象的原因可能是這2個區(qū)域中的斷裂分布較為復(fù)雜。由于斷層錯動的應(yīng)力釋放作用對地震的發(fā)生起到了主導(dǎo)作用，故剔除了3次大震的破裂跡線附近區(qū)域發(fā)生的地震后，分別對5個時段對比圖中正負(fù)應(yīng)力區(qū)域地震分布個數(shù)進(jìn)行統(tǒng)計（表3）。

表3 1.0級及以上地震正負(fù)應(yīng)力區(qū)域分布個數(shù)統(tǒng)計表

圖4 地震分布同庫侖破裂應(yīng)力分布對比圖

從3次大震產(chǎn)生的庫侖應(yīng)力分布情況演變過程中可以看出，4個小研究區(qū)域拼接得到的整個研究區(qū)域中庫侖破裂應(yīng)力情況較為連續(xù)?？傃芯繀^(qū)域的庫侖破裂應(yīng)力較清晰地展示出：分布大致呈正負(fù)應(yīng)力交替出現(xiàn)的蝴蝶形特征，在震后10年內(nèi)粘性因素造成的庫侖破裂應(yīng)力變化較小，在發(fā)震10年后粘性因素造成的應(yīng)力變化趨于明顯，尤其是庫侖應(yīng)力為正的區(qū)域在逐漸擴(kuò)大，這一過程與麥克斯威爾體的粘彈性介質(zhì)松弛過程一致。從地震分布演變過程可以發(fā)現(xiàn)，在大地震發(fā)生后的10年內(nèi)，地震主要集中在3次大地震的破裂跡線周圍。而此時間段內(nèi)庫侖破裂應(yīng)力的變化也較為緩慢，且地震的發(fā)生可能與斷層的應(yīng)力釋放更相關(guān)，同張素欣等對唐山地區(qū)的時空演化特征中同時間段為正常的衰減過程結(jié)論一致[20]。從1986年開始正應(yīng)力區(qū)域開始向外擴(kuò)張，其中西北和東南的正應(yīng)力花瓣變化較為明顯，1996年后西北和東南的正應(yīng)力花瓣中正應(yīng)力的變化速率更加明顯。根據(jù)前人對“庫侖破裂應(yīng)力對地震觸發(fā)作用”的研究，觸發(fā)閾值為0.01 MPa[22-24]，在庫侖破裂應(yīng)力變化的過程中處于受觸發(fā)區(qū)域地震發(fā)生的個數(shù)也明顯多于影區(qū)數(shù)目。在1986—2015年間，隨著正應(yīng)力區(qū)域的擴(kuò)大，觸發(fā)區(qū)域地震發(fā)生的個數(shù)表現(xiàn)為顯著增多特征，這說明歷史大震對唐山地區(qū)后續(xù)地震的發(fā)生具有一定的觸發(fā)作用。值得注意的是，在北西觸發(fā)區(qū)域的地震要多于東南區(qū)域，這與萬永革等唐山大地震對唐山地區(qū)余震觸發(fā)研究中給出的“大地震有利于唐山地區(qū)余震地震發(fā)生且在東南和西北花瓣區(qū)域余震數(shù)目更多”[2]的結(jié)論一致。2016—2021年4月地震分布圖結(jié)合表3的統(tǒng)計結(jié)果可以看出，在研究區(qū)域北西和東南的正應(yīng)力分布區(qū)域中仍有較多地震發(fā)生，但是發(fā)生在影區(qū)內(nèi)的地震比例較1986—2015年有了明顯的上升，該現(xiàn)象的出現(xiàn)可能預(yù)示著歷史大震對唐山地區(qū)地震發(fā)生的影響在逐漸減弱。

2.2 大地震后部分中強地震受3次大震觸發(fā)情況分析

基于較為平均的唐山地區(qū)斷層接收參數(shù)計算庫侖破裂應(yīng)力結(jié)果可能不夠精細(xì)，而使用前人計算出的單個地震的震源機(jī)制參數(shù)來計算其受庫侖破裂應(yīng)力情況則更加精確。萬永革提供了其對研究區(qū)域2002—2020年部分地震震源機(jī)制解計算結(jié)果。通過對所有震源機(jī)制解2個節(jié)面受庫侖破裂應(yīng)力情況進(jìn)行計算后，首先篩選出了2個節(jié)面受庫侖破裂應(yīng)力作用正負(fù)相同的結(jié)果；然后對受庫侖破裂應(yīng)力正負(fù)不同的兩節(jié)面，參考萬永革等[25]利用小震分布計算出的對應(yīng)區(qū)域可能的斷層面參數(shù)（表4），保留與其所得斷層面參數(shù)較為接近的節(jié)面計算結(jié)果。經(jīng)篩選后的結(jié)果見表5，受正負(fù)應(yīng)力作用地震分布見圖5。

表4 利用小震分布和區(qū)域應(yīng)力場確定唐山及鄰區(qū)的斷層面參數(shù)[25]

圖5 受正負(fù)應(yīng)力作用的地震分布圖

由圖5可以看出，收集到的震源機(jī)制解結(jié)果全部分布在研究區(qū)域的III和IV區(qū)域，主要集中在唐山斷裂、薊運河斷裂周圍。分布在III區(qū)域的地震，受庫侖破裂正應(yīng)力個數(shù)較多；分布在IV區(qū)域的地震多發(fā)生在唐山斷裂的北段，收集到的中強地震受歷史大震庫侖破裂應(yīng)力作用情況見表5。

結(jié)合圖5和表5結(jié)果分析，在III區(qū)域薊運河斷裂附近，受庫侖破裂應(yīng)力觸發(fā)的地震較多，而在IV區(qū)域特別是唐山斷裂北段受正應(yīng)力作用個數(shù)要少于位于影區(qū)的地震個數(shù)。其中在2002—2009年內(nèi)正應(yīng)力區(qū)域地震略多于應(yīng)力影區(qū)地震，在2010—2020年間應(yīng)力影區(qū)的地震遠(yuǎn)大于正應(yīng)力區(qū)域地震。造成該現(xiàn)象的原因可能是IV區(qū)域中斷裂分布較為復(fù)雜，且唐山7.8級地震震中附近因多次大震余震發(fā)生導(dǎo)致該區(qū)域地下介質(zhì)較為破碎，不易積累大的應(yīng)力，同時破碎會發(fā)生蠕滑消耗能量，所以斷層活動的作用力向NE方向傳遞，進(jìn)一步引起了該區(qū)域新的活動[20]，故該區(qū)域受到歷史大震庫侖破裂應(yīng)力的調(diào)制作用較小。

表5 部分中強地震受歷史大震庫侖破裂應(yīng)力作用情況

2.3 古冶5.1、灤州4.3級地震受庫侖破裂應(yīng)力影響分析

2020年7月12日發(fā)生了古冶5.1級地震，2021年4月16日發(fā)生了灤州4.3級地震。本文選取萬永革等綜合各家結(jié)果得到的震源機(jī)制中心解作為接收斷層參數(shù)，即：古冶地震（走向：238.96°，傾向：76.76°，滑動角：-172.58°）、灤州地震（走向：295°，傾向：75°，滑動角：-22°），進(jìn)行其所受庫侖破裂應(yīng)力計算，得到結(jié)果為古冶地震受歷史大震庫侖破裂的應(yīng)力作用，約為-0.05 MPa（圖6），灤州地震受歷史大震庫侖破裂應(yīng)力作用約為-2.0 MPa（圖7），計算結(jié)果顯示這2次地震受到了抑制作用。

圖6 古冶地震受大震庫侖破裂應(yīng)力圖

圖7 灤州地震受大震庫侖破裂應(yīng)力圖

3 討論和結(jié)論

3.1 討論

破裂模型的選取是計算庫侖破裂應(yīng)力中的一個較為關(guān)鍵的環(huán)節(jié)。本文在庫侖破裂應(yīng)力計算中選取的破裂模型將唐山地區(qū)的破裂結(jié)構(gòu)分為了5 km×5 km的共計108個子破裂分布，較前人的破裂模型精細(xì)了許多。萬永革等2008年對唐山大震庫侖破裂應(yīng)力的計算中選取的破裂模型較前人的計算引入了10 km范圍內(nèi)外的修正[2]，已將結(jié)果精確性進(jìn)行了提升。而本文中選取的模型子破裂結(jié)構(gòu)的范圍寬度為5 km，遠(yuǎn)小于其對斷層附近的修正范圍，故本文中的計算結(jié)果可能更接近真實情況。

基于較為精細(xì)的破裂模型計算出的庫侖破裂應(yīng)力分布結(jié)果同重定位后的地震分布對比，理論上可以更加清晰地展示出地震分布同庫侖破裂應(yīng)力分布的對應(yīng)關(guān)系。但就所得結(jié)果可以看出，地震分布同庫侖破裂正應(yīng)力分布的大方向上有一定的對應(yīng)關(guān)系，而在應(yīng)力影區(qū)的地震數(shù)目卻仍然較大，尤其是地震在歷史大震的破裂跡線周圍更加集中。在該地震集中區(qū)域內(nèi)，張素欣等[20，26]研究發(fā)現(xiàn)在1996年之前余震分布較均勻，1996年以來余震在唐山-古冶斷裂的東北端點附近較集中，但是庫侖破裂應(yīng)力的變化與這些特征對應(yīng)關(guān)系不夠明顯。這可能是唐山斷裂附近區(qū)域斷裂分布多且復(fù)雜，斷裂間相互的吸納能力和阻隔作用及地殼介質(zhì)較破碎發(fā)生蠕滑對能量積累消耗產(chǎn)生影響[20]造成的。

通過進(jìn)一步對2002—2020年部分中強地震受庫侖破裂應(yīng)力作用計算結(jié)果分析，得到了歷史大震對唐山斷裂及薊運河斷裂附近區(qū)域的觸發(fā)情況。由于部分震源機(jī)制解結(jié)果因2個節(jié)面受力的正負(fù)情況不一致或與計算所得斷層面結(jié)果不相近而進(jìn)行剔除，因此用于計算的震源機(jī)制解數(shù)量較少，所得結(jié)果有誤差存在。就該區(qū)域庫侖破裂應(yīng)力作用情況，后續(xù)需綜合更精確的斷層接收參數(shù)、斷裂間的相互作用、地介質(zhì)能量變化等因素進(jìn)行計算。

3.2 結(jié)論

1）地震分布同庫侖破裂應(yīng)力分布對比分析顯示，在正負(fù)應(yīng)力花瓣分布區(qū)域中正應(yīng)力區(qū)域地震分布的地震較為密集且數(shù)量較多，地震分布與正應(yīng)力花瓣形分布方向有一定的相關(guān)性，在北西和東南方向觸發(fā)作用較為明顯。從時間上分析，歷史大震發(fā)生10年后對唐山及鄰近地區(qū)地震發(fā)生的觸發(fā)作用更加凸顯。近年來發(fā)生在庫侖破裂應(yīng)力影區(qū)的地震占比顯著增大，可能預(yù)示著歷史大震對唐山地區(qū)地震發(fā)生的影響在逐漸減弱。

2）通過對2002—2020年部分地震震源機(jī)制結(jié)果進(jìn)行受庫侖破裂應(yīng)力作用計算，認(rèn)為歷史大震可能對唐山斷裂北段地區(qū)地震發(fā)生的觸發(fā)作用較小，對薊運河斷裂附近區(qū)域地震發(fā)生的觸發(fā)作用較大。

3）通過對2020年7月12日古冶5.2級地震、2021年4月16日灤州4.3級地震受歷史大震庫侖破裂應(yīng)力作用計算分析，認(rèn)為這2次地震受到了歷史大震的抑制作用。

致謝審稿專家提出了寶貴意見，本文使用了Waldhauser F.和Ellsworth L提供的雙差定位程序，Wessel P.和W.H.F.Smith提供的GMT4繪圖程序，在此一并致謝。