南迦巴瓦地区地应力场估算与构造稳定性探讨

丰成君李 滨李 惠周铭辉张 鹏朱思雨任雅哲戚帮申王苗苗谭成轩陈群策

1.中国地质科学院地质力学研究所,北京 100081;

2.自然资源部活动构造与地质安全重点实验室,北京 100081;

3.中国地质调查局,北京 100037;

4.山东省地震局菏泽地震监测中心站,山东 菏泽 274020;

5.中国地质大学(北京),北京 100083;

6.长安大学,陕西 西安 710064

0 引言

东喜马拉雅构造结(简称东构造结)指喜马拉雅弧形造山带东部弧顶部位,位于青藏高原东南段、印度次大陆的东北和缅甸西北角的地带,是欧亚板块与印度板块碰撞、汇聚的南缘地带(滕吉文等,2006)。始新世,印度板块和欧亚板块碰撞挤压,东构造结作为两大板块汇聚挤压应力最为集中的部位,地壳物质强烈变形,发生了大规模的北东—南西向水平缩短和垂直增厚(England and Molnar, 1997; 王二七等,2001);上新世,东构造结快速抬升,逐渐形成了以南迦巴瓦峰(海拔7782 m)为隆升中心的“拇指状”南迦巴瓦构造结(丁林等,1995;张进江等,2003)。

地质构造上,东构造结由拉萨地块、印度板块和缅甸地块构成,区域内分布南喜马拉雅逆冲推覆褶皱带、藏南拆离系、雅鲁藏布江断裂带、嘉黎断裂带及阿帕龙断裂带等区域性活动构造带(Schelling and Arita, 1991; Burchfiel et al., 1992; Yin and Harrison, 2000; Ding et al., 2001)。南迦巴瓦周缘还发育东久-米林、西兴拉、墨脱等第四纪地震断裂(董汉文等,2018;唐方头等,2019),地震活动强烈,1970—2016年间共发生3级以上地震2270次,4.5级以上地震234次,6级以上地震27次,大多集中在雅鲁藏布江大拐弯顶端,此外,墨脱断裂带西南段与喜马拉雅主中央断裂的桑构造结结顶端地震活动也频发(杨建亚等,2017)。20世纪以来,最大强震为1950年察隅8.6级地震,震中位于阿帕龙断裂带与喜马拉雅主前缘断裂的交汇部位,最大地震烈度Ⅻ度,地震滑坡导致雅鲁藏布江下游堵江,产生大量堰塞湖(詹慧丽等,2023)。

地壳深部动力过程与浅表层构造响应及产生的自然灾害,如地震、火山、地震滑坡、岩爆及软岩大变形等,都与地壳应力密切相关(王思敬,2002;彭建兵,2006;何满潮,2014)。区域构造尺度上,地震发生是在关键构造部位地壳应力长期积累、集中、加强并最终导致应变能突然释放的过程(石耀霖等,2018)。局部山体尺度上,大型山体变形、斜坡失稳等均与山体或坡体应力调整直接相关(殷跃平,2008)。工程场区尺度上,深埋隧道面临的岩爆、软岩大变形等工程灾害一般多在高地应力环境下发生(冯夏庭等,2008)。“十四五”期间,随着东构造结及其周边地区的工程规划建设,全面掌握南迦巴瓦地区深浅部地应力场无疑对于区域地壳稳定性研究、交通廊道地质安全风险防控具有重要的科学意义。

关于南迦巴瓦地区地应力场的研究,交通廊道规划建设之前,主要依靠大地构造解析(吴中海等,2015;谢超等,2017)、震源机制反演(王晓楠等,2018;杨帆等,2019)、GPS观测(唐方头等,2010; Zhang et al., 2020)、数值模拟(宋键等,2011;孙玉军等,2017)等方法,重点对区域构造应力场最大主应力方向进行详细刻画,基本形成了较确切的认识。交通廊道规划建设之后,逐步开展的原位地应力测量初步揭示了南迦巴瓦周边地区重要桥隧区地应力大小和方向,主要集中在波密至林芝地区(张鹏等,2017;陈兴强等,2022;田朝阳等,2022;张重远等,2022;张宁等,2022)。然而,在新构造活动最为强烈的大拐弯段及雅鲁藏布江下游墨脱断裂带沿线,有关地应力场的研究成果几乎为空白,当前难以有效支撑国家交通廊道规划选址与安全运行对地质安全风险评估的新需求。

