豐成君李 濱李 惠周銘輝張 鵬朱思雨任雅哲戚幫申王苗苗譚成軒陳群策

1.中國地質(zhì)科學院地質(zhì)力學研究所,北京 100081;

2.自然資源部活動構(gòu)造與地質(zhì)安全重點實驗室,北京 100081;

3.中國地質(zhì)調(diào)查局,北京 100037;

4.山東省地震局菏澤地震監(jiān)測中心站,山東 菏澤 274020;

5.中國地質(zhì)大學(北京),北京 100083;

6.長安大學,陜西 西安 710064

0 引言

東喜馬拉雅構(gòu)造結(jié)(簡稱東構(gòu)造結(jié))指喜馬拉雅弧形造山帶東部弧頂部位,位于青藏高原東南段、印度次大陸的東北和緬甸西北角的地帶,是歐亞板塊與印度板塊碰撞、匯聚的南緣地帶(滕吉文等,2006)。始新世,印度板塊和歐亞板塊碰撞擠壓,東構(gòu)造結(jié)作為兩大板塊匯聚擠壓應力最為集中的部位,地殼物質(zhì)強烈變形,發(fā)生了大規(guī)模的北東—南西向水平縮短和垂直增厚(England and Molnar, 1997; 王二七等,2001);上新世,東構(gòu)造結(jié)快速抬升,逐漸形成了以南迦巴瓦峰(海拔7782 m)為隆升中心的“拇指狀”南迦巴瓦構(gòu)造結(jié)(丁林等,1995;張進江等,2003)。

地質(zhì)構(gòu)造上,東構(gòu)造結(jié)由拉薩地塊、印度板塊和緬甸地塊構(gòu)成,區(qū)域內(nèi)分布南喜馬拉雅逆沖推覆褶皺帶、藏南拆離系、雅魯藏布江斷裂帶、嘉黎斷裂帶及阿帕龍斷裂帶等區(qū)域性活動構(gòu)造帶(Schelling and Arita, 1991; Burchfiel et al., 1992; Yin and Harrison, 2000; Ding et al., 2001)。南迦巴瓦周緣還發(fā)育東久-米林、西興拉、墨脫等第四紀地震斷裂(董漢文等,2018;唐方頭等,2019),地震活動強烈,1970—2016年間共發(fā)生3級以上地震2270次,4.5級以上地震234次,6級以上地震27次,大多集中在雅魯藏布江大拐彎頂端,此外,墨脫斷裂帶西南段與喜馬拉雅主中央斷裂的桑構(gòu)造結(jié)結(jié)頂端地震活動也頻發(fā)(楊建亞等,2017)。20世紀以來,最大強震為1950年察隅8.6級地震,震中位于阿帕龍斷裂帶與喜馬拉雅主前緣斷裂的交匯部位,最大地震烈度Ⅻ度,地震滑坡導致雅魯藏布江下游堵江,產(chǎn)生大量堰塞湖(詹慧麗等,2023)。

地殼深部動力過程與淺表層構(gòu)造響應及產(chǎn)生的自然災害,如地震、火山、地震滑坡、巖爆及軟巖大變形等,都與地殼應力密切相關(guān)(王思敬,2002;彭建兵,2006;何滿潮,2014)。區(qū)域構(gòu)造尺度上,地震發(fā)生是在關(guān)鍵構(gòu)造部位地殼應力長期積累、集中、加強并最終導致應變能突然釋放的過程(石耀霖等,2018)。局部山體尺度上,大型山體變形、斜坡失穩(wěn)等均與山體或坡體應力調(diào)整直接相關(guān)(殷躍平,2008)。工程場區(qū)尺度上,深埋隧道面臨的巖爆、軟巖大變形等工程災害一般多在高地應力環(huán)境下發(fā)生(馮夏庭等,2008)。“十四五”期間,隨著東構(gòu)造結(jié)及其周邊地區(qū)的工程規(guī)劃建設(shè),全面掌握南迦巴瓦地區(qū)深淺部地應力場無疑對于區(qū)域地殼穩(wěn)定性研究、交通廊道地質(zhì)安全風險防控具有重要的科學意義。

關(guān)于南迦巴瓦地區(qū)地應力場的研究,交通廊道規(guī)劃建設(shè)之前,主要依靠大地構(gòu)造解析(吳中海等,2015;謝超等,2017)、震源機制反演(王曉楠等,2018;楊帆等,2019)、GPS觀測(唐方頭等,2010; Zhang et al., 2020)、數(shù)值模擬(宋鍵等,2011;孫玉軍等,2017)等方法,重點對區(qū)域構(gòu)造應力場最大主應力方向進行詳細刻畫,基本形成了較確切的認識。交通廊道規(guī)劃建設(shè)之后,逐步開展的原位地應力測量初步揭示了南迦巴瓦周邊地區(qū)重要橋隧區(qū)地應力大小和方向,主要集中在波密至林芝地區(qū)(張鵬等,2017;陳興強等,2022;田朝陽等,2022;張重遠等,2022;張寧等,2022)。然而,在新構(gòu)造活動最為強烈的大拐彎段及雅魯藏布江下游墨脫斷裂帶沿線,有關(guān)地應力場的研究成果幾乎為空白,當前難以有效支撐國家交通廊道規(guī)劃選址與安全運行對地質(zhì)安全風險評估的新需求。

