馮淦, 萬永革,2*, 許鑫, 李梟
1 防災科技學院, 河北三河 065201 2 河北省地震動力學重點實驗室, 河北三河 065201
根據中國地震臺網中心測定,2021年5月22日2時4分,青海果洛州瑪多縣發生MS7.4地震,震中位于北緯34.59°,東經98.34°,震源深度17 km.截至5月23日10時0分,青海地震臺網已記錄到余震共1169條,3級以上余震共36條,其中最大震級為5.1級(北緯34.85°,東經97.5°),距離主震約82 km.
青藏高原是我國現代構造活動與地震活動最強烈的地區(鄧起東等,2002),從動力學角度來看,青藏高原是由于印度洋洋脊擴張使印度板塊向北推擠與歐亞板塊碰撞形成,同時受到塔里木盆地的阻擋,從而導致青藏高原物質東流,極易誘發地震.鄧起東等(1994)根據斷塊的構造活動特征把青藏斷塊分為了8個Ⅱ級斷塊,本次地震發生在巴顏喀拉斷塊,該斷塊四周邊界分別為鮮水河斷裂帶、龍門山斷裂帶、東昆侖斷裂帶、甘孜—玉樹斷裂帶、瑪爾蓋茶卡—若拉崗日斷裂帶(鄧起東等,2010).鄧起東等(2014)收集了1900年以來青藏斷塊7級以上地震,根據地震的時空分布提出了三個地震系列,其中昆侖—汶川地震系列的主體位于巴顏喀拉斷塊.該地震系列的主體包括1997年西藏瑪尼MS7.5地震、2001年昆侖山MS8.1地震、2008年新疆于田MS7.3地震、2008年四川汶川MS8.0地震、2010年青海玉樹MS7.1地震和2017年九寨溝MS7.0地震.為了更好了解青海瑪多MS7.4地震與上述地震的關聯性,我們利用多個機構和作者給出的震源機制解求出青海瑪多MS7.4地震的中心震源機制解(萬永革,2019),結果顯示此次地震為走滑型地震,與上述地震主體以走滑為主的特征相符,同時程佳等(2011)通過計算上述強震產生的庫侖應力變化發現部分強震之間存在一定的觸發作用,而本次瑪多地震的發生會再次引起巴顏喀拉塊體斷裂帶的庫侖應力變化,所以研究該次地震對青海瑪多地區周圍的主要斷層造成的靜態庫侖應力變化有重要意義.
大量的震例研究表明,地震引起的靜態庫侖應力變化,會對周圍地區后續地震的發生起到提前或延遲的作用(Nalbant et al.,1998;萬永革等,2009,2010;McCloskey et al.,2005).如McCloskey等(2005)通過2004年蘇門答臘島MW9.0地震對周圍斷層的靜態應力觸發研究后發現,Nias-Simeulue段是未來可能發生地震的危險區域,2005年發生在Nias的MW8.7地震證實了他們的預測.萬永革等(2009,2010)分別給出了2008年汶川MS8.0地震和于田MS7.3地震在周圍斷層的不同斷層段上產生的庫侖應力變化,并指出部分斷層上的大震的發震被提前或滯后,2013年4月20日發生在汶川地震南部的蘆山7.0級地震和2014年2月12日發生在2008年于田地震震中東北部的7.3級地震就驗證了他們的研究成果.Jia等(2014,2018)計算了2008年汶川MS8.0產生的同震(彈性)與震后流變作用兩部分的影響,結果表明汶川地震分別造成蘆山地震的提前和九寨溝地震的滯后.萬永革等(2015)計算了2015年尼泊爾強震序列在中國大陸地區造成的靜態庫侖應力變化,結果表明該次強震序列產生的應力加載區主要集中在青藏斷塊西南地區的逆沖斷裂和大部分拉張斷裂.靳志同等(2019)研究了2017年九寨溝MS7.0地震對周圍斷層的應力影響,發現該次地震的應力加載區集中在虎牙斷裂中段和南段、岷江斷裂北段南部、塔藏斷裂西段.
青海瑪多MS7.4地震是我國近4年來第一次超過7級的地震,上一次還要追溯到2017年九寨溝MS7.0地震.地震發生后,華俊等(2021)和USGS(2021)分別給出了此次地震的破裂模型,但是USGS模型所采用的震源機制解(決定了破裂面的走向和空間展布)與發震斷層、余震精定位的結果等存在較大的偏差,因此可以認為是非常不可靠的破裂模型,華俊等(2021)基于InSAR數據和余震精定位反演了此次地震的同震滑動分布,更加符合實際破裂情況,所以本研究采用華俊等(2021)提供的破裂模型,基于Okada(1992)給出的適用于半彈性空間的位錯理論,計算了瑪多地震在周圍地區的主要斷層面上產生的庫侖破裂應力變化.其次,為了在沒有準確斷層參數的情況下理解該地震對附近區域的影響,本研究給出了該地震在地震發生深度平面上的面應力和主應力的方向及大小,并結合背景構造應力場以此來分析該次地震對其鄰區的動力學過程造成的影響.
地震的發生是由地下巖石的破裂造成的,Coulomb認為,巖石的抗剪強度受到作用在節面上的剪應力和正應力的影響(Jaeger et al. ,2007),地震專家和學者基于地震活動性觀測和應力變化,提出一次地震的發生會對后續的地震產生影響(應力觸發)(Chinnery,1963;Smith and Van De Lindt,1969;Rybicki,1973).在計算靜態庫侖應力時,主要采用基于彈性半空間的解析表達式(Okada,1992).如果已知發震斷層面的幾何模型和滑動量,則可求出彈性體產生的應力變化張量(Chinnery,1963;陳運泰等,1975),從而把應力變化張量投影到接收斷層面的滑動方向和法向方向,求得庫侖破裂應力變化:
Δσf=Δτ+μ′Δσ,
(1)
式中,Δτ為沿接收斷層面滑動方向的剪切應力變化(為正時,Δτ與滑動方向一致);Δσ為沿接收斷層面法向的正應力變化(拉張為正);μ′為視摩擦系數,由于孔隙流體也會引起應力變化,同時考慮到斷層面上介質特性,一般μ′∈[0.2,0.8](Stein, 1999; Cotton and Coutant, 1997; Wan and Shen, 2010).本研究延續前人的經驗(King et al., 1994; Wan et al., 2003, 2004),μ′的取值為0.4.當Δσf為正時,說明該次地震促進了研究區域斷層的破裂,未來地震發生在研究區域的時間將會提前;當Δσf為負時,說明該次地震抑制了研究區域斷層的破裂(也稱應力影區),未來地震發生在研究區域的時間將滯后.同時研究(Harris, 1998, 2000; Freed, 2005)表明,庫侖應力變化大于0.01 MPa(靜態應力觸發閾值)時能有效地對后續地震的時空分布造成影響.
在計算庫侖破裂應力時,需要接收斷層的精確參數(走向、傾角、滑動角).活動斷層的走向可以根據其在地表的幾何展布獲取,結果較為精確;但是活動斷層的傾角通常只給出一個大概的范圍,精度無法滿足計算的要求;對于活動斷層的滑動角,在地質調查中只給出基礎的滑動性質描述(正(逆)斷層或左旋(右旋)走滑斷層),導致無法獲取一個精確的值.由于以上參數的精確度問題,在計算沿接收斷層面滑動方向的剪應力變化和法向方向的正應力變化會產生誤差,從而導致庫侖破裂應力變化也存在誤差.另外,在研究地震對周圍地區的影響時,水平應力作用(應力變化張量水平分量)也是一個重要因素,這樣可以直觀、簡捷的看出該地震在地圖上的動力學作用特征.
在只考慮應力變化張量水平分量時,應力變化張量可改寫成如下形式:
(2)
面應力則可以表示為
σHc=σHmax+σHmin,
(3)
式中,σHmax表示水平面最大主應力;σHmin表示水平面最小主應力;σHc大于0表示處于膨脹區,反之,表示處于壓縮區.
最大主應力與x軸的夾角:
(4)
由上節內容可知,在研究地震產生的同震庫侖應力變化時,必須先確定位錯理論、發震斷層的幾何模型和滑動量、接收斷層面參數(走向、傾角、滑動角).本研究基于Okada(1992)給出的彈性半空間的位錯理論,采用華俊等(2021)給出的發震斷層面幾何模型和滑動量,接收斷層面參數根據鄧起東等(2002)給出的中國大陸斷層的幾何形狀和運動特征得到(表1).以此來計算該次地震對周圍斷層造成的影響,同時根據地震復發周期和離逝時間為劃分地震危險區給出依據.
由圖1可見:瑪多地震主要造成東昆侖斷裂東段、昆侖山口—江錯斷裂西段、瑪多—甘德斷裂中段、達日斷裂西段西部的庫侖破裂應力增加.其中,在東昆侖斷裂東段西部的庫侖應力增加是本次地震最顯著的區域,最大達171200 Pa(約0.171 MPa),遠超靜態應力觸發閾值(0.01 MPa),東昆侖斷裂東段東部的庫侖應力增加較小,其最東段僅為百帕量級,遠低于靜態應力觸發閾值;昆侖山口—江錯斷裂西段的庫侖應力增加,最大達22820 Pa(約0.023 MPa),超過靜態應力觸發閾值,本次地震就發生在該斷裂上;瑪多—甘德斷裂中段的庫侖應力增加較為顯著,最大達98600 Pa(約0.099 MPa),大于靜態應力觸發閾值;其他斷層的庫侖應力變化詳見表1.