在收集整理南迦巴瓦周边地区已有地震震源机制解数据的基础上,通过应力张量反演,揭示南迦巴瓦构造结不同区段构造应力场最大主应力方向,依据不同活动类型断层失稳临界地应力条件,联合应力形因子和断层摩擦系数反演,估算南迦巴瓦周边不同区域地应力场大小,并与原位实测结果进行合理性比较,在此基础上,采用基于摩尔-库仑准则的断层滑动失稳概率分析方法,重点探讨东久-米林断裂带、墨脱断裂带及嘉黎断裂带中段滑动失稳危险性。文章研究方法可为基于深部震源机制解数据估算浅层地应力绝对大小提供有效借鉴,研究成果为南迦巴瓦地区未开展原位地应力测量区域提供地应力参考,并为交通廊道和清洁能源基地规划选址开展区域构造稳定性评价提供地应力支撑。

1 活动断裂与第四纪活动性

文中研究区主要位于喜马拉雅东构造结南迦巴瓦周边地区(图1),区域内分布北东向喜马拉雅主前缘断裂(F3)、主边界断裂(F4)、主中央断裂(F5),近东西向的雅鲁藏布江断裂带(F10),以及北西西向的嘉黎断裂带(F1)和阿帕龙断裂带(F2)等区域性大型活动断裂带(Yin and Harrison, 2000; Ding et al., 2001)。已有研究表明,喜马拉雅主前缘断裂(MFT)、主边界断裂(MBT)、主中央断裂(MCT)均为向北倾的低角度(25°～35°)韧 性 逆 冲 断 层,其 中,MCT自10 Ma以来仍有活动性(尹安,2006),MBT全新世以来活动,沿断裂两侧时有地震发生(Nakata, 1989),MFT为全新世活动断裂,1934年尼泊尔8.2级地震发生在该断裂上,并产生了明显的地表破裂带(Kumar et al., 2010; Sapkota et al., 2013)。雅鲁藏布江断裂带,倾向北北西向,倾角在55°～80°之间,逆冲为主兼走滑运动特征,其东段最新活动时代为晚更新世(武长得等,1990;彭小龙和王道永,2013)。嘉黎断裂带中段(易贡-波密)、南段(波密-察隅)晚第四纪以来表现为右旋走滑兼逆冲活动,具有低活动速率、长复发周期地震特征,发生M7.0级以上地震危险性高(钟宁等,2021)。阿帕龙断裂带,倾向北东向,全新世活动断层,具有右旋走滑运动特征,很大可能为1950年察隅8.6级地震的发震断层(白玲等,2019)。

此外,南迦巴瓦周缘还分布北东向的东久-米林断裂带(F9)、墨脱断裂带(F6)和北西向的西兴拉断裂(F7)等主要活动断裂(图1)。东久-米林断裂带为南迦巴瓦构造结西边界断层,从米林以南向东北方向延伸至通麦,东北端与嘉黎断裂带交汇,以西倾为主,倾角为60°～70°,全新世以来以左旋走滑兼逆冲活动为主,2017年米林M6.9级地震发生在断裂带北段与北西向的西兴拉断裂交汇部位(叶进等,2020)。墨脱断裂带为南迦巴瓦构造结东边界断裂,展布在雅鲁藏布江东、西两侧,由多条次级断层斜列而成,倾向南东向,倾角为60°～70°,全新世左旋走滑活动断裂,其西南段局部还表现出逆冲、正断运动特征(谢超等,2016;董汉文等,2018)。西兴拉断裂,倾向南西向,倾角60°～70°,全新世活动断裂,右旋走滑为主兼具逆冲活动,1950年以来,沿该断裂发生过2次6.0级以上地震和多次4.0级地震(韦伟等,2018)。