在收集整理南迦巴瓦周邊地區(qū)已有地震震源機制解數(shù)據(jù)的基礎(chǔ)上,通過應力張量反演,揭示南迦巴瓦構(gòu)造結(jié)不同區(qū)段構(gòu)造應力場最大主應力方向,依據(jù)不同活動類型斷層失穩(wěn)臨界地應力條件,聯(lián)合應力形因子和斷層摩擦系數(shù)反演,估算南迦巴瓦周邊不同區(qū)域地應力場大小,并與原位實測結(jié)果進行合理性比較,在此基礎(chǔ)上,采用基于摩爾-庫侖準則的斷層滑動失穩(wěn)概率分析方法,重點探討東久-米林斷裂帶、墨脫斷裂帶及嘉黎斷裂帶中段滑動失穩(wěn)危險性。文章研究方法可為基于深部震源機制解數(shù)據(jù)估算淺層地應力絕對大小提供有效借鑒,研究成果為南迦巴瓦地區(qū)未開展原位地應力測量區(qū)域提供地應力參考,并為交通廊道和清潔能源基地規(guī)劃選址開展區(qū)域構(gòu)造穩(wěn)定性評價提供地應力支撐。

1 活動斷裂與第四紀活動性

文中研究區(qū)主要位于喜馬拉雅東構(gòu)造結(jié)南迦巴瓦周邊地區(qū)(圖1),區(qū)域內(nèi)分布北東向喜馬拉雅主前緣斷裂(F3)、主邊界斷裂(F4)、主中央斷裂(F5),近東西向的雅魯藏布江斷裂帶(F10),以及北西西向的嘉黎斷裂帶(F1)和阿帕龍斷裂帶(F2)等區(qū)域性大型活動斷裂帶(Yin and Harrison, 2000; Ding et al., 2001)。已有研究表明,喜馬拉雅主前緣斷裂(MFT)、主邊界斷裂(MBT)、主中央斷裂(MCT)均為向北傾的低角度(25°～35°)韌 性 逆 沖 斷 層,其 中,MCT自10 Ma以來仍有活動性(尹安,2006),MBT全新世以來活動,沿斷裂兩側(cè)時有地震發(fā)生(Nakata, 1989),MFT為全新世活動斷裂,1934年尼泊爾8.2級地震發(fā)生在該斷裂上,并產(chǎn)生了明顯的地表破裂帶(Kumar et al., 2010; Sapkota et al., 2013)。雅魯藏布江斷裂帶,傾向北北西向,傾角在55°～80°之間,逆沖為主兼走滑運動特征,其東段最新活動時代為晚更新世(武長得等,1990;彭小龍和王道永,2013)。嘉黎斷裂帶中段(易貢-波密)、南段(波密-察隅)晚第四紀以來表現(xiàn)為右旋走滑兼逆沖活動,具有低活動速率、長復發(fā)周期地震特征,發(fā)生M7.0級以上地震危險性高(鐘寧等,2021)。阿帕龍斷裂帶,傾向北東向,全新世活動斷層,具有右旋走滑運動特征,很大可能為1950年察隅8.6級地震的發(fā)震斷層(白玲等,2019)。

此外,南迦巴瓦周緣還分布北東向的東久-米林斷裂帶(F9)、墨脫斷裂帶(F6)和北西向的西興拉斷裂(F7)等主要活動斷裂(圖1)。東久-米林斷裂帶為南迦巴瓦構(gòu)造結(jié)西邊界斷層,從米林以南向東北方向延伸至通麥,東北端與嘉黎斷裂帶交匯,以西傾為主,傾角為60°～70°,全新世以來以左旋走滑兼逆沖活動為主,2017年米林M6.9級地震發(fā)生在斷裂帶北段與北西向的西興拉斷裂交匯部位(葉進等,2020)。墨脫斷裂帶為南迦巴瓦構(gòu)造結(jié)東邊界斷裂,展布在雅魯藏布江東、西兩側(cè),由多條次級斷層斜列而成,傾向南東向,傾角為60°～70°,全新世左旋走滑活動斷裂,其西南段局部還表現(xiàn)出逆沖、正斷運動特征(謝超等,2016;董漢文等,2018)。西興拉斷裂,傾向南西向,傾角60°～70°,全新世活動斷裂,右旋走滑為主兼具逆沖活動,1950年以來,沿該斷裂發(fā)生過2次6.0級以上地震和多次4.0級地震(韋偉等,2018)。

圖1 南迦巴瓦周邊地區(qū)主要活動斷裂帶與中強地震活動分布(1951—2019年)Fig.1 The main active fault belts surrounding the Namcha Barwa syntaxis and the distribution of moderate-strong earthquakes from 1951 to 2019

2 地應力估算模型與方法

對于任意先存斷層(或節(jié)理、裂隙等結(jié)構(gòu)面)(圖2a),在不考慮內(nèi)聚力的前提下,依據(jù)摩爾-庫侖破裂準則可知,當斷層面上的剪應力大于或等于其抗滑阻力時,理論上斷層就會發(fā)生剪切滑動失穩(wěn)(圖2b;Healy et al., 1968; Raleigh et al., 1976):

公式中:σn為斷層面上的正應力;τf為斷層面上的剪應力;μ為斷層面摩擦系數(shù);μσn為斷層面上的抗滑阻力。

依據(jù)斷層滑動失穩(wěn)時的臨界應力摩爾圓(圖2c),可得到斷層面滑動失穩(wěn)狀態(tài)時的剪應力、正應力和摩擦系數(shù)表達式:

圖2 斷層滑動失穩(wěn)力學解析示意圖Fig.2 Sketch map of mechanical analysis on fault slip instability(a) Two-dimensional stress state of different fractures; (b) Frictional sliding on an optimally oriented fault; (c) Mohr diagram corresponding to arbitrary sliding fault σ1-The maximum principal stress; σ3-The minimum principal stress; σn-The normal stress on the fault plane; τf-The shear stress on the fault plane; φ-Internal frictional angel of the fault; β-Angle between the outward normal vector of the fault plane and the direction of σ1

公式中:σ1和σ3分別代表最大和最小主應力大小;φ為斷層內(nèi)摩擦角;其他含義同上。

由公式(2)中斷層面摩擦系數(shù)表達式還可得:

聯(lián)合公式(1)—(3),可得到斷層剪切滑動失穩(wěn)時臨界地應力條件和摩擦系數(shù)的相互關(guān)系:

引入有效應力概念,可將公式(4)中最大、最小主應力替換為最大、最小有效主應力,此時斷層剪切滑動失穩(wěn)時的臨界地應力條件可進一步改寫為(Zoback and Healy, 1992):

公式中:Pp為地層孔隙水壓力(近似等于靜水壓力)。

區(qū)域構(gòu)造應力場通常以水平或近水平作用為主導,在三維主應力空間直角坐標系中,不妨假設(shè)兩個主應力為水平,另外一個主應力垂直且近似等于上覆巖層的重量。依據(jù)Anderson主應力大小和斷層類型分類建議(表1;Zoback, 2007),由公式(5)可得到正斷層、逆沖斷層、走滑斷層滑動失穩(wěn)時的臨界地應力條件,見公式(6)—(8)。

表1 主應力大小與斷層分類Table 1 Principal stress magnitudes and faulting regimes

正斷層滑動失穩(wěn)時的臨界應力狀態(tài):

逆沖斷層滑動失穩(wěn)時的臨界應力狀態(tài):

走滑斷層滑動失穩(wěn)時的臨界應力狀態(tài):

公式中變量含義同表1。

對于單一活動類型或者以某類活動性為主的斷層,當其滑動失穩(wěn)發(fā)生地震時,理論上,在給定斷層臨界摩擦系數(shù),就能限定其滑動失穩(wěn)時的地應力條件。

然而,由公式(6)—(8)僅能估算出有效主應力的比值,暫不能得到最大、最小水平主應力的絕對大小。鑒于此,繼續(xù)引入應力形因子(R)這一關(guān)鍵地應力參數(shù)。在三維主應力直角坐標系下,當斷層面發(fā)生剪切滑動破壞時,垂直于滑動方向的平面內(nèi)剪應力為0,應力形因子可表示為(Gephart and Forsyth, 1984):

公式中:σ1、σ2和σ3分別為最大、中間和最小主應力大小。

由公式(9)進一步得到不同應力類型作用下斷層活動的應力形因子:

① 正斷層狀態(tài)應力形因子

② 逆沖斷層狀態(tài)應力形因子

③ 走滑斷層狀態(tài)應力形因子

確定應力形因子R后,再聯(lián)合公式(6)—(8)就可估算得到斷層滑動失穩(wěn)時的最大、最小水平主應力絕對大小。應力形因子參數(shù)R可通過地震震源機制解數(shù)據(jù)反演得到(萬永革,2015;王璞等,2019)。

3 南迦巴瓦地區(qū)地應力場估測

3.1 震源機制解數(shù)據(jù)與應力參數(shù)反演方法

當前,南迦巴瓦地區(qū)地應力實測數(shù)據(jù)較缺乏,但該區(qū)中小地震頻發(fā),已公開發(fā)表的地震震源機制解數(shù)據(jù)相對充分。利用震源機制解反演地殼應力場已在地質(zhì)、地震、地球物理、地球動力學等領(lǐng)域研究中得到普遍應用。世界地應力圖中基于震源機制解反演得到的應力數(shù)據(jù)占比約54%(Zoback, 1992)。中國大陸及鄰區(qū)現(xiàn)代構(gòu)造應力場基礎(chǔ)數(shù)據(jù)庫中,震源機制解資料占比近70%(謝富仁等,2004)。震源機制解反演和原地應力測量得到的最大主應力方向一致性較好,總體能反映出穩(wěn)定的地殼深部區(qū)域構(gòu)造應力場(Feng et al., 2020)。

參考南迦巴瓦周邊地區(qū)已公開發(fā)表的中小地震震源機制解數(shù)據(jù)(邵翠茹,2009;王曉楠等,2018)、2017年11月18日米林M6.9級(白玲等,2017)、2019年4月23日米林M6.3級等強震數(shù)據(jù)(李國輝等,2020),收集得到研究區(qū)M≥2.0級地震震源機制解數(shù)據(jù)78條(1985年7月至2019年7月)(圖3),其中,M≥2.5級地震的震源機制解數(shù)據(jù)48條,且主要分布在:

圖3 南迦巴瓦周邊地區(qū)M≥2.0級地震震源機制解分布(1985—2019年)Fig.3 Focal mechanism solutions of the M≥2.0 earthquakes surrounding the Namcha Barwa syntaxis from 1985 to 2019

(1)南迦巴瓦峰東北側(cè)的“派鎮(zhèn)-魯朗-易貢-古鄉(xiāng)”區(qū),位于東久-米林斷裂帶和西興拉斷裂交匯部位,2017年11月18日的米林M6.9地震為逆沖型,其余中小地震多以走滑、逆走滑型為主;

(2)南迦巴瓦峰東南側(cè)的“都登-幫辛-布宗-阿帕龍”區(qū),其西邊界為墨脫斷裂帶,北邊界為嘉黎斷裂帶,南邊界為阿帕龍斷裂帶,東南部可延伸至1950年8月15日察隅8.6級地震位置,震源機制解以走滑型為主;

(3)南迦巴瓦峰西南側(cè)的“梅楚卡-都登-格刀”區(qū),位于墨脫斷裂帶與阿帕龍斷裂帶交匯部位,震源機制解以走滑型、逆走滑型為主。

基于震源機制解數(shù)據(jù)反演應力張量時假設(shè):地震滑移沿已知的斷層面發(fā)生,且斷層面為平面;斷層滑動方向(滑動角)與斷層面上的剪應力方向一致;地質(zhì)體為各向同性;區(qū)域地殼應力場相對穩(wěn)定,且參與反演的震源機制解為相互獨立事件。國外學者提出并發(fā)展了幾種基于震源機制解數(shù)據(jù)反演應力張量的方法(Michael, 1984; Gephart and Forsyth, 1984; Lund and Slunga, 1999; Hardebeck and Michael, 2006)。目 前,基 于Michael方法開發(fā)的MSATSI反演程序被廣泛應用于應力張量反演,代入震源機制解樣本數(shù)據(jù),采用網(wǎng)格化的阻尼最小二乘法,可以反演得到不同網(wǎng)格內(nèi)3個主應力的方向、傾角和應力形因子等參數(shù)(Martínez-Garzón et al., 2014)。