圖1 瑪多地震12 km深度處在周圍主要斷裂上產生的庫侖破裂應力變化填充在斷層上的顏色表示庫侖破裂應力變化的大小,紅色表示增加,藍色表示減小,大小如色標所示.Fig.1 The influence of Coulomb stress change on surrounding faults at 12 km depth of the 2021 Madoi earthquakeThe colors filled on the faults indicate the change of Coulomb stress, which are between red and green for increase, others indicate decrease.

表1 本研究所用的斷層面性質和庫侖破裂應力變化Table 1 Fault segment properties and the change of Coulomb stress
該地震導致東昆侖斷裂中段、昆侖山口—江錯斷裂東段、瑪多—甘德斷裂西段和東段、達日斷裂的庫侖破裂應力降低.其中,東昆侖斷裂中段東部庫侖破裂應力卸載最大達180300 Pa(約0.180 MPa),東昆侖斷裂中段西部庫侖破裂應力卸載僅達到千帕量級;昆侖山口—江錯斷裂東段庫侖破裂應力卸載最為顯著,最大達423500 Pa(約0.424 MPa),本次地震的震中就位于該段;瑪多—甘德斷裂西段庫侖破裂應力卸載最大達296600 Pa(約0.297 MPa),瑪多—甘德斷裂東段庫侖破裂應力卸載由于距震中較遠,僅達到百帕量級;達日斷裂除西段西部外均處于庫侖破裂應力卸載區,斷裂的最東端僅達到百帕量級,其余段均達到104帕量級,最大達52910 Pa(約0.053 MPa).
東昆侖斷裂東段西部庫侖應力增加已超過靜態應力觸發閾值(0.01 MPa),約為0.171 MPa,使得該區域的庫侖應力得到了很大的積累,同時東昆侖破裂帶是一條分割巴顏喀拉斷塊和柴達木斷塊的大型左旋走滑斷裂帶,此外,它還是一條與強震活動關系密切的主要地震構造(劉光勛,1996),但是蔡瑤瑤和張軍龍(2018)通過研究斷裂帶古地震,指出該段的強震復發周期為455±69 a,距上一次強震(1937年7.5級地震)的發生只有84 a,所以東昆侖斷裂東段西部的地震危險性相對較低,東昆侖斷裂東段東部的庫侖破裂應力增加平均為千帕量級,但是根據Ziv和Rubin(2000)的研究,較小的應力變化也會影響地震活動性,邵志剛等(2010)和靳志同等(2019)分別研究了2008汶川地震和2017年九寨溝地震對東昆侖斷裂帶東段的影響,研究結果表明兩次地震均使得該段的庫侖應力增加,同時該段屬于“瑪曲空區”(Wen et al., 2007)范圍,蔡瑤瑤和張軍龍(2018)研究表明該段的強震復發周期為1350±123 a,距上一次強震的發生已1293 a,已經接近強震復發周期,地震危險性相對較高.昆侖山口—江錯斷裂西段的部分區域庫侖應力增加已經超過閾值,但是相關地震地質研究較少,且通過查詢歷史強震目錄發現該斷裂帶從未發生過7級以上的地震,地震危險性需要進一步的研究分析.熊仁偉等(2010)通過野外地質調查發現,瑪多—甘德斷裂自晚第四紀以來可能有過強烈的地質活動且活躍至今,同時任金衛和王敏(2005)通過GPS數據分析出該斷裂正處于應力積累階段,而瑪多—甘德斷裂中段的庫侖應力增加是該次地震較為顯著的地區,最大達0.099 MPa,該段的地震危險性值得重視.
本節基于第1節的方法,計算并討論瑪多地震產生的水平面應力變化對周圍地區的影響.瑪多地震在12 km深度處對周圍地區產生的水平面應力變化如圖2所示,由圖2可見此次地震產生的面應力變化呈現出明顯的四象限分布,地震發生時,震源處會釋放大量的應變能,震中的面應力值下降,導致震中處于面擠壓應力區.由于破裂不均勻,震中西部水平面應力變化絕對值大于東部.震中東北-西南兩側屬于面膨脹應力區,震中附近面膨脹應力最大,約2.63 MPa,并且在向外擴展的過程中數值逐漸減小,到達震中北東面的莊浪河斷裂中段,以及震中南面的怒江斷裂西段時,數值已減小到千帕量級.另外,在震中及震中西北-東南兩側屬于面擠壓應力區,震中附近的面擠壓應力最大,約為3.37 MPa,沿著兩側向外擴展,面擠壓應力逐漸減小,到達震中西北面的大柴旦—宗務隆山斷裂,以及震中東南面的玉樹—瑪曲斷裂東段時,數值已減小到數千帕.

