1998年以來伽師地震(MW≥6.0)應力相互作用與強震活動的關系

周 云 潘正洋 王衛(wèi)民 何建坤 王 洵 李國輝

1)中國地震局地球物理研究所，北京 100081 2)中國科學院青藏高原研究所，大陸碰撞與高原隆升重點實驗室，北京 100101 3)中國地震局地震預測研究所，北京 100036

0 引言

據中國地震臺網中心正式測定，2020年1月19日在新疆喀什地區(qū)伽師縣發(fā)生MS6.4地震(下文簡稱2020年地震)，震源深度為16km。根據USGS地震目錄可知，本次地震是繼1996—1998年一系列地震以及2003年MW6.3地震(下文簡稱2003年地震)以來，發(fā)生在伽師地區(qū)的又一次MW>6的地震事件。基于InSAR資料得到的發(fā)震斷層模型，結合破裂運動學特征推測，這次伽師地震的發(fā)震構造是柯坪塔格褶皺帶山前出露的柯坪塔格逆斷裂(李成龍等，2020)。伽師地區(qū)強震活動頻繁，與所處的大地構造環(huán)境密切相關。該區(qū)位于天山、帕米爾和塔里木塊體3大構造的接觸部位，新生代以來，受印度-歐亞板塊遠程碰撞作用的影響，天山造山帶重新活動并再次隆升，成為歐亞大陸內部規(guī)模最大的再生造山帶(Tapponnieretal.，1977；Avouacetal.，1993)。同時，天山造山帶夾持在古老而穩(wěn)定的準噶爾地塊和塔里木地塊之間，在持續(xù)隆升的同時向兩側擴展，向S逆沖于塔里木地塊之上，在山前和盆地北緣形成了逆沖斷層和大規(guī)模褶皺帶(張培震等，2003)，這導致伽師地區(qū)成為世界范圍內破壞性地震頻發(fā)、致災程度嚴重的強震活動區(qū)。

地震應力觸發(fā)理論認為，當一個地震發(fā)生后，它將調節(jié)周邊活動斷裂的應力狀態(tài)，進而觸發(fā)或抑制周邊潛在斷裂的地震危險性(Dengetal.，1997；Linetal.，2004)。這種調節(jié)作用包括彈性同震應力階變和震后黏彈性松弛導致的應力場變化。在早期研究中，人們普遍采用彈性半空間位錯模型研究大地震與后續(xù)地震之間的應力觸發(fā)關系。如1992年美國加州MW7.3 Landers地震，King等(1994)計算發(fā)現地震使距震中30km處Big Bear地區(qū)的同震庫侖應力值增加了0.3MPa，進而觸發(fā)了隨后的MW6.3地震；針對1857年Fort Tejon大地震的研究發(fā)現，1857—1907年間的13次南加州5.5級以上的地震中，至少有11次與1857年的地震有應力觸發(fā)關系(Harrisetal.，1996)。然而，僅用同震庫侖應力難以解釋所有的地震觸發(fā)作用，這是由于巖石圈的下地殼和上地幔并非是完全彈性而是黏彈性的，在震后較長的時間里(幾十a甚至幾百a)黏滯松弛效應可將應力傳遞到孕震層，從而改變斷層的應力狀態(tài)(Pollitz，1992；Lorenzo-Martínetal.，2006)。如對1992年MW7.3 Landers地震和1999年MS7.1 Hector Mine地震的研究表明，如果僅考慮同震應力的作用，Hector Mine地震位于Landers地震的應力影響區(qū)(Parsonsetal.，2000)；但如果同時考慮同震應力和震后黏彈性效應，則可認為Landers地震觸發(fā)了后續(xù)地震(Freedetal.，2001)。基于黏彈性層狀介質模型，沈正康等(2003)模擬了1937年以來發(fā)生在東昆侖斷裂帶上的5個M>7地震斷層錯動導致的應力演化過程，結果表明前面4個地震對2001年可可西里MW7.8地震有觸發(fā)作用。對紅河斷裂帶庫侖應力變化的研究發(fā)現，1833年發(fā)生的嵩明M8.0地震至今仍然影響著紅河斷裂帶的地震活動(尹鳳玲等，2017)。這些研究結果均表明了考慮黏彈性應力松弛過程在估計地震危險性及分析地震應力觸發(fā)時的重要性。因此，在計算地震導致的長期應力場變化時，同震應力和震后黏彈性效應都不可忽視。