然而,采用MSATSI程序代入震源機制解節(jié)面參數(shù)反演時,需要輸入與發(fā)震斷層產(chǎn)狀相一致的主斷面參數(shù),以確保得到相對精確的主應力方向和應力形因子。事實上,僅根據(jù)震源機制解難以限定發(fā)震斷層的優(yōu)勢面,往往需要結(jié)合地質(zhì)調(diào)查、地震觀測或地球物理探測等資料綜合判定,在不能準確輸入優(yōu)勢斷層面參數(shù)情況下,反演得到的主應力方向與區(qū)域構(gòu)造應力場相差甚遠(Feng et al., 2020)。

對此,基于Michael方法,捷克學者Vavrycuk在反演過程中增加了斷層失穩(wěn)參數(shù)I(0≤I≤1;公式(13))用于識別主斷層面,當I=0,斷層相對穩(wěn)定;I=1時,斷層處于滑動失穩(wěn)狀態(tài),此時反演得到的節(jié)理面即為優(yōu)勢發(fā)震斷層面(Vavrycˇuk, 2011):

公式中:σn和τf分別為任意走向斷層面上的正應力和剪應力(圖4中黑圈);σc和τc分別為主斷層面上的正應力和剪應力(圖4中紅圈);μ為摩擦系數(shù);其他含義同上。

圖4 斷層失穩(wěn)參數(shù)在應力莫爾圓中的定義Fig.4 Definition of the fault instability parameter in the Mohr’s diagram σ1-The maximum principal stress;

斷層失穩(wěn)參數(shù)I不依賴于主應力的絕對大小,反演過程中設(shè)定最大主應力σ1為單位應力大小“-1”(負代表壓應力)、最小主應力σ3為“1”(正代表拉應力),由公式(9)得到中間主應力σ2=2R-1。聯(lián)合公式(2)和(3)進一步得到主斷層面上滑動失穩(wěn)時的正應力σc和剪應力τc表達式:

將公式(14)代入(13)中,斷層失穩(wěn)參數(shù)I可改寫為: 任意斷層面上的正應力σn和剪應力τf計算公式為:

公式中:n1、n2、n3分別為斷層面外法線和3個主應力軸夾角的方向余弦,R為應力形因子。

斷層失穩(wěn)參數(shù)對識別發(fā)震斷層面具有重要的物理意義,給定任意斷層面參數(shù)(走向、傾角和滑動角),理論上,均可計算其I值,I值接近“1”的節(jié)面即為發(fā)震斷層面?；谏鲜隼碚?Vavrycˇuk開發(fā)了基于Matlab平臺的聯(lián)合迭代線性反演程序STRESSINVERSE,并以捷克2008年West Bohemia地區(qū)167個中小地震震源機制解為例,比較了與Michael方法反演結(jié)果的差別,結(jié)果顯示,STRESSINVERSE程序反演得到的主應力方向、主斷層面參數(shù)和應力形因子更加穩(wěn)定(Vavry cˇuk, 2014)。此外,在輸入震 源機制解參數(shù)反 演時,也不需要輸入主斷層面參數(shù)。文章將采用STRESSINVERSE程序,反演南迦巴瓦地區(qū)最大主應力方向、斷層摩擦系數(shù)和應力形因子等關(guān)鍵地應力參數(shù)。

3.2 南迦巴瓦周邊地區(qū)最大主應力方向

綜合考慮地震震源機制解分布、區(qū)域內(nèi)北東向和北西西向活動斷裂規(guī)模、活動強度及平面切割關(guān)系等因素,將南迦巴瓦周邊地區(qū)進一步劃分為6個研究區(qū)(圖5):Ⅰ區(qū),“易貢-魯朗-派鎮(zhèn)-古鄉(xiāng)”封閉區(qū)域,位于南迦巴瓦峰西北角;Ⅱ區(qū),“傾多-古鄉(xiāng)-幫辛-那母隆巴-松宗-波密”封閉區(qū)域,位于南迦巴瓦峰東北角;Ⅲ區(qū),“幫辛-都登-格刀-阿帕龍-布宗-那母隆巴”封閉區(qū)域,位于南迦巴瓦峰東南側(cè);Ⅳ區(qū),“格刀-阿帕龍-納姆卡姆”封閉區(qū)域,位于Ⅲ區(qū)南側(cè);Ⅴ區(qū),“都登-米林-梅楚卡-格刀”封閉區(qū)域,位于南迦巴瓦峰西南側(cè);Ⅵ區(qū),由喜馬拉雅主邊界斷裂、主中央斷裂及其他次級斷層組成的圈閉區(qū)域。

圖5 南迦巴瓦周邊地區(qū)震源機制解分區(qū)特征Fig.5 Subarea characteristics of focal mechanism solutions surrounding the Namcha Barwa syntaxis

采用STRESSINVERSE程序,分別將各區(qū)內(nèi)震源機制解數(shù)據(jù)進行整體反演,得到每個應力分區(qū)內(nèi)最大主應力方向(圖6),結(jié)果顯示,Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ區(qū)內(nèi)最大主應力優(yōu)勢方向分別為N28°E、N15°W、N43°E、N35°E、N34°E、N62°W,其中,南迦巴瓦峰北側(cè)的Ⅰ區(qū)和Ⅱ區(qū)為北北東向—近南北向,南側(cè)的Ⅲ、Ⅳ、Ⅴ區(qū)為北東—北北東向。