图1 南迦巴瓦周边地区主要活动断裂带与中强地震活动分布(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 地应力估算模型与方法

对于任意先存断层(或节理、裂隙等结构面)(图2a),在不考虑内聚力的前提下,依据摩尔-库仑破裂准则可知,当断层面上的剪应力大于或等于其抗滑阻力时,理论上断层就会发生剪切滑动失稳(图2b;Healy et al., 1968; Raleigh et al., 1976):

公式中:σn为断层面上的正应力;τf为断层面上的剪应力;μ为断层面摩擦系数;μσn为断层面上的抗滑阻力。

依据断层滑动失稳时的临界应力摩尔圆(图2c),可得到断层面滑动失稳状态时的剪应力、正应力和摩擦系数表达式:

图2 断层滑动失稳力学解析示意图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分别代表最大和最小主应力大小;φ为断层内摩擦角;其他含义同上。

由公式(2)中断层面摩擦系数表达式还可得:

联合公式(1)—(3),可得到断层剪切滑动失稳时临界地应力条件和摩擦系数的相互关系:

引入有效应力概念,可将公式(4)中最大、最小主应力替换为最大、最小有效主应力,此时断层剪切滑动失稳时的临界地应力条件可进一步改写为(Zoback and Healy, 1992):

公式中:Pp为地层孔隙水压力(近似等于静水压力)。

区域构造应力场通常以水平或近水平作用为主导,在三维主应力空间直角坐标系中,不妨假设两个主应力为水平,另外一个主应力垂直且近似等于上覆岩层的重量。依据Anderson主应力大小和断层类型分类建议(表1;Zoback, 2007),由公式(5)可得到正断层、逆冲断层、走滑断层滑动失稳时的临界地应力条件,见公式(6)—(8)。

表1 主应力大小与断层分类Table 1 Principal stress magnitudes and faulting regimes

正断层滑动失稳时的临界应力状态:

逆冲断层滑动失稳时的临界应力状态:

走滑断层滑动失稳时的临界应力状态:

公式中变量含义同表1。

对于单一活动类型或者以某类活动性为主的断层,当其滑动失稳发生地震时,理论上,在给定断层临界摩擦系数,就能限定其滑动失稳时的地应力条件。

然而,由公式(6)—(8)仅能估算出有效主应力的比值,暂不能得到最大、最小水平主应力的绝对大小。鉴于此,继续引入应力形因子(R)这一关键地应力参数。在三维主应力直角坐标系下,当断层面发生剪切滑动破坏时,垂直于滑动方向的平面内剪应力为0,应力形因子可表示为(Gephart and Forsyth, 1984):

公式中:σ1、σ2和σ3分别为最大、中间和最小主应力大小。

由公式(9)进一步得到不同应力类型作用下断层活动的应力形因子:

① 正断层状态应力形因子

② 逆冲断层状态应力形因子

③ 走滑断层状态应力形因子

确定应力形因子R后,再联合公式(6)—(8)就可估算得到断层滑动失稳时的最大、最小水平主应力绝对大小。应力形因子参数R可通过地震震源机制解数据反演得到(万永革,2015;王璞等,2019)。

3 南迦巴瓦地区地应力场估测

3.1 震源机制解数据与应力参数反演方法

当前,南迦巴瓦地区地应力实测数据较缺乏,但该区中小地震频发,已公开发表的地震震源机制解数据相对充分。利用震源机制解反演地壳应力场已在地质、地震、地球物理、地球动力学等领域研究中得到普遍应用。世界地应力图中基于震源机制解反演得到的应力数据占比约54%(Zoback, 1992)。中国大陆及邻区现代构造应力场基础数据库中,震源机制解资料占比近70%(谢富仁等,2004)。震源机制解反演和原地应力测量得到的最大主应力方向一致性较好,总体能反映出稳定的地壳深部区域构造应力场(Feng et al., 2020)。

参考南迦巴瓦周边地区已公开发表的中小地震震源机制解数据(邵翠茹,2009;王晓楠等,2018)、2017年11月18日米林M6.9级(白玲等,2017)、2019年4月23日米林M6.3级等强震数据(李国辉等,2020),收集得到研究区M≥2.0级地震震源机制解数据78条(1985年7月至2019年7月)(图3),其中,M≥2.5级地震的震源机制解数据48条,且主要分布在:

图3 南迦巴瓦周边地区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)南迦巴瓦峰东北侧的“派镇-鲁朗-易贡-古乡”区,位于东久-米林断裂带和西兴拉断裂交汇部位,2017年11月18日的米林M6.9地震为逆冲型,其余中小地震多以走滑、逆走滑型为主;

(2)南迦巴瓦峰东南侧的“都登-帮辛-布宗-阿帕龙”区,其西边界为墨脱断裂带,北边界为嘉黎断裂带,南边界为阿帕龙断裂带,东南部可延伸至1950年8月15日察隅8.6级地震位置,震源机制解以走滑型为主;

(3)南迦巴瓦峰西南侧的“梅楚卡-都登-格刀”区,位于墨脱断裂带与阿帕龙断裂带交汇部位,震源机制解以走滑型、逆走滑型为主。

基于震源机制解数据反演应力张量时假设:地震滑移沿已知的断层面发生,且断层面为平面;断层滑动方向(滑动角)与断层面上的剪应力方向一致;地质体为各向同性;区域地壳应力场相对稳定,且参与反演的震源机制解为相互独立事件。国外学者提出并发展了几种基于震源机制解数据反演应力张量的方法(Michael, 1984; Gephart and Forsyth, 1984; Lund and Slunga, 1999; Hardebeck and Michael, 2006)。目 前,基 于Michael方法开发的MSATSI反演程序被广泛应用于应力张量反演,代入震源机制解样本数据,采用网格化的阻尼最小二乘法,可以反演得到不同网格内3个主应力的方向、倾角和应力形因子等参数(Martínez-Garzón et al., 2014)。