圖2 瑪多地震12 km深度造成的水平面應力變化黑色箭頭和白色箭頭分別表示水平面最小主應力和水平面最大主應力的方向;底色表示面應力,拉張為正;圖中黑色曲線代表水平面應力等值線(單位: MPa);紅色線段代表斷層:F1:東昆侖斷裂,F2:昆侖山口—江錯斷裂,F3:瑪多—甘德斷裂,F4:達日斷裂,F5:怒江斷裂,F6:玉樹—瑪曲斷裂,F7:莊浪河斷裂,F8:大柴旦—宗務隆山斷裂.Fig.2 Horizontal stress changes at 12 km depth generated by the 2021 Madoi earthquakeThe principal compressive stress and principal extensional stress are indicated by the black arrows and white arrows respectively. The background colors represents areal stress. The black curves in the figure represent the areal stress contour (unit: MPa). The red lines represent faults, F1: East Kunlun fault, F2: Kunlunshankou-Jangco fault, F3: Madoi-Gadê fault, F4: Dari fault, F5: Nujiang fault, F6: Yushu-Maqu fault, F7: Zhuanglanghe fault, F8: Da Qaidam-Jun Ul Shan fault.
從水平面最大主應力和最小主應力方向來看,震中附近區域的水平面最小主應力為近東南向,水平面最大主應力方向為近西南向,而Wan(2010)給出的中國現代應力場在34°N,97°E點的主壓應力走向和傾伏角分別為52°和 7°,主張應力走向和傾伏角分別為315°和 44°,基本與本次地震產生的,在震中周圍的水平面最小主應力和水平面最大主應力方向基本相反,這說明本次地震產生的應力場在一定程度上抵消了該區域構造應力場,所以該地震是在區域構造應力場背景下的一次正常應變能釋放.按照萬永革(2020)提出的模擬方法計算瑪多地震在該地區構造應力場中的相對剪應力和正應力(圖3),由圖3可見在該次地震的發震斷層面上相對剪應力和正應力分別為0.831、-0.296,這說明構造應力場在該斷層面的滑動方向上產生了較大的剪應力,此地震的發生釋放了背景構造應力場,使震源附近的背景應力水平降低.