已有原位地應力測試結(jié)果表明,通麥隧道(位于Ⅰ區(qū)內(nèi))最大水平主應力方向為北東向(陳興強,2022),色季拉山隧道(魯朗南側(cè))為北東—北北東向(田朝陽等,2022),林芝-通麥沿線拉月隧道(位于Ⅰ區(qū)內(nèi))為北東—北東東向(張重遠等,2022),林芝八一鎮(zhèn)為北北東向(平均N26°E)(張鵬等,2017)。文中反演得到Ⅰ區(qū)內(nèi)地殼深部最大主應力方向與淺層原位地應力測量結(jié)果一致性較好。

喜馬拉雅東構(gòu)造結(jié)及周邊地區(qū)GPS應變場和震源機制解反演結(jié)果顯示,喜馬拉雅構(gòu)造帶現(xiàn)今地殼形變以南北向劇烈壓縮和東西向伸展變形為主要特征,南迦巴瓦構(gòu)造結(jié)及周邊地區(qū)最大主壓應力方向自西向東由北北東向偏轉(zhuǎn)至北東向(楊帆等,2019)。文中結(jié)果顯示,從南迦巴瓦峰西北側(cè)的Ⅰ區(qū)至東南側(cè)的Ⅴ區(qū)最大主應力方向基本由北北東向順時針偏轉(zhuǎn)至北東向(圖6),與相關(guān)學者關(guān)于東構(gòu)造結(jié)區(qū)域構(gòu)造應力場的研究結(jié)論相符合。南迦巴瓦周邊地區(qū)北北東—近南北向的最大主壓應力方向反映了印度板塊北向俯沖青藏高原的深部動力過程(張培震等,2002),而東側(cè)逐漸偏轉(zhuǎn)的北東向最大主應力方向一定程度上揭示了青藏高原內(nèi)部地殼物質(zhì)向東南不均勻擠出,并圍繞東構(gòu)造結(jié)產(chǎn)生順時針旋轉(zhuǎn)的運動模式(曹建玲等,2009;吳中海等,2015)。

圖6 南迦巴瓦周邊地區(qū)不同應力區(qū)主應力方向反演結(jié)果Fig.6 Principal stress directions retrieved by the STRESSINVERSE program in different stress areas surrounding the Namcha Barwa syntaxis

反演結(jié)果中,Ⅵ區(qū)內(nèi)最大主應力方向為北北西向或北西西向,與區(qū)域構(gòu)造應力場主壓應力方向不一致。該區(qū)內(nèi)震源機制解樣本數(shù)僅有3個,反演結(jié)果相對離散,可靠性有所降低。區(qū)域內(nèi)最大地震為2019年4月24日的米林M6.3級地震,震源機制解走向/傾角/滑動角分別為202°/17°/20°(李國輝等,2020),基于該斷層面參數(shù)反演得到的主壓應力方向為N16°W(近南北向),由圖3可以看出,該地震發(fā)生在喜馬拉雅主邊界斷裂近似90°大拐彎地區(qū),現(xiàn)今北北西向的最大主壓應力方向可能反映了桑構(gòu)造結(jié)西向俯沖南迦巴瓦構(gòu)造結(jié)和北向推擠的動力過程(白玲等,2019)。

3.3 南迦巴瓦周邊地區(qū)應力形因子與斷層摩擦系數(shù)

Ⅰ—Ⅵ應力分區(qū)內(nèi)應力形因子(R)反演結(jié)果如圖7所示,整體上,南迦巴瓦周邊地區(qū)應力形因子位于0.37～0.75之間,南迦巴瓦峰西側(cè)(0.53～0.59)低于東側(cè)(0.62～0.75)。楊帆等(2019)利用102個M4.7～6.5地震震源機制解反演得到南迦巴瓦構(gòu)造結(jié)周邊地區(qū)應力形因子為0.4～0.8,且顯示該值自西北向東南呈增加趨勢。文中反演得到的南迦巴瓦地區(qū)應力形因子結(jié)果與楊帆等(2019)研究成果基本一致。Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ區(qū)內(nèi)優(yōu)勢應力形因子R值分別為0.59、0.53、0.62、0.70、0.75和0.37。依據(jù)林芝-通麥段地應力測量結(jié)果擬合得到的主應力隨深度變化關(guān)系式(張重遠等,2022),估算得到在1000～2000 m深度范圍內(nèi)應力形因子為0.44～0.65,平均為0.53,與文中Ⅰ區(qū)內(nèi)反演結(jié)果0.59比較接近,一定程度上也表明地殼深淺部應力形因子具有較好的關(guān)聯(lián)性。

圖7 南迦巴瓦周邊地區(qū)不同應力區(qū)應力形因子反演結(jié)果Fig.7 Stress shape ratio retrieved by the STRESSINVERSE program in different stress areas surrounding the Namcha Barwa syntaxis

將Ⅰ—Ⅵ區(qū)內(nèi)每個震源機制解樣本進行單獨反演,還可以得到優(yōu)勢斷層面滑動失穩(wěn)時的臨界摩擦系數(shù)μ值。Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ區(qū)內(nèi)μ值分 別 為0.13～0.32、0.12～0.27、0.12～0.33、0.12～0.32、0.14～0.32、0.12～0.34。潘 前(2021)基于83個地應力測量數(shù)據(jù)計算得到重大鐵路工程拉林段沿線1200 m深度范圍應力積累水平表征參數(shù)μm值[(σ1-σ3)/(σ1+σ3)],其結(jié)果顯示,μm值為0.23～0.30。張寧等(2022)基于青藏高原東南緣900組地應力實測數(shù)據(jù),計算得到墨脫-昌都應力區(qū)內(nèi)μm值為0.28。

由圖2c應力摩爾圓可以看出,μm可以用公式(17)表達:

由公式(17)可以計算得到,川藏鐵路拉林段沿線地殼淺層μ為0.24～0.31,墨脫-昌都應力區(qū)內(nèi)μ為0.29,上述結(jié)果與文中估算得到的南迦巴瓦周邊墨脫-波密地區(qū)(Ⅱ、Ⅲ區(qū))μ最大值0.33基本一致,充分表明文中反演結(jié)果相對合理、可靠。

分析認為,南迦巴瓦周邊地區(qū)主要發(fā)震斷層相對低的摩擦系數(shù)可能與斷裂帶物質(zhì)含較高的黑云母有關(guān)。以往研究表明,黑云母礦物在400℃剪切條件下的平均摩擦系數(shù)一般在0.25左右,最大不超過0.36(任鳳文和何昌榮,2014;路珍和何昌榮,2014)?；A(chǔ)地質(zhì)資料顯示,南迦巴瓦周邊地區(qū)片狀分布含黑云母片巖、片麻巖等構(gòu)造混雜巖帶(潘桂棠等,2020;涂繼耀等,2021)。

3.4 南迦巴瓦周邊地區(qū)地應力大小

綜合考慮南迦巴瓦周邊地區(qū)交通廊道和清潔能源基地規(guī)劃選址和建設(shè)需求,將重點對Ⅰ區(qū)、Ⅱ區(qū)、Ⅲ區(qū)和Ⅴ區(qū)內(nèi)的地應力大小進行估算。

3.4.1 Ⅰ區(qū)地應力大小

Ⅰ區(qū)位于北東向東久-米林斷裂帶(北段)和北西向的西興拉斷裂(西段)交匯部位,第四紀以來,前者以左旋走滑活動為主,后者表現(xiàn)為右旋走滑運動方式,該區(qū)北部還展布北西西向嘉黎斷裂帶,晚第四紀以來表現(xiàn)出右旋走滑活動特征。此外,Ⅰ區(qū)內(nèi)震源機制解類型除2017年11月18日的M6.9級地震為逆沖型外,大多中小地震以走滑型為主,與區(qū)域內(nèi)主要發(fā)震斷裂運動學特征相符合。鑒于此,選取走滑斷層失穩(wěn)臨界應力條件和應力形因子計算公式進行估算。

將Ⅰ區(qū)內(nèi)發(fā)震斷層臨界摩擦系數(shù)的最大值0.32代入公式(8),并將該區(qū)優(yōu)勢應力形因子R=0.59代入公式(12),孔隙水壓力Pp近似等于靜水壓力P0。對于花崗巖、閃長巖和片麻巖,其平均天然密度取2.72 g/cm3,進而估算得到Ⅰ區(qū)內(nèi)最大水平主應力值(σH)線性增加梯度為0.0337 MPa/m,最小水平主應力值(σh)線性增加梯度為0.0227 MPa/m,垂向應力線性增加梯度為0.0272 MPa/m(圖8)。

圖8 Ⅰ區(qū)主應力大小隨深度變化特征Fig.8 Variation characteristics of principal stresses with depth in region Ⅰ

陳興強(2022)基于通麥隧道15個鉆孔地應力測量結(jié)果,擬合得到相對完整巖體內(nèi)最大、最小水平主應力線性增加梯度分別為0.0336 MPa/m 和0.0222 MPa/m,與文中Ⅰ區(qū)估算結(jié)果吻合度高。依據(jù)張重遠等(2022)擬合得到的林芝-通麥段主應力隨深度變化關(guān)系式,估算該區(qū)1000 m深度最大、最小水平主應力值分別為32.95 MPa和23.02 MPa,由圖8擬合關(guān)系式估算相同深度最大、最小水平主應力值分別為33.7 MPa和22.7 MPa,兩者結(jié)果也很相近?？傮w上,文中估算得到的Ⅰ區(qū)內(nèi)地應力大小結(jié)果是合理的。

3.4.2 Ⅱ區(qū)地應力大小

Ⅱ區(qū)內(nèi)主要分布北西西向的嘉黎斷裂帶(中段),晚第四紀以來右旋走滑活動,區(qū)域內(nèi)震源機制解均為走滑型。分別將μ最大值0.27代入公式(8)、R=0.59代入公式(12),估算得到Ⅱ區(qū)主應力大小隨深度變化規(guī)律(圖9),最大水平主應力值(σH)線性增加梯度為0.032 MPa/m,最小水平主應力值(σh)線性增加梯度為0.0229 MPa/m,垂向應力線性增加梯度為0.0272 MPa/m。位于區(qū)內(nèi)的波密縣嘎隆拉隧道鉆孔應力解除法測量結(jié)果顯示,在800 m埋深最大水平主應力值為25.3 MPa(孫曉光等,2011),由圖9估算得到該深度最大水平主應力為25.6 MPa,與嘎隆拉隧道地應力實測結(jié)果基本相符。

圖9 Ⅱ區(qū)主應力大小隨深度變化特征Fig.9 Variation characteristics of principal stresses with depth in region Ⅱ

3.4.3 Ⅲ區(qū)地應力大小

Ⅲ區(qū)內(nèi)主要分布北東向的墨脫斷裂帶,晚第四紀以來以左旋走滑運動為主,該區(qū)內(nèi)震源機制解多表現(xiàn)為走滑型。分別將μ最大值0.33代入公式(8)、R=0.62代入公式(12),估算得到Ⅲ區(qū)主應力大小隨深度變化規(guī)律(圖10),最大水平主應力值(σH)線性增加梯度為0.0344 MPa/m,最小水平主應力值(σh)線性增加梯度為0.0227 MPa/m,垂向應力線性增加梯度為0.0272 MPa/m。根據(jù)張寧等(2022)擬合得到的墨脫-昌都應力區(qū)主應力隨深度變化關(guān)系式,在1000 m深度最大、最小水平主應力值分別為29.89 MPa和17.34 MPa,由圖10中關(guān)系式估算得到該深度最大、最小水平主應力值分別為33.4 MPa和22.7 MPa,該區(qū)水平主應力估算結(jié)果普遍高于墨脫-昌都區(qū)平均應力水平約4～5 MPa。分析認為,墨脫-昌都區(qū)覆蓋面積近6×104km2,將整個區(qū)域內(nèi)的地應力實測數(shù)據(jù)進行統(tǒng)一擬合,往往忽略地應力場的非均勻特征,部分高應力區(qū)難以得到體現(xiàn)。