1998年以來，伽師地區(qū)發(fā)生了3次MW>6的地震，這些地震相互之間的距離均<50km(圖1)。那么，這些地震之間是否也有著相互觸發(fā)關系？如果有，同震應力階變和震后黏彈性松弛，哪個發(fā)揮了主要作用？為了解決這一問題，本文首先利用地震波形資料，反演了前2次地震的震源機制與破裂過程；然后基于更符合實際的分層黏彈性地球介質模型和得到的破裂滑動分布結果，采用已有的研究方法(Wangetal.，2006)計算了1998年和2003年2次地震對伽師地區(qū)造成的應力場變化及其對2020年地震造成的影響，以研究三者的相互觸發(fā)作用。

1 地震震源特征

震源機制解顯示，2003年地震的震源深度為15.2km，節(jié)面1和節(jié)面2的走向、傾角和滑動角分別為78°、73°、56°和293°、20°、124°。考慮到該地震有一定的逆沖分量，而節(jié)面1的傾角比較大，同時地震考察結果表明發(fā)震斷層是一條近EW走向且向N傾的逆斷層(沈軍等，2006)，因此選擇節(jié)面2為斷層破裂面。將震源機制反演計算的理論地震圖與觀測波形資料進行擬合，結果顯示波形的相關性較好(圖2)。表1 為2次地震的震源機制解。

圖2 1998年和2003年2次伽師地震的震源機制解Fig.2 Focal mechanisms of two Jiashi earthquakes in 1998 and 2003.采用下半球投影，同時給出了點源模型的P波垂向位移理論圖(紅線)與資料(黑線)的擬合情況。圖形下方給出了2組節(jié)面解(圖a中，λ、δ、θ、h分別表示斷層滑動角、斷層傾角、斷層走向、震源深度)和點源模型的震源時間函數(b)

表 1 1998年和2003年2次MW>6地震震源及斷層參數Table 1 Source and fault parameters of the two earthquakes(MW>6)in 1998 and 2003

根據有限斷層反演方法(王衛(wèi)民等，2008；周云等，2015)，結合2次地震的點源模型反演結果，選用沿走向長51km、沿傾向寬18km的斷層面，并將其劃分成3km×2km的153個子斷層，使用P波和SH波數據進行斷層面時空破裂過程反演。計算同震位移響應采用基于分層介質模型的廣義反射透射系數矩陣方法(謝小碧等，1989)，震源區(qū)的地殼結構模型獲取自Crust2.0。反演使用非線性最優(yōu)化的模擬退火技術，目標函數由采用波形相關系數與同震位移矢量差的歸一化加權疊加定義。圖3 為反演得到的2次伽師地震的斷層滑動分布。結果顯示，2次地震均為MW6.3地震，1998年地震的破裂較為集中，最大滑動量為60cm，標量地震矩為3.22×1018N·m，該結果與趙翠萍(2006)的滑動分布結果類似： 最大滑動量為63cm，呈雙側破裂特征。2003年地震的破裂較為分散，最大滑動量為56cm，標量地震矩為3.15×1018N·m。其中2003年的地震以SE向單側破裂為主，這與沈軍等(2006)野外考察得出的該地震為N傾逆沖斷層自NW向SE由深至淺破裂的結果一致。

圖3 1998年和2003年2次伽師地震斷層破裂滑動分布Fig.3 Fault slip distribution of Jiashi earthquakes in 1998 and 2003.a 1998年8月27日MW6.3地震斷層滑動分布；b 2003年2月24日MW6.3地震斷層滑動分布