圖3 瑪多地震在區域構造應力場中的相對剪應力(a)和正應力(b)(a)、(b)的底色分別表示相對剪應力和正應力,大小如色標所示. 震源機制類型圖例如子圖(b)下方所示.Fig.3 The relative shear stress and relative normal stress on the tectonic stress field of the 2021 Madoi earthquakeThe background colors of (a) and (b) represent the relative shear stress and normal stress, respectively. The legend of type of focal mechanism is shown in the bottom of (b).
在面膨脹應力區,水平面最大主應力以震中為中心呈輻射狀以向外擴散,而水平面最小主應力以震中為中心呈同心圓分布;在面擠壓應力區,水平面最小主應力以震中為中心呈輻射狀向外擴散,而水平面最大主應力以震中為中心呈同心圓分布,從整體上來看,水平面最大主應力和水平面最小主應力的分布類似于磁場線的分布,并且兩者相互垂直.物質在面擠壓應力區的運動與水平面最小主應力方向大體一致,而在面膨脹應力區的運動與水平面最大主應力方向大體一致,通過分析水平面最大主應力和水平面最小主應力的方向可以大致判斷出物質在研究區域的運動方向,與圖4的結果一致.
瑪多地震發生后不久,通過收集不同機構和作者給出的震源機制解求出中心震源機制解(萬永革,2019),結果顯示該次地震的發震斷層面走向為101.68°,傾角為81.62°,滑動角為-4.19°,同時根據Zoback(1992)給出的震源機制劃分表判斷該次地震為左旋走滑型地震.從圖2中可以看出面擠壓應力區與面膨脹應力區的分界線與震源機制兩個節面的走向較為一致,而面膨脹應力區與面擠壓應力區是震源機制伸張區與壓縮區的較好延伸,說明瑪多地震的水平應力變化符合左旋走滑型地震的特征.
為了能夠更好的理解瑪多地震產生的應力場結果,本文計算了瑪多地震在12 km深度產生的位移場.結果如圖4,圖中可以看出震中西南和東北兩側物質向震中匯聚,且距震中越遠位移量越小,說明區域內拉張作用明顯,從垂直位移來看,震中西南和東北兩側表現出明顯的隆升.同理,震中東南和西北兩側物質向外流出,距震中越遠位移量越小,說明區域內擠壓作用明顯,從垂直位移來看,震中東南和西北兩側表現出明顯的沉降.同時水平位移大于垂直位移,與左旋走滑型地震產生的位移場吻合.自該位移場的分布也可以看出與圖2給出的該地震產生的主應力方向和面應力分布一致.