然而,采用MSATSI程序代入震源机制解节面参数反演时,需要输入与发震断层产状相一致的主断面参数,以确保得到相对精确的主应力方向和应力形因子。事实上,仅根据震源机制解难以限定发震断层的优势面,往往需要结合地质调查、地震观测或地球物理探测等资料综合判定,在不能准确输入优势断层面参数情况下,反演得到的主应力方向与区域构造应力场相差甚远(Feng et al., 2020)。

对此,基于Michael方法,捷克学者Vavrycuk在反演过程中增加了断层失稳参数I(0≤I≤1;公式(13))用于识别主断层面,当I=0,断层相对稳定;I=1时,断层处于滑动失稳状态,此时反演得到的节理面即为优势发震断层面(Vavrycˇuk, 2011):

公式中:σn和τf分别为任意走向断层面上的正应力和剪应力(图4中黑圈);σc和τc分别为主断层面上的正应力和剪应力(图4中红圈);μ为摩擦系数;其他含义同上。

图4 断层失稳参数在应力莫尔圆中的定义Fig.4 Definition of the fault instability parameter in the Mohr’s diagram σ1-The maximum principal stress;

断层失稳参数I不依赖于主应力的绝对大小,反演过程中设定最大主应力σ1为单位应力大小“-1”(负代表压应力)、最小主应力σ3为“1”(正代表拉应力),由公式(9)得到中间主应力σ2=2R-1。联合公式(2)和(3)进一步得到主断层面上滑动失稳时的正应力σc和剪应力τc表达式:

将公式(14)代入(13)中,断层失稳参数I可改写为: 任意断层面上的正应力σn和剪应力τf计算公式为:

公式中:n1、n2、n3分别为断层面外法线和3个主应力轴夹角的方向余弦,R为应力形因子。

断层失稳参数对识别发震断层面具有重要的物理意义,给定任意断层面参数(走向、倾角和滑动角),理论上,均可计算其I值,I值接近“1”的节面即为发震断层面。基于上述理论,Vavrycˇuk开发了基于Matlab平台的联合迭代线性反演程序STRESSINVERSE,并以捷克2008年West Bohemia地区167个中小地震震源机制解为例,比较了与Michael方法反演结果的差别,结果显示,STRESSINVERSE程序反演得到的主应力方向、主断层面参数和应力形因子更加稳定(Vavry cˇuk, 2014)。此外,在输入震 源机制解参数反 演时,也不需要输入主断层面参数。文章将采用STRESSINVERSE程序,反演南迦巴瓦地区最大主应力方向、断层摩擦系数和应力形因子等关键地应力参数。

3.2 南迦巴瓦周边地区最大主应力方向

综合考虑地震震源机制解分布、区域内北东向和北西西向活动断裂规模、活动强度及平面切割关系等因素,将南迦巴瓦周边地区进一步划分为6个研究区(图5):Ⅰ区,“易贡-鲁朗-派镇-古乡”封闭区域,位于南迦巴瓦峰西北角;Ⅱ区,“倾多-古乡-帮辛-那母隆巴-松宗-波密”封闭区域,位于南迦巴瓦峰东北角;Ⅲ区,“帮辛-都登-格刀-阿帕龙-布宗-那母隆巴”封闭区域,位于南迦巴瓦峰东南侧;Ⅳ区,“格刀-阿帕龙-纳姆卡姆”封闭区域,位于Ⅲ区南侧;Ⅴ区,“都登-米林-梅楚卡-格刀”封闭区域,位于南迦巴瓦峰西南侧;Ⅵ区,由喜马拉雅主边界断裂、主中央断裂及其他次级断层组成的圈闭区域。

图5 南迦巴瓦周边地区震源机制解分区特征Fig.5 Subarea characteristics of focal mechanism solutions surrounding the Namcha Barwa syntaxis

采用STRESSINVERSE程序,分别将各区内震源机制解数据进行整体反演,得到每个应力分区内最大主应力方向(图6),结果显示,Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ区内最大主应力优势方向分别为N28°E、N15°W、N43°E、N35°E、N34°E、N62°W,其中,南迦巴瓦峰北侧的Ⅰ区和Ⅱ区为北北东向—近南北向,南侧的Ⅲ、Ⅳ、Ⅴ区为北东—北北东向。