2 應力計算方法及參數

2.1 研究方法

地震發(fā)生后，同震和震后應力調整均會導致斷層面上庫侖應力變化。根據庫侖破裂準則，接收斷層的靜態(tài)庫侖應力變化可以表示為

ΔCFS=Δτ+μΔσ

(1)

式中，ΔCFS為斷層面上庫侖應力的變化量，Δτ為斷層面上剪切應力的變化量，μ為斷層的有效摩擦系數，Δσ為斷層面上正應力的變化量，以拉伸為正。在不同的研究中，μ的取值也有所不同，但不會改變庫侖應力變化的整體分布形態(tài)(張國宏等，2008；萬永革等，2009)。在本文的計算中，取μ=0.4。因本次研究主要計算2003年和2020年地震發(fā)震斷層面受到的影響，參考前文反演結果和王衛(wèi)民等(1)http: ∥www.itpcas.ac.cnkycgyjcg/202001/t20200122_5494135.html。的結果，選取地下15km為計算庫侖應力變化的深度。

在計算庫侖應力時，本文使用了德國波茨坦地球科學中心(GFZ)Wang等(2006)開發(fā)的PSGRN/PSCMP程序。該軟件主要包括2個部分，PSGRN用來計算給定分層模型下基本位錯源的時變格林函數，所得結果為第2個程序PSCMP的基礎數據，之后通過線性疊加方法對點位錯進行計算，得到同震和震后形變的結果。計算時重力和分層介質影響都被考慮在內，通過使用基于反混淆技術的快速傅里葉變換得到格林函數，保證了結果的穩(wěn)定性。

2.2 模型參數

本文的模型計算范圍為(38°～41°N，76°～79°E)。在進行格林函數計算時，需要先建立地下速度介質模型。選擇何玉梅等(2001)以及趙翠萍(2006)研究伽師地震群采用的速度模型(表2)，該模型綜合了上地殼速度反演結果和劉啟元等(2000)利用接收函數得到的深部速度結構。石耀霖等(2008)認為，塔里木塊體中地殼的等效黏滯系數為1021～1023Pa·s，下地殼的黏滯系數為1021～1022Pa·s。本模型中彈性層厚度為16km，黏彈性中下地殼采用Burgers體模型。中地殼的穩(wěn)態(tài)黏滯系數取1.0×1022Pa·s，下地殼穩(wěn)態(tài)黏滯系數取2.0×1021Pa·s，瞬時與穩(wěn)態(tài)黏滯系數一般差1個數量級(Ryderetal.，2011)。地幔黏滯系數比下地殼大得多，取為2.0×1022Pa·s。

表 2 伽師地區(qū)巖石圈地層介質模型Table 2 Layered lithosphere medium model in Jiashi region

2003年伽師MW6.3地震的發(fā)震斷層(圖1 中的F1)為隱伏斷裂，未能很好地確定其斷層參數，因此在設置接收斷層參數時，我們參考了該地震的震源機制解。2020年MW6.0地震的發(fā)震斷層(圖1 中的F2)為柯坪塔格逆斷裂的一部分，我們根據周德敏(2013)的結果給出了斷層滑動參數，如表3 所示。

表 3 2003年和2020年地震的發(fā)震斷層參數Table 3 Fault parameters of earthquakes in 2003 and 2020

3 計算結果

3.1 對發(fā)震斷層應力的影響

圖4 發(fā)震斷層上的庫侖應力變化Fig.4 Coulomb stress change on the seismogenic faults.a 1998年地震的同震應力；b 2003年地震發(fā)生之前的應力分布；c 2003年地震發(fā)生之后的應力分布，黑色小圓圈為2003年2月24日—2003年4月30日M>3.5的余震(黃媛等，2006)；d 2020年地震發(fā)生之前的應力分布