圖4 瑪多地震在12 km深度產生的位移場圖中紅色箭頭表示位移場水平分量,底色表示位移場垂直分量.Fig.4 The displacement field at 12 km depth generated by the Madoi earthquakeThe red arrows show the horizontal displacement field and the background colors show the vertical displacement field.
本研究基于彈性半空間位錯理論,采用華俊等(2021)提供的瑪多地震破裂模型,計算該地震在周圍地區主要斷層產生的同震庫侖應力變化以及在周圍地區產生的水平應力變化,得到以下結論:
(1) 2021年青海瑪多地震造成東昆侖斷裂東段西部、瑪多—甘德斷裂中段、昆侖山口—江錯斷裂西段庫侖破裂應力增加,最大值分別為0.171 MPa、0.099 MPa、0.023 MPa,均超過靜態應力觸發閾值(0.01 MPa),起到明顯的觸發作用,另外造成東昆侖斷裂中段東部、昆侖山口—江錯斷裂東段、瑪多—甘德斷裂西段、達日斷裂庫侖破裂應力下降,最大值分別為0.180 MPa、0.424 MPa、0.297 MPa、0.053 MPa,應力卸載效果明顯.雖然東昆侖斷裂東段東部平均庫侖破裂應力增加僅為數千帕,但是該段屬于‘瑪曲地震空區’,同時受到2008年汶川地震和2017年九寨溝地震的影響(邵志剛,2010;靳志同,2019),并且接近地震復發周期,地震危險性也同樣值得注意.
(2) 2021年青海瑪多地震造成東昆侖—柴達木斷塊西部和巴顏喀拉斷塊東部的面應力上升,而使得東昆侖—柴達木斷塊東部和巴顏喀拉斷塊西部的面應力下降,在地圖上呈現出明顯的四象限分布,但是由于破裂不均勻,震中西部面應力變化的絕對值大于東部.從主張應力和主壓應力的方向來看,震中附近區域的主壓應力為近東南向,主張應力方向為近西南向,與Wan(2010)得到的中國現代應力場的主壓應力和主張應力的走向基本相反,根據萬永革(2020)提出的方法計算了瑪多地震在該地區構造應力場中的相對剪應力(0.831)和正應力(-0.296),結果表明本次地震基本是在背景構造應力場作用下的一次正常釋放,其動力來源主要是印度洋板塊的NE擠壓和太平洋板塊的西向俯沖.在面膨脹應力區,主張應力以震中為中心呈輻射狀向外擴散,而主壓應力以震中為中心呈磁場線分布,在面擠壓應力區與之相反.通過以上特征可以判斷出物質在面膨脹應力區的運動近北東向,而在面擠壓應力區的運動近西南向.
本文在計算瑪多地震產生的水平面應力變化時基于半無限均勻彈性空間模型,但是在計算中是不考慮介質體存在背景構造應力場的.所以本文根據Wan(2010)給出的中國現代應力場的主張應力軸、中間應力軸、主壓應力軸的方向和R值(取值范圍:32°N—38°N,95°E—103°E),假設主張應力大小為5 MPa,利用matlab中的interp2函數進行二維數據內插(插值方法采用雙三次插值:bicubic),計算研究區域的構造應力場(圖5),并以此討論本次地震產生的水平面應力變化對構造應力場的影響.