已有原位地应力测试结果表明,通麦隧道(位于Ⅰ区内)最大水平主应力方向为北东向(陈兴强,2022),色季拉山隧道(鲁朗南侧)为北东—北北东向(田朝阳等,2022),林芝-通麦沿线拉月隧道(位于Ⅰ区内)为北东—北东东向(张重远等,2022),林芝八一镇为北北东向(平均N26°E)(张鹏等,2017)。文中反演得到Ⅰ区内地壳深部最大主应力方向与浅层原位地应力测量结果一致性较好。

喜马拉雅东构造结及周边地区GPS应变场和震源机制解反演结果显示,喜马拉雅构造带现今地壳形变以南北向剧烈压缩和东西向伸展变形为主要特征,南迦巴瓦构造结及周边地区最大主压应力方向自西向东由北北东向偏转至北东向(杨帆等,2019)。文中结果显示,从南迦巴瓦峰西北侧的Ⅰ区至东南侧的Ⅴ区最大主应力方向基本由北北东向顺时针偏转至北东向(图6),与相关学者关于东构造结区域构造应力场的研究结论相符合。南迦巴瓦周边地区北北东—近南北向的最大主压应力方向反映了印度板块北向俯冲青藏高原的深部动力过程(张培震等,2002),而东侧逐渐偏转的北东向最大主应力方向一定程度上揭示了青藏高原内部地壳物质向东南不均匀挤出,并围绕东构造结产生顺时针旋转的运动模式(曹建玲等,2009;吴中海等,2015)。

图6 南迦巴瓦周边地区不同应力区主应力方向反演结果Fig.6 Principal stress directions retrieved by the STRESSINVERSE program in different stress areas surrounding the Namcha Barwa syntaxis

反演结果中,Ⅵ区内最大主应力方向为北北西向或北西西向,与区域构造应力场主压应力方向不一致。该区内震源机制解样本数仅有3个,反演结果相对离散,可靠性有所降低。区域内最大地震为2019年4月24日的米林M6.3级地震,震源机制解走向/倾角/滑动角分别为202°/17°/20°(李国辉等,2020),基于该断层面参数反演得到的主压应力方向为N16°W(近南北向),由图3可以看出,该地震发生在喜马拉雅主边界断裂近似90°大拐弯地区,现今北北西向的最大主压应力方向可能反映了桑构造结西向俯冲南迦巴瓦构造结和北向推挤的动力过程(白玲等,2019)。

3.3 南迦巴瓦周边地区应力形因子与断层摩擦系数

Ⅰ—Ⅵ应力分区内应力形因子(R)反演结果如图7所示,整体上,南迦巴瓦周边地区应力形因子位于0.37～0.75之间,南迦巴瓦峰西侧(0.53～0.59)低于东侧(0.62～0.75)。杨帆等(2019)利用102个M4.7～6.5地震震源机制解反演得到南迦巴瓦构造结周边地区应力形因子为0.4～0.8,且显示该值自西北向东南呈增加趋势。文中反演得到的南迦巴瓦地区应力形因子结果与杨帆等(2019)研究成果基本一致。Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ区内优势应力形因子R值分别为0.59、0.53、0.62、0.70、0.75和0.37。依据林芝-通麦段地应力测量结果拟合得到的主应力随深度变化关系式(张重远等,2022),估算得到在1000～2000 m深度范围内应力形因子为0.44～0.65,平均为0.53,与文中Ⅰ区内反演结果0.59比较接近,一定程度上也表明地壳深浅部应力形因子具有较好的关联性。

图7 南迦巴瓦周边地区不同应力区应力形因子反演结果Fig.7 Stress shape ratio retrieved by the STRESSINVERSE program in different stress areas surrounding the Namcha Barwa syntaxis

将Ⅰ—Ⅵ区内每个震源机制解样本进行单独反演,还可以得到优势断层面滑动失稳时的临界摩擦系数μ值。Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ区内μ值分 别 为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个地应力测量数据计算得到重大铁路工程拉林段沿线1200 m深度范围应力积累水平表征参数μm值[(σ1-σ3)/(σ1+σ3)],其结果显示,μm值为0.23～0.30。张宁等(2022)基于青藏高原东南缘900组地应力实测数据,计算得到墨脱-昌都应力区内μm值为0.28。

由图2c应力摩尔圆可以看出,μm可以用公式(17)表达:

由公式(17)可以计算得到,川藏铁路拉林段沿线地壳浅层μ为0.24～0.31,墨脱-昌都应力区内μ为0.29,上述结果与文中估算得到的南迦巴瓦周边墨脱-波密地区(Ⅱ、Ⅲ区)μ最大值0.33基本一致,充分表明文中反演结果相对合理、可靠。