為了進一步研究1998年以來的3次MW≥6地震之間的相互關系，我們計算了同震及黏彈性松弛導致的2003年2月24日MW6.3地震的發(fā)震斷層F1和2020年1月19日MW6.0地震的發(fā)震斷層F2上的庫侖應力變化，如圖4 所示。圖4a為1998年地震導致的發(fā)震斷層上的同震庫侖應力變化，可以看到本次地震使斷層F1西段的應力值增加，東段的應力值減小，變化量基本都超過了0.01MPa；同時，該地震使F2西段應力減小，東段應力增加。經過4.5a的調整后，圖4b顯示應力分布和應力值并沒有明顯的變化，2003年地震位于應力增加區(qū)范圍內，應力增加值為0.01～1MPa，達到了0.04MPa。圖4c為在原有的應力基礎上增加了2003年地震同震應力的情況，本次地震明顯調整了斷層F1的應力狀態(tài)，卸載了F1西段的應力，增加了東段的應力，同時也擴張了斷層F2東段的應力增強區(qū)。黃媛等(2006)重定位后的余震主要集中在斷層的應力增強區(qū)，同時余震的分布也支持發(fā)震斷層N傾的結論。在2020年MW6.0地震發(fā)生之前，經過近20a的震后調整，斷層面上應力變化不大，2020年地震位于應力增強區(qū)內，應力增加值介于0.003～0.006MPa之間，為0.005MPa。

3.2 改變黏滯系數對斷層庫侖應力分布的影響

圖5 調整黏滯系數后發(fā)震斷層上的庫侖應力變化Fig.5 Coulomb stress change on the seismogenic faults after adjusting viscosity coefficient.a 1998年地震的同震應力；b 2003年地震發(fā)生之前的應力分布；c 2003年地震發(fā)生之后的應力分布，黑色小圓圈為2003年2月24日—4月30日M>3.5的余震(黃媛等，2006)；d 2020年地震發(fā)生之前的應力分布

前文的計算結果顯示，震后調整對應力分布的作用并不明顯。在計算震后的應力場演化時，地層黏滯系數的選取對結果可能會有一定影響。我們主要參考了石耀霖等(2008)給出的塔里木盆地的黏滯系數參數，由于發(fā)震區(qū)位于盆山交會處，而山體與盆地的黏滯系數一般并不相等，因此，本文另外選取了一組黏滯性參數進行對比。由于研究區(qū)域西南天山地殼和地幔黏滯系數尚無較好的約束，故參考了青藏高原的參數數值。石耀霖等(2008)認為青藏高原下地殼等效黏滯系數較低，約為1019～1020Pa·s，因此將中地殼和下地殼的黏滯系數調低為5.0×1020Pa·s 和5.0×1019Pa·s，上地幔黏滯系數設為5.0×1020Pa·s，計算結果如圖5 所示。可以看到，與圖4 相比，經4.5a的調整后，1998年地震導致F1西段的應力有所減弱；在2020年MW6.0地震發(fā)生之前，經過近20a的震后調節(jié)，F1中段的應力有所減小。震后下地殼和上地幔中的同震應力變化，向上傳遞到上地殼孕震層的過程中，黏滯系數越低，應力傳遞越快(Shanetal.，2013)。而天山的下地殼和上地幔黏滯系數可能比青藏高原更大(盧雙疆等，2013)，因此在本次伽師地震應力計算中，改變黏滯系數對應力分布并無明顯影響。

3.3 改變接收斷層參數對斷層庫侖應力分布的影響