圖5 瑪多地震震前(a)及震后(b)構造應力場的水平分量黑色箭頭和白色箭頭分別表示水平面最小主應力和最大主應力的方向;底色表示面應力.Fig.5 The horizontal component of tectonic stress field before (a) and after (b) the 2021 Madoi earthquakeThe principal compressive stress and principal extensional stress are indicated by the black arrows and white arrows respectively. The background colors represents areal stress.
圖5中可以看出在考慮本次地震產生的應力變化后,構造應力場的水平面應力僅在震中附近(34°N—35°N,97.6°E—99.2°E)產生了明顯的變化,而在遠離震中處引起的變化幾乎可以忽略不計.從構造應力場水平分量的最大主應力和最小主應力的方向來看,最大偏轉量僅達2°(34.6°N,98.6°E),這與萬永革等(2006)給出的走滑大地震造成的構造應力場的應力軸偏轉的推測一致.同時需要注意的是,在計算研究區域構造應力時,主張應力假設為5 MPa,隨著主張應力值的逐漸擴大(比如數十兆帕),本次地震對構造應力場的影響會逐漸減弱.這從另一方面說明背景構造應力場是相對非常穩定的,即使是大地震,也只能改變地震震源附近的小范圍內的背景構造應力場,不能改變構造應力場的大趨勢分布.
本研究采用彈性半空間模型對青海瑪多地震震中附近主要斷層的庫侖破裂應力變化及水平面應力的變化進行了研究,但實際的地球地下介質不是均勻的、各向同性的完全彈性體,計算結果會產生一定的誤差.另外,大地震后,上地幔和下地殼的黏彈性松弛效應會產生應變的擴散(Robinson and Zhou, 2005; 萬永革等, 2007, 2008),但是這種黏彈性松弛效應在震后較短時間內的影響可以忽略不計,只對時間跨度較大的地震起作用.震中附近的1937年花石峽地震,1947年的達日地震便可通過計算黏彈性庫侖應力變化判斷是否對本次地震起到觸發作用,但是相關資料較少,無法達到精確計算的要求.
在研究靜態應力觸發的過程中,主震對余震的應力觸發得到了廣泛的關注,Toda等(2011)通過計算2011年日本地震的靜態庫侖應力變化,發現大部分的余震位于主震的庫侖應力增加區,但是并不是所有的結果都盡如人意,如Astiz等(2000)研究了1990年Upland地震對余震的觸發,發現余震的分布與主震的觸發區域并不是很吻合.另外,在計算主震對余震的應力觸發時,通常假設余震發生在最優破裂面上,但是最優破裂面的選取缺少嚴格的論證,同時由于缺少本次地震的余震精確定位和震源機制,所以本研究并未討論該地震對余震的觸發作用.
盡管有上述不確定性,但本文在一階近似下給出了該地震破裂對周圍斷層的影響,并計算了地震發生深度平面上的面應力和主應力的方向及大小.這些對于該地區的地震動力學研究是有意義的.
致謝中國地震局地質研究所華俊等(2021)為本研究提供了瑪多地震的破裂模型,本研究繪圖采用GMT軟件(Wessel et al., 1995)繪制,3位審稿專家為本文提出了建設性修改建議,明顯增加了本文的邏輯性和完整性.特此致謝!