分析认为,南迦巴瓦周边地区主要发震断层相对低的摩擦系数可能与断裂带物质含较高的黑云母有关。以往研究表明,黑云母矿物在400℃剪切条件下的平均摩擦系数一般在0.25左右,最大不超过0.36(任凤文和何昌荣,2014;路珍和何昌荣,2014)。基础地质资料显示,南迦巴瓦周边地区片状分布含黑云母片岩、片麻岩等构造混杂岩带(潘桂棠等,2020;涂继耀等,2021)。

3.4 南迦巴瓦周边地区地应力大小

综合考虑南迦巴瓦周边地区交通廊道和清洁能源基地规划选址和建设需求,将重点对Ⅰ区、Ⅱ区、Ⅲ区和Ⅴ区内的地应力大小进行估算。

3.4.1 Ⅰ区地应力大小

Ⅰ区位于北东向东久-米林断裂带(北段)和北西向的西兴拉断裂(西段)交汇部位,第四纪以来,前者以左旋走滑活动为主,后者表现为右旋走滑运动方式,该区北部还展布北西西向嘉黎断裂带,晚第四纪以来表现出右旋走滑活动特征。此外,Ⅰ区内震源机制解类型除2017年11月18日的M6.9级地震为逆冲型外,大多中小地震以走滑型为主,与区域内主要发震断裂运动学特征相符合。鉴于此,选取走滑断层失稳临界应力条件和应力形因子计算公式进行估算。

将Ⅰ区内发震断层临界摩擦系数的最大值0.32代入公式(8),并将该区优势应力形因子R=0.59代入公式(12),孔隙水压力Pp近似等于静水压力P0。对于花岗岩、闪长岩和片麻岩,其平均天然密度取2.72 g/cm3,进而估算得到Ⅰ区内最大水平主应力值(σH)线性增加梯度为0.0337 MPa/m,最小水平主应力值(σh)线性增加梯度为0.0227 MPa/m,垂向应力线性增加梯度为0.0272 MPa/m(图8)。

图8 Ⅰ区主应力大小随深度变化特征Fig.8 Variation characteristics of principal stresses with depth in region Ⅰ

陈兴强(2022)基于通麦隧道15个钻孔地应力测量结果,拟合得到相对完整岩体内最大、最小水平主应力线性增加梯度分别为0.0336 MPa/m 和0.0222 MPa/m,与文中Ⅰ区估算结果吻合度高。依据张重远等(2022)拟合得到的林芝-通麦段主应力随深度变化关系式,估算该区1000 m深度最大、最小水平主应力值分别为32.95 MPa和23.02 MPa,由图8拟合关系式估算相同深度最大、最小水平主应力值分别为33.7 MPa和22.7 MPa,两者结果也很相近。总体上,文中估算得到的Ⅰ区内地应力大小结果是合理的。

3.4.2 Ⅱ区地应力大小

Ⅱ区内主要分布北西西向的嘉黎断裂带(中段),晚第四纪以来右旋走滑活动,区域内震源机制解均为走滑型。分别将μ最大值0.27代入公式(8)、R=0.59代入公式(12),估算得到Ⅱ区主应力大小随深度变化规律(图9),最大水平主应力值(σH)线性增加梯度为0.032 MPa/m,最小水平主应力值(σh)线性增加梯度为0.0229 MPa/m,垂向应力线性增加梯度为0.0272 MPa/m。位于区内的波密县嘎隆拉隧道钻孔应力解除法测量结果显示,在800 m埋深最大水平主应力值为25.3 MPa(孙晓光等,2011),由图9估算得到该深度最大水平主应力为25.6 MPa,与嘎隆拉隧道地应力实测结果基本相符。

图9 Ⅱ区主应力大小随深度变化特征Fig.9 Variation characteristics of principal stresses with depth in region Ⅱ

3.4.3 Ⅲ区地应力大小

Ⅲ区内主要分布北东向的墨脱断裂带,晚第四纪以来以左旋走滑运动为主,该区内震源机制解多表现为走滑型。分别将μ最大值0.33代入公式(8)、R=0.62代入公式(12),估算得到Ⅲ区主应力大小随深度变化规律(图10),最大水平主应力值(σH)线性增加梯度为0.0344 MPa/m,最小水平主应力值(σh)线性增加梯度为0.0227 MPa/m,垂向应力线性增加梯度为0.0272 MPa/m。根据张宁等(2022)拟合得到的墨脱-昌都应力区主应力随深度变化关系式,在1000 m深度最大、最小水平主应力值分别为29.89 MPa和17.34 MPa,由图10中关系式估算得到该深度最大、最小水平主应力值分别为33.4 MPa和22.7 MPa,该区水平主应力估算结果普遍高于墨脱-昌都区平均应力水平约4～5 MPa。分析认为,墨脱-昌都区覆盖面积近6×104km2,将整个区域内的地应力实测数据进行统一拟合,往往忽略地应力场的非均匀特征,部分高应力区难以得到体现。