由于未能很好地界定發(fā)震斷層F1和F2的斷層參數，而不同的接收斷層參數會對應力分布的計算結果產生影響，因此，我們參考了前人的反演結果對斷層參數進行了設置。前人利用CDSN的數字波形資料反演得到了2003年地震的震源機制解(走向306°、傾角50°、滑動角129°)(徐錫偉等，2006)；而根據InSAR反演獲得的2020年地震的斷層走向為270°、傾角為15°、滑動角為85°(李成龍等，2020)。以3.2節(jié)中的參數作為地殼和地幔的黏滯系數，計算結果如圖6 所示。從圖中可知，1998年地震加載了F1西段的應力，2003年地震震源處的應力增加>0.01MPa，同時卸載了東段的應力。2003年地震卸載了F1西段的應力，加載了東段的應力，與余震分布的相關性較好；同時，加載了F2東段的應力，在2020年地震震源處的應力改變≤0.01MPa。總體而言，與圖4 和圖5 相比，接收斷層參數的改變對應力分布的影響不大。

圖6 調整斷層參數后發(fā)震斷層上庫侖應力變化Fig.6 Coulomb stress change on the seismogenic faults after adjusting fault parameter.a 1998年地震的同震應力；b 2003年地震發(fā)生之前的應力分布；c 2003年地震發(fā)生之后的應力分布，黑色小圓圈為2003年2月24日—4月30日M>3.5的余震(黃媛等，2006)；d 2020年地震發(fā)生之前的應力分布

3.4 改變地球模型參數對斷層庫侖應力分布的影響

不同的地球模型參數會對應力分布的計算結果產生影響。我們參考了周德敏(2013)在研究柯坪-阿克蘇地震帶時使用的介質模型。由于文中只給出了P波速度，故根據波速比1.73計算S波速度，得到的速度模型如表4 所示。計算應力時采用3.2節(jié)中的黏滯系數。

表 4 柯坪-阿克蘇巖石圈的地層介質模型Table 4 Layered lithosphere medium model in Aksu region

圖7 調整地球模型參數后發(fā)震斷層上的庫侖應力變化Fig.7 Coulomb stress change on the seismogenic faults after adjusting Earth model parameter.a 1998年地震的同震應力；b 2003年地震發(fā)生之前的應力分布；c 2003年地震發(fā)生之后的應力分布，黑色小圓圈為2003年2月24日—4月30日M>3.5的余震(黃媛等，2006)；d 2020年地震發(fā)生前的應力分布

計算結果如圖7 所示。相比于圖4、圖5 和圖6，更改地球模型參數沒有改變應力的整體分布。1998年地震使得2003年地震震源處的庫侖應力增加值>0.01MPa；前2次地震對2020年地震震源處的應力影響介于0.006～0.01MPa之間。震后的黏彈性松弛作用依然不明顯。

4 討論

當斷層面上的庫侖應力增量>0.01MPa，就有可能會觸發(fā)地震，即認為0.01MPa為 “觸發(fā)閾值”(Kingetal.，1994)。本文的計算結果表明，1998年MW6.3地震對斷層F1在2003年地震震中所在區(qū)域附近產生的應力增量為0.01～1MPa，超過了觸發(fā)閾值，表明兩者之間存在觸發(fā)作用。王瓊等(2007)認為1998年地震對2003年地震沒有觸發(fā)作用，產生的應力增加近似為0。兩者的差異可能是斷層破裂模型不同所導致的。王瓊等(2007)基于斷層長度估算公式，根據地震震級構建了破裂模型，并認為滑動是均勻分布的。這種處理方法計算的應力結果對較遠的斷層影響較小。然而2003年地震的發(fā)震斷層F1距離1998年的地震震中很近(圖4)，這樣得到的應力結果誤差可能比較大。本文根據基于地震波形反演得到的破裂滑動分布計算庫侖應力，可相對更真實地反映斷層的應力狀態(tài)。本文的計算結果表明，2020年地震也位于庫侖應力增加區(qū)內，但應力增加值為0.003～0.006MPa，<0.01MPa，說明前2個地震對后一個地震并沒有明顯的觸發(fā)作用，本次地震的發(fā)生可能為其他因素所導致，如區(qū)域應力加載或者以前大地震的黏彈性觸發(fā)等。同時計算結果也顯示，經歷長期應力演化后，2003年地震的發(fā)震斷層F1東段的應力加載依然>0.01MPa，因此應加強對該段的地震危險性監(jiān)測工作。