圖10 Ⅲ區(qū)主應力大小隨深度變化特征Fig.10 Variation characteristics of principal stresses with depth in region Ⅲ

3.4.4 Ⅴ區(qū)地應力大小

Ⅴ區(qū)位于北東向墨脫斷裂帶的西南端,區(qū)域內(nèi)震源機制解有正斷型、逆沖型和走滑型、逆走滑型等,其中走滑型地震占比54%、逆沖型占比20%、正斷型占比20%,綜合上述考慮,該區(qū)仍然選取走滑斷層滑動失穩(wěn)臨界條件進行地應力估算。分別將μ最大值0.32代入公式(8)、R=0.75代入公式(12),估算得到Ⅴ區(qū)主應力大小隨深度變化規(guī)律(圖11),最大水平主應力值(σH)線性增加梯度為0.0355 MPa/m,最小水平主應力值(σh)線性增加梯度為0.0236 MPa/m,垂向應力線性增加梯度為0.0272 MPa/m。由于該區(qū)內(nèi)目前未見相關(guān)地應力測試數(shù)據(jù)公開,暫不能與實測結(jié)果進行可靠性比較。

圖11 Ⅴ區(qū)主應力大小隨深度變化特征Fig.11 Variation characteristics of principal stresses with depth in region Ⅴ

4 南迦巴瓦地區(qū)主要活動斷裂構(gòu)造穩(wěn)定性探討

南迦巴瓦地區(qū)位于印度板塊與歐亞板塊碰撞最強烈的前緣地帶,即喜馬拉雅東構(gòu)造結(jié)雅魯藏布江大拐彎附近,構(gòu)造變形顯著,新構(gòu)造活動強烈,活動斷裂發(fā)育,地震活動頻發(fā)。1950年察隅8.6級地震后,該區(qū)內(nèi)未發(fā)生M7.0級以上強震,多以中小地震為主(楊建亞等,2017),其中最大地震為2017年11月18日米林M6.9級地震,震中烈度達到Ⅷ度,Ⅵ度以上區(qū)域面積近1.3×104km2,涉及巴宜、米林、墨脫、波密及工布江達5個縣(區(qū))(唐方頭等,2019)。

研究表明,1950年察隅8.6級和2017年米林6.9級地震導致周邊區(qū)域應力場重新分布,對喜馬拉雅東構(gòu)造結(jié)周邊活動斷層產(chǎn)生較大影響,庫侖應力變化達到數(shù)兆帕量級(尹鳳玲等,2018;李振月等,2020)。南迦巴瓦構(gòu)造結(jié)周緣主要活動斷裂在現(xiàn)今構(gòu)造應力背景下仍可能產(chǎn)生強震、錯斷等地質(zhì)安全問題,影響高鐵、水電等交通廊道地質(zhì)安全風險防控(蘭恒星等,2021;祁生文等,2022)。相關(guān)案例也屢見不鮮,如2022年1月8日青藏高原東北緣門源縣M6.6級地震,沿冷龍嶺斷裂產(chǎn)生明顯地表破裂帶,地表最大同震位移達3.1 m,發(fā)震斷裂的左旋走滑剪切作用直接導致蘭新高鐵大梁隧道、硫磺溝大橋等設(shè)施嚴重損壞(張威等,2022;韓帥等,2022)。

文章聚焦活動斷裂構(gòu)造穩(wěn)定性問題,采用基于摩爾-庫侖破裂準則的斷層滑動失穩(wěn)概率分析方法(Snee and Zoback, 2018),將上文地應力估算結(jié)果作為邊界條件,定量評價南迦巴瓦周邊Ⅰ區(qū)、Ⅱ區(qū)、Ⅲ區(qū)和Ⅴ區(qū)內(nèi)主要活動斷裂滑動失穩(wěn)風險,以期為交通廊道和清潔能源基地活動斷裂地質(zhì)安全風險評價提供參考。

為便于計算,根據(jù)不同斷裂走向,將4個應力區(qū)內(nèi)的主要活動斷裂簡化為不同長度的直線段(圖12)。參考地質(zhì)、地球物理等資料,合理確定各分段斷層的走向、傾角、長度等參數(shù),同時考慮不同參考資料之間的差異性,分別將各分段斷層的走向和傾角設(shè)定為±5°和±10°誤差范圍(表2)。如上所述,4個區(qū)域內(nèi)反演得到的斷層摩擦系數(shù)與其他研究成果相比,最大殘差為0.08,故將各分段斷層摩擦系數(shù)的誤差范圍設(shè)定為±0.08?；谀?庫侖準則的斷層滑動失穩(wěn)概率分析方法,可詳見作者團隊已發(fā)表論文(朱思雨等,2022;Zhu et al., 2022; Feng et al., 2022; Fan et al., 2022)。

表2 南迦巴瓦周邊地區(qū)主要活動斷裂帶屬性參數(shù)Table 2 Characteristic parameters of the main active fault belts surrounding the Namcha Barwa syntaxis

圖12 南迦巴瓦周邊地區(qū)主要活動斷裂帶簡化分段特征Fig.12 Simplified segments of the main active fault belts surrounding the Namcha Barwa syntaxis

采用美國斯坦福大學開發(fā)的斷層滑動失穩(wěn)概率計算程序FSP1.0(Walsh Ⅲ and Zoback, 2016),定量估算Ⅰ區(qū)、Ⅱ區(qū)、Ⅲ區(qū)和Ⅴ區(qū)內(nèi)主要活動斷裂當前滑動失穩(wěn)概率(圖13)。