图10 Ⅲ区主应力大小随深度变化特征Fig.10 Variation characteristics of principal stresses with depth in region Ⅲ

3.4.4 Ⅴ区地应力大小

Ⅴ区位于北东向墨脱断裂带的西南端,区域内震源机制解有正断型、逆冲型和走滑型、逆走滑型等,其中走滑型地震占比54%、逆冲型占比20%、正断型占比20%,综合上述考虑,该区仍然选取走滑断层滑动失稳临界条件进行地应力估算。分别将μ最大值0.32代入公式(8)、R=0.75代入公式(12),估算得到Ⅴ区主应力大小随深度变化规律(图11),最大水平主应力值(σH)线性增加梯度为0.0355 MPa/m,最小水平主应力值(σh)线性增加梯度为0.0236 MPa/m,垂向应力线性增加梯度为0.0272 MPa/m。由于该区内目前未见相关地应力测试数据公开,暂不能与实测结果进行可靠性比较。

图11 Ⅴ区主应力大小随深度变化特征Fig.11 Variation characteristics of principal stresses with depth in region Ⅴ

4 南迦巴瓦地区主要活动断裂构造稳定性探讨

南迦巴瓦地区位于印度板块与欧亚板块碰撞最强烈的前缘地带,即喜马拉雅东构造结雅鲁藏布江大拐弯附近,构造变形显著,新构造活动强烈,活动断裂发育,地震活动频发。1950年察隅8.6级地震后,该区内未发生M7.0级以上强震,多以中小地震为主(杨建亚等,2017),其中最大地震为2017年11月18日米林M6.9级地震,震中烈度达到Ⅷ度,Ⅵ度以上区域面积近1.3×104km2,涉及巴宜、米林、墨脱、波密及工布江达5个县(区)(唐方头等,2019)。

研究表明,1950年察隅8.6级和2017年米林6.9级地震导致周边区域应力场重新分布,对喜马拉雅东构造结周边活动断层产生较大影响,库仑应力变化达到数兆帕量级(尹凤玲等,2018;李振月等,2020)。南迦巴瓦构造结周缘主要活动断裂在现今构造应力背景下仍可能产生强震、错断等地质安全问题,影响高铁、水电等交通廊道地质安全风险防控(兰恒星等,2021;祁生文等,2022)。相关案例也屡见不鲜,如2022年1月8日青藏高原东北缘门源县M6.6级地震,沿冷龙岭断裂产生明显地表破裂带,地表最大同震位移达3.1 m,发震断裂的左旋走滑剪切作用直接导致兰新高铁大梁隧道、硫磺沟大桥等设施严重损坏(张威等,2022;韩帅等,2022)。

文章聚焦活动断裂构造稳定性问题,采用基于摩尔-库仑破裂准则的断层滑动失稳概率分析方法(Snee and Zoback, 2018),将上文地应力估算结果作为边界条件,定量评价南迦巴瓦周边Ⅰ区、Ⅱ区、Ⅲ区和Ⅴ区内主要活动断裂滑动失稳风险,以期为交通廊道和清洁能源基地活动断裂地质安全风险评价提供参考。

为便于计算,根据不同断裂走向,将4个应力区内的主要活动断裂简化为不同长度的直线段(图12)。参考地质、地球物理等资料,合理确定各分段断层的走向、倾角、长度等参数,同时考虑不同参考资料之间的差异性,分别将各分段断层的走向和倾角设定为±5°和±10°误差范围(表2)。如上所述,4个区域内反演得到的断层摩擦系数与其他研究成果相比,最大残差为0.08,故将各分段断层摩擦系数的误差范围设定为±0.08。基于摩尔-库仑准则的断层滑动失稳概率分析方法,可详见作者团队已发表论文(朱思雨等,2022;Zhu et al., 2022; Feng et al., 2022; Fan et al., 2022)。

表2 南迦巴瓦周边地区主要活动断裂带属性参数Table 2 Characteristic parameters of the main active fault belts surrounding the Namcha Barwa syntaxis

图12 南迦巴瓦周边地区主要活动断裂带简化分段特征Fig.12 Simplified segments of the main active fault belts surrounding the Namcha Barwa syntaxis