影響斷層上庫侖應力分布的因素較多，例如地球模型參數、接收斷層參數和地殼黏彈性系數等。本文考慮了不同的影響因素，對應力的分布情況進行了研究。結果顯示，斷層面上的應力在經過長時間的黏彈性松弛后并沒有發(fā)生明顯的變化，這可能是由于地震震級較小導致的。從反演結果來看，前2次地震都為MW6.3地震。根據斷層傾角的計算結果，2個發(fā)震斷層的破裂最大深度均≤20km，且破裂主要發(fā)生在淺部，最大滑動量≤60cm，在斷層模型深部幾乎沒有明顯破裂(圖3)，說明這2個地震對下地殼和地幔的應力擾動有限。黏彈性松弛理論認為，熱的下地殼和地幔的流變作用，使得在震后的數年至數百年時間里，下地殼和上地幔的同震應力場變化傳遞到上地殼的孕震層中，進而影響斷層的力學性質(Pollitz，1992；Marsanetal.，2003；Freedetal.，2007)。由于2個地震在深部產生的同震應力較小，故黏彈性松弛作用也不明顯。因此，在對后續(xù)地震的影響中，由同震產生的應力階變依然占據主導地位。

1997年1月—2003年3月，新疆伽師地區(qū)發(fā)生了一系列強震活動。尤其是自1997年1月21日起，在短短4個月的時間內，該區(qū)域發(fā)生了7次MS≥6.0的大地震，這種現象在中國大陸甚至世界范圍內都是極為罕見的。張竹琪等(2008)計算了1997年強震間的庫侖應力作用，認為這些地震之間存在觸發(fā)關系，且伽師地震的正斷層活動在地震發(fā)生中起到了重要的促進作用，這在一定程度上解釋了伽師地震群的發(fā)震機制。然而，依然有一些問題沒有得到很好解決，例如： 為什么是短時間內發(fā)生多次零散的5級、6級地震，而不是發(fā)生1次7級地震使得能量一次性全部釋放出來？根據應力觸發(fā)理論，這些地震是否是一個早期特大地震(比如1902年阿圖什M8地震)導致的庫侖應力變化后的連續(xù)觸發(fā)？黏彈性應力轉化在這一過程中起到了多大的作用？以上這些問題有待于未來更進一步的研究。

5 結論

本文基于遠場波形反演，重新構建了1998年8月27日與2003年2月24日2次MW6.3地震的破裂滑動分布，并基于震后黏彈性松弛的方法計算了2次地震與2020年MW6.0地震相互之間的應力作用，得到的結論如下：

(1)1998年地震為典型的左旋走滑型地震，震源深度為11.5km，破裂比較集中，斷層的走向為57°，傾角為81°；2003年地震為逆沖兼走滑型地震，震源深度為15.2km，發(fā)震斷層為一條向N低角度傾斜的逆斷裂，斷層的走向為293°，傾角為20°，地震以SE單方向破裂為主。

(2)1998年地震經4.5a的應力調整之后，導致2003年地震發(fā)震斷層的西部應力增加，東部應力卸載，使得2003年地震震源處附近應力增加了0.01～1MPa，前者對后者具有明顯的觸發(fā)作用。前2個地震導致2020年地震的發(fā)震斷層東部應力增加，但震中位置處的應力增加≤0.006MPa，因此對2020年地震不具有明顯的觸發(fā)作用，2020年地震的發(fā)生主要由其他因素所致。在這一過程中，震后黏彈性松弛作用并不明顯，同震應力階變占據了主導地位。

致謝本文利用IRIS地震波形資料進行計算；圖件由GMT軟件生成；使用汪榮江老師的PSGRN/PSCMP程序計算庫侖應力。在此一并表示感謝!