圖13 南迦巴瓦周邊地區(qū)主要活動斷裂滑動失穩(wěn)概率Fig.13 Fault slip probability in the main active fault belts surrounding the Namcha Barwa syntaxis

(1)Ⅰ區(qū)內(nèi),嘉黎斷裂帶西北段(F1-4)、西興拉斷裂帶西支(F7-1、F7-2)和東久-米林斷裂帶東北支(F9-3)在當前地應力作用下滑動失穩(wěn)概率很低,一般不超過10%,斷層相對穩(wěn)定,而東久-米林斷裂帶東北段(F9-2、F9-1)滑動失穩(wěn)概率較高,分別為25%和39%,需重點關(guān)注F9-2段未來地震活動危險性。

(2)Ⅱ區(qū)內(nèi),嘉黎斷裂帶中段(F1-1、F1-2和F1-3)滑動失穩(wěn)概率較高,分別為24%、38%和39%,地震活動危險性較高,而墨脫斷裂帶的東北端(F6-9)滑動失穩(wěn)概率很低,小于5%,斷層活動相對穩(wěn)定。

(3)Ⅲ區(qū)內(nèi),墨脫斷裂帶中段(F6-4、F6-5、F6-6、F6-7、F6-8)整體滑動失穩(wěn)概率偏高,位于16%～45%之間,其中,背崩鄉(xiāng)至幫辛鄉(xiāng)段斷層滑動失穩(wěn)概率普遍大于30%,最高達45%(F6-8),需要重點加強對墨脫斷裂帶中段地震活動性的持續(xù)跟蹤與動態(tài)監(jiān)測。阿帕龍斷裂帶滑動失穩(wěn)概率較高,為36%～39%,地震活動危險性高,已有GPS觀測數(shù)據(jù)顯示,阿帕龍斷裂帶可能處于閉鎖階段,為強震危險區(qū)。

(4)Ⅴ區(qū)內(nèi),墨脫斷裂帶南段(F6-1、F6-2、F6-3)當前滑動失穩(wěn)概率很低,均小于8%,相比之下,喜馬拉雅主邊界斷裂北段(F4-1)滑動失穩(wěn)概率很高,最高可達50%,2019年4月24日M6.3級地震就發(fā)生在F4-1斷層南側(cè)約20 km位置,該區(qū)未來仍是中強地震頻發(fā)區(qū)。

此外,從交通廊道地質(zhì)安全風險防控角度來看,交通廊道波密至林芝段在古鄉(xiāng)附近轉(zhuǎn)至與嘉黎斷裂帶(F1-4)平行,由于嘉黎斷裂帶(F1-4)當前滑動失穩(wěn)概率很低,斷層活動穩(wěn)定,可能不會對工程建設(shè)和安全運行產(chǎn)生嚴重影響,比較而言,東久-米林斷裂東北段(F9-2)具有較高的滑動失穩(wěn)概率,后期需深入研究該條斷裂發(fā)生強震時對交通廊道產(chǎn)生的工程錯斷效應。另外,對于區(qū)域內(nèi)水電工程規(guī)劃選址而言,則需要重點評估雅魯藏布江下游墨脫斷裂帶的地震活動危險性及其可能產(chǎn)生的地質(zhì)安全問題。

5 結(jié)論與建議

(1)南迦巴瓦周邊地區(qū)現(xiàn)今構(gòu)造應力場最大主應力方向為北東—北北東向,反映了印度板塊北向俯沖青藏高原的動力背景。

(2)南迦巴瓦周邊地區(qū)地應力大小表現(xiàn)出非均勻特征,最大、最小水平主應力線性增加梯度在西北部的“易貢-魯朗-派鎮(zhèn)-古鄉(xiāng)”區(qū)為0.0337 MPa/m、0.0227 MPa/m,東北部的“古鄉(xiāng)-幫辛-那母隆巴-松宗-波密”區(qū)為0.032 MPa/m、0.0229 MPa/m,東南部的“幫辛-都登-阿帕龍-布宗”區(qū)為0.0344 MPa/m、0.0227 MPa/m,西南部的“都登-米林-梅楚卡-格刀”區(qū)為0.0355 MPa/m、0.0236 MPa/m。

(3)在當前地應力環(huán)境下,南迦巴瓦周邊地區(qū)主要斷裂活動危險性也具有差異性,東久-米林斷裂帶東北段(派鎮(zhèn)至古鄉(xiāng))滑動失穩(wěn)概率為25%～39%,嘉黎斷裂帶中段(古鄉(xiāng)至松宗)為24%～39%,墨脫斷裂帶東北段(背崩至幫辛)為30%～45%,阿帕龍斷裂帶西北段(都登至阿帕龍)為30%～45%,以上區(qū)段均具有較高的地震活動性。水電工程規(guī)劃選址階段,需重點關(guān)注墨脫斷裂帶未來大概率發(fā)生強震并可能引起工程錯斷等地質(zhì)安全問題。

(4)依據(jù)斷層滑動失穩(wěn)臨界應力條件,結(jié)合基于震源機制解數(shù)據(jù)的應力形因子、斷層摩擦系數(shù)等參數(shù)反演,可估算地應力場的優(yōu)勢方向和絕對大小,估算與實測結(jié)果相似性較好,一定程度上可為未開展原位地應力測量的區(qū)域提供補充。估算過程中,設(shè)定的區(qū)域面積、震源機制解樣本數(shù)量、地震震級與類型等因素一定程度上會影響應力形因子和斷層摩擦系數(shù)結(jié)果,因此,文中估算得到的南迦巴瓦周邊地區(qū)地應力場,仍需要結(jié)合區(qū)域內(nèi)原位地應力實測數(shù)據(jù)數(shù)據(jù)進行檢驗和完善。