采用美国斯坦福大学开发的断层滑动失稳概率计算程序FSP1.0(Walsh Ⅲ and Zoback, 2016),定量估算Ⅰ区、Ⅱ区、Ⅲ区和Ⅴ区内主要活动断裂当前滑动失稳概率(图13)。

图13 南迦巴瓦周边地区主要活动断裂滑动失稳概率Fig.13 Fault slip probability in the main active fault belts surrounding the Namcha Barwa syntaxis

(1)Ⅰ区内,嘉黎断裂带西北段(F1-4)、西兴拉断裂带西支(F7-1、F7-2)和东久-米林断裂带东北支(F9-3)在当前地应力作用下滑动失稳概率很低,一般不超过10%,断层相对稳定,而东久-米林断裂带东北段(F9-2、F9-1)滑动失稳概率较高,分别为25%和39%,需重点关注F9-2段未来地震活动危险性。

(2)Ⅱ区内,嘉黎断裂带中段(F1-1、F1-2和F1-3)滑动失稳概率较高,分别为24%、38%和39%,地震活动危险性较高,而墨脱断裂带的东北端(F6-9)滑动失稳概率很低,小于5%,断层活动相对稳定。

(3)Ⅲ区内,墨脱断裂带中段(F6-4、F6-5、F6-6、F6-7、F6-8)整体滑动失稳概率偏高,位于16%～45%之间,其中,背崩乡至帮辛乡段断层滑动失稳概率普遍大于30%,最高达45%(F6-8),需要重点加强对墨脱断裂带中段地震活动性的持续跟踪与动态监测。阿帕龙断裂带滑动失稳概率较高,为36%～39%,地震活动危险性高,已有GPS观测数据显示,阿帕龙断裂带可能处于闭锁阶段,为强震危险区。

(4)Ⅴ区内,墨脱断裂带南段(F6-1、F6-2、F6-3)当前滑动失稳概率很低,均小于8%,相比之下,喜马拉雅主边界断裂北段(F4-1)滑动失稳概率很高,最高可达50%,2019年4月24日M6.3级地震就发生在F4-1断层南侧约20 km位置,该区未来仍是中强地震频发区。

此外,从交通廊道地质安全风险防控角度来看,交通廊道波密至林芝段在古乡附近转至与嘉黎断裂带(F1-4)平行,由于嘉黎断裂带(F1-4)当前滑动失稳概率很低,断层活动稳定,可能不会对工程建设和安全运行产生严重影响,比较而言,东久-米林断裂东北段(F9-2)具有较高的滑动失稳概率,后期需深入研究该条断裂发生强震时对交通廊道产生的工程错断效应。另外,对于区域内水电工程规划选址而言,则需要重点评估雅鲁藏布江下游墨脱断裂带的地震活动危险性及其可能产生的地质安全问题。

5 结论与建议

(1)南迦巴瓦周边地区现今构造应力场最大主应力方向为北东—北北东向,反映了印度板块北向俯冲青藏高原的动力背景。

(2)南迦巴瓦周边地区地应力大小表现出非均匀特征,最大、最小水平主应力线性增加梯度在西北部的“易贡-鲁朗-派镇-古乡”区为0.0337 MPa/m、0.0227 MPa/m,东北部的“古乡-帮辛-那母隆巴-松宗-波密”区为0.032 MPa/m、0.0229 MPa/m,东南部的“帮辛-都登-阿帕龙-布宗”区为0.0344 MPa/m、0.0227 MPa/m,西南部的“都登-米林-梅楚卡-格刀”区为0.0355 MPa/m、0.0236 MPa/m。

(3)在当前地应力环境下,南迦巴瓦周边地区主要断裂活动危险性也具有差异性,东久-米林断裂带东北段(派镇至古乡)滑动失稳概率为25%～39%,嘉黎断裂带中段(古乡至松宗)为24%～39%,墨脱断裂带东北段(背崩至帮辛)为30%～45%,阿帕龙断裂带西北段(都登至阿帕龙)为30%～45%,以上区段均具有较高的地震活动性。水电工程规划选址阶段,需重点关注墨脱断裂带未来大概率发生强震并可能引起工程错断等地质安全问题。

(4)依据断层滑动失稳临界应力条件,结合基于震源机制解数据的应力形因子、断层摩擦系数等参数反演,可估算地应力场的优势方向和绝对大小,估算与实测结果相似性较好,一定程度上可为未开展原位地应力测量的区域提供补充。估算过程中,设定的区域面积、震源机制解样本数量、地震震级与类型等因素一定程度上会影响应力形因子和断层摩擦系数结果,因此,文中估算得到的南迦巴瓦周边地区地应力场,仍需要结合区域内原位地应力实测数据数据进行检验和完善。