2021年3月19日那曲M6.1地震震源区应力应变特征及其动力学意义

李玉江，邵志刚，胡幸平，石富强，刘皓晴，陈连旺

1 应急管理部国家自然灾害防治研究院，北京 100085 2 中国地震局地震预测研究所，北京 100036 3 陕西省地震局，西安 710068

0 引言

据中国地震台网中心测定，2021年3月19日在西藏那曲比如县发生M6.1地震，震中位于92.74°E，31.94°N，震源深度10 km. 地震发生在印度板块与欧亚板块碰撞下的强烈变形地带，区域内活动断裂发育，历史强震活跃(图1)，曾发生1951年崩错M8.0和1952年当雄北M7.5两次强震(吴章明和邓起东, 1989; 吴章明和曹忠权, 1991). 那曲地震发生后，不同研究机构利用全球或区域地震台网数据反演给出此次地震的震源机制解(走向/倾角/滑动角)(表1).多数研究结果显示，该地震是发生在近北东走向断层上的一次以正断活动为主兼左旋走滑运动性质的地震.

地震的活动性质与区域构造应力场结构密切相关(Zoback, 1992; Hu et al., 2017; 李泽潇等, 2020).程惠红等(2014)利用有限元模型给出新疆于田地区的主应力场特征，发现应力结构的差异较好地解释两次不同类型于田地震的发震机制.胡幸平等(2021)利用长宁地区精细的地壳应力场结构，较好地解释区域内中强震震源机制的空间差异. 世界应力图项目综合各类地应力数据，给出全球范围内的应力场结构(Heidbach et al., 2018)，为理解不同构造部位强震发生的动力学过程提供力学依据.那么，那曲地震震中区构造应力场如何分布以及它是如何控制那曲地震以正断运动为主的活动性质？亟待深入研究.

表1 那曲地震震源机制解Table 1 Focal mechanism solution of Nakchu earthquake

地震的发生除受到区域构造应力场控制外(Heidbach and Ben-Avraham, 2007)，还会受到周边强震的影响(Harris, 1998; 石耀霖和曹建玲, 2010; He et al., 2011; Wang et al., 2014).从历史强震活动性来看，在那曲地震周边曾发生过1411年当雄M72/3、1934年申扎M7.0、1950年察隅M8.7等强震，包括在地理空间上距离较近的1951年崩错M8.0和1952年当雄北M7.5两次强震.这些历史强震的发生会对那曲地震产生什么样的影响，值得深入分析.库仑破裂应力作为探讨强震间相互作用机制的一种手段(Stein, 1999; He et al., 2011; Shan et al., 2013)，为认识强震的发震机理提供理论依据.

为认识上述问题，本文利用“中国大陆地壳应力环境基础数据库”中的应力数据，进行数据平滑计算分析，给出震中周边构造应力场背景特征.基于球面小波多尺度方法，解算GPS数据获得区域应变率场特征.考虑周缘历史强震的同震位错与震后黏弹性松弛效应，分析那曲地震震源处应力场演化与应力变化.综合应力应变场结果，探讨那曲地震发生的力学环境与发震机理.

1 计算方法与模型建立

1.1 计算方法

本文从“中国大陆地壳应力环境基础数据库”中搜集震中周边区域的各类地应力数据，按照世界应力图的标准对这些应力数据进行质量评级(Heidbach et al., 2018)，从中挑选出质量可靠的A-C级应力数据，并依据数据质量赋予不同的质量权重系数(A级1.0，B级0.75，C级0.5)，进行应力方向的平滑计算分析(Hansen and Mount, 1990; Heidbach et al, 2010)，给出区域应力场的分布变化规律，具体平滑分析方法见Hu 等 (2017).

基于高斯差分的球面小波多尺度方法(Tape et al., 2009; Li et al., 2019)，解算震中周边不同空间尺度的应变率场.该方法的优点是，可以通过确定小波基函数的阈值，实现自动根据GPS观测点密度差异在不同区域使用不同的内插尺度，获得反映不同空间尺度下的地壳形变特征(苏小宁等, 2016).高斯差分球面小波函数描述为，在半径为1的单位球，球面上任何一点x处的小波函数表达式为

(1)

式中，γ为球面坐标下，观测点位矢量x与球面小波中心极之间的夹角，取值范围为0≤γ≤180°；a=2-q，q表示尺度，值越大，尺度越小，q=0,1,2…；α>1，用于调节函数形状.

基于库仑破裂准则的库仑应力变化的表达式为(King et al., 1994)

Δσf=Δτ+μ′Δσn,

(2)

式中，Δτ为断层面上剪应力变化，与断层滑动方向一致时为正；Δσn为断层面上的正应力变化，定义张性为正；μ′为断层面介质的有效摩擦系数.如果Δσf>0，则有利于后续地震的提前发生.

参考前人的做法(King et al., 1994)，本研究中有效摩擦系数μ′取0.4.使用弹性位错模型计算同震库仑应力变化，同时考虑到Kelvin体和Maxwell体在表征瞬态变形和长期稳态变形方面的缺陷，采用能够综合协调瞬态变形和长期稳态变形的Burgers体来模拟震后黏弹性松弛效应(Shao et al., 2016; Li et al., 2020a).Burgers体的本构关系为

(3)

其中，σ和ε分别表示应力和应变；k1和k2分别为Kelvin体和Maxwell体的弹性模量；η1和η2分别表示瞬态和稳态粘滞系数.

采用基于分层黏弹性地球模型的PSGRN/PSCMP程序(Wang et al., 2006)，计算强震的同震位错与震后黏弹性松弛效应引起的应力变化.根据地震震源深度的结果(表1)，库仑应力变化计算深度为10 km.

1.2 介质模型

参考青藏高原地壳上地幔速度结构(滕吉文等, 2012)，以及2015年尼泊尔地震震源区下地壳和地幔岩石圈流变性质的研究成果(Zhao et al., 2017; Wang and Fialko, 2015)，确定介质模型的速度结构与中-下地壳、上地幔黏滞系数.选择Burgers体模型模拟中-下地壳、上地幔的黏弹性流变特性(Pollitz et al., 2001; Wang et al., 2012).其中，上地壳为弹性介质，深度为20 km；中、下地壳和上地幔为黏弹介质. 具体介质参数见图2.

1.3 震源模型

依据地震地质、震源破裂过程反演等方面的研究结果(Armijo et al., 1989; 吴章明和邓起东, 1989; 吴章明等, 1990; 吴章明和曹忠权, 1991; Ambraseys and Douglas, 2004; Bilham, 2004;吴中海等, 2008; 张勇等, 2010; 李保昆等, 2014)，确定历史强震的震源模型参数.其中，1947年朗县M7.7地震由于缺少直接的观测与反演结果，其破裂长度和滑动量依据经验公式给出(Wells and Coppersmith, 1994)；对于1950年察隅地震，基于余震分布和滑坡资料的研究认为，此次地震存在两个明显的破裂面(Coudurier-Curveur et al., 2020)，本文依据其结果确定地震的震源机制参数.同时在“讨论”中对不同震源模型参数下的结果开展对比分析(Ben-Menahem et al., 1974; 李保昆等, 2015).震源模型参数具体见表2.

图2 黏弹性介质分层模型参数VP表示P波速度，VS表示S波速度，ρ表示密度；η1、η2分别表示瞬态和稳态粘滞系数.Fig.2 The parameters of the stratified viscoelastic relaxation modelVP indicates the velocity of P wave, VS the velocity of S wave, and ρ the density; η1 and η2 represent the transient and steady-state viscosity coefficients, respectively.

表2 历史强震震源模型参数Table 2 Parameters of source model for strong historic earthquakes

2 结果

2.1 区域应力应变环境

图3为本文从“中国大陆地壳应力环境基础数据库”中搜集出的震中及周边的各类应力数据(图3a)，以及考虑数据质量及距离加权的水平最大主应力方向平滑结果(图3b).可以看出，由于地处青藏高原腹地，研究区内应力数据主要为震源机制解，这些震源机制解的水平最大主压应力方向具有较好的一致性，在研究区内基本都呈现为北北东-南南西向.同时，利用搜集到的51个M≥4.3历史地震的震源机制解结果，采用计算综合震源机制解的方法(许忠淮等, 1983)，对该地区的构造应力场进行了反演，结果也表明区域应力场的最大主压应力轴(σ1)方位为北北东-南南西向(图3c).在此应力控制下，对于走向为近北东-南西的发震断层，具备了发生以正断为主兼具左旋走滑运动的动力学特征.结合不同研究机构给出的那曲地震的震源机制解(表1)，可以看出，发震节面走向与最大主压应力轴方位近乎平行，以正断性质为主的那曲地震就是在这样的构造应力环境下发生.

同时，选取1991—2016年GPS观测资料(Wang and Shen, 2020)，基于球面小波多尺度方法，选取q=2～7尺度因子，解算得到震中周边最大剪应变率、面应变率和主应变率方位(图4).结果显示，在2～7尺度因子下反演的速度场结果与GPS观测结果较为吻合(图4a)，两者之间的速度场残差绝大多数<1.0 mm·a-1(图4b).震中区最大主压应变率的方位为北北东向，与最大主压应力轴方位基本一致(图4c)，均和震源机制解显示的发震节面近乎平行.同时，地震发生在最大剪应变率的梯度带上.

2.2 震源处应力演化

基于强震同震位错模型和岩石圈分层模型，利用PSGRN/PSCMP程序(Wang et al., 2006)计算历史强震的同震位错和震后黏弹性松弛效应引起的应力张量变化，并将应力张量变化投影到那曲地震的发震断层面上.由于地震的发震断层面参数存在差异，本文选取5组以正断性质为主的震源机制解，对发震断层面分别进行应力张量投影，给出μ′=0.4时震源处库仑应力演化和应力变化.

图3 那曲地震震中周边构造应力环境(a) 震中周边应力数据； (b) 平滑给出的水平最大主应力方向； (c) 基于震中周边历史地震震源机制解的区域应力场反演结果.Fig.3 Tectonic stress environment around the Nakchu earthquake epicenterPanels (a), (b) and (c) represent stress data around the epicenter, orientation of the horizontal maximum principal stress after smoothing, and the inversed regional stress field based on historical focal mechanism solution in surrounding areas, respectively.

计算结果显示，使用5组震源机制参数分别进行应力张量投影，所计算出的库仑应力变化尽管存在微弱差异，但应力变化的极性是稳定的(图5).以中国地震局地球物理研究所(IG)给出的震源机制参数为例(下同)，1411年当雄南M72/3以来9次历史强震的发生造成那曲地震震源处同震和累积库仑应力变化分别为1.63×104Pa和9.06×104Pa.其中，1950年察隅M8.7地震的发生造成那曲地震震源处明显的库仑应力增加，同震和震后应力变化分别为1.91×104Pa和8.27×104Pa.分析地震引起的断层面上正应力和剪应力变化可以看出，察隅地震的发生造成那曲地震震源处正应力发生显著的张性变化、剪应力发生与断层滑动方向一致的变化(图6a，6b)，进而引起库仑应力的增加(图6c).同时，来自中-下地壳、上地幔的震后黏弹性松弛效应亦造成震源处正应力、剪应力和库仑应力的显著增加(图6e—6f).相对于察隅地震引起的较为明显的库仑应力变化，其它历史强震的影响相对较小，各历史强震引起的那曲地震震源处同震和震后应力变化具体见表3.

表3 历史强震引起的那曲地震震源处库仑应力变化(以IG震源机制参数为例)Table 3 Coulomb stress change at the hypocenter of the Nakchu earthquake caused by historic strong earthquakes

图4 那曲地震震中周边应变率场特征(a) 速度场实测值(箭头)与模拟值(背景等值线)； (b) 速度场残差； (c) 最大剪应变率与主应变率(箭头).Fig.4 Characteristics of the strain rates around the epicenter of the Nakchu earthquake(a) GPS velocity field (black arrows), the background contours show the magnitude of the modeled horizontal velocity; (b) the residual of velocities, and (c) the maximum shear strain rate and principal strain rate (black arrows).

图5 不同震源机制参数下那曲地震震源处库仑应力演化(a) 同震库仑应力变化； (b) 累积库仑应力变化 (同震+震后).Fig.5 Coulomb stress evolution at the hypocenter of the Nakchu earthquake with different focal mechanism parametersPanel (a) denotes the co-seismic Coulomb stress change, while (b) the combined Coulomb stress change (co- and post-seismic).

图6 察隅地震引起的同震与震后应力变化(a)—(c)分别表示同震正应力、剪应力和库仑应力变化；(d)—(f)分别表示震后正应力、剪应力和库仑应力变化.Fig.6 Co- and post-seismic stress changes associated with the Chayu earthquakePanels (a)—(c) denote the co-seismic normal stress change, shear stress change, and the Coulomb stress change, respectively. Panels (d)—(f) denote the post-seismic normal stress change, shear stress change, and the Coulomb stress change, respectively.

3 讨论

强震的孕育、发生和活动性质不仅受区域构造应力应变场环境的控制(Heidbach et al., 2010)，也与断裂带应力变化密切相关(Heidbach and Ben-Avraham, 2007; Shan et al., 2013)，而断裂带应力变化的计算又受到震源破裂模型以及摩擦系数等参数的影响(Xiong et al., 2010; 徐杜远等, 2020).基于InSAR资料反演给出的那曲地震同震破裂模型的研究结果显示，那曲地震的震源机制参数为237°/69°/-70°(走向/倾角/滑动角)(李永生，个人通讯)，与CGMT、GFZ、中国地震局地球物理研究所等机构给出的结果基本一致，均反映出该地震是一次以正断活动为主兼左旋走滑运动性质的地震.

3.1 震源参数差异对应力变化的影响

1950年察隅地震发生在喜马拉雅东构造结地区，是有地震记录以来中国大陆内部最强的地震(Coudurier-Curveur et al., 2020).受当时监测技术和条件等限制，很难获取准确的断层面解，造成此次地震的活动机制存在右旋走滑型(Ben-Menahem et al., 1974)和低角度逆冲型(Chen and Molnar, 1977)两种不同认识.基于P波初动、余震重定位和地质特征的研究结果显示(李保昆等, 2015)，察隅地震的震源机制解与Ben-Menahem等(1974)的结果基本一致.而基于余震分布和滑坡资料的最新研究认为，察隅地震是一次以逆冲为主兼具走滑性质的地震(Coudurier-Curveur et al., 2020).为此，本文使用走滑型和逆冲兼走滑型两组震源模型分别计算应力变化，进一步讨论震源模型参数选取对应力变化结果的影响.结果显示，在察隅地震为右旋走滑型时，地震的同震与震后黏弹性松弛效应则分别造成那曲地震震源处应力减小1.03×104Pa和7.04×104Pa(图7)，察隅地震引起的卸载效应起到了延缓那曲地震发生的作用；而在察隅地震为低角度逆冲兼走滑型时，地震的同震与震后黏弹性松弛效应则分别造成那曲地震震源处应力增加1.91×104Pa和8.27×104Pa，察隅地震引起的加载效应造成了那曲地震的提前发生.

从同震与震后位移场特征来看(图8)，在察隅地震为低角度逆冲兼走滑型机制时(Coudurier-Curveur et al., 2020)，同震与震后附加位移场方向均为南东方向，与那曲地震发震节面的方向近垂直.这种附加位移场沿南东方向的逐渐增大，意味着那曲地震震源区受到南东优势方向上的拉张作用，有利于那曲地震的发震断层发生正断性质的活动，与库仑应力增加所反映的特征一致(图7).而在察隅地震为右旋走滑型时(Ben-Menahem et al., 1974)，产生了近西向的同震位移场，且由东南向西北逐渐减小(图8a)，意味着那曲震源区受到东西优势方向上的挤压作用，对于走向为北东的断层而言，有利于发生右旋走滑运动，与那曲地震左旋走滑的发震机制相反，与同震库仑应力减小所反映的特征一致.同时，震后黏弹性松弛效应引起的位移场方向为北西向(图8b)，虽然也与那曲地震发震节面的方向近垂直，但是这种位移场在北西方向的减小意味着那曲地震震源区受到北西优势方向上的挤压作用，不利于那曲地震的发震断层发生正断性质的活动，与库仑应力减小所反映的特征一致(图7).

图7 察隅地震发生引起的那曲地震震源处库仑应力演化Fig.7 Coulomb stress evolution at the hypocenter of the Nakchu earthquake caused by the Chayu earthquake

另外，对于1411年当雄南等历史强震，由于发震时间较早使得难以合理地回溯其破裂的基本参数.为深入认识震源模型参数的选取对应力计算结果的影响，我们基于震级与破裂长度、滑动量之间的经验公式(Wells and Coppersmith, 1994)，给出历史强震的破裂长度和滑动量(表4)，并计算在新的震源参数下那曲地震震中的应力演化(图9).可以看出，两组参数下震源区应力演化的计算结果是稳定的，但存在个别地震引起的应力变化在总的应力变化中的权重显著增加或减小的现象.比如，在新的震源参数约束下，1951年崩错地震同震与震后效应引起的应力变化相对图5中的结果明显要大.因此，在计算断裂带应力变化和探讨强震间相互作用机理时，应尽可能使用与实际接近的震源破裂模型.

表4 基于经验公式的历史强震震源模型参数

图8 基于两组不同震源机制解的同震与震后位移场其中，红色和蓝色箭头分别表示察隅地震为右旋走滑型机制和低角度逆冲兼走滑型机制下产生的位移场.(a) 同震位移场；(b) 震后黏弹性松弛效应引起的位移场.Fig.8 Co- and post-seismic displacement fields based on two different focal mechanism solutionsIn which, the red and blue arrows indicate the displacement fields generated by the dextral strike-slip and the lower-angle thrust with strike-slip motion of the Chayu earthquake, respectively. (a) Co-seismic displacement; (b) Displacement induced by post-seismic viscoelastic relaxation effect.

图9 两组震源模型下那曲地震震源处库仑应力演化其中，模型1中的参数见表2；模型2中的参数见表4.Fig.9 Coulomb stress evolution at the hypocenter of the Nakchu earthquake with two different source modelsIn which,source parameters in model 1 and model 2 are shown in Table 2 and Table 4, respectively.

3.2 摩擦系数选取对应力变化的影响

库仑应力变化是研究强震间相互作用的重要参考，但其受断层有效摩擦系数取值的影响(Xiong et al., 2010).由公式(2)可知，在库仑应力计算中，有效摩擦系数影响着正应力与剪应力的权重.为分析有效摩擦系数选取对库仑应力计算结果的影响，以中国地震局地球物理研究所给出的震源机制参数(231°/55°/-65°)为例，选取不同的有效摩擦系数(μ′=0.0，0.2，0.4，0.6，0.8)，进一步探讨有效摩擦系数选取对震源处库仑应力变化计算结果的影响.由图10可知，对于那曲地震而言，有效摩擦系数的取值会改变震源处库仑应力变化量的大小，但应力变化的极性一致.在1950年察隅地震发生前，断层面库仑应力变化随着摩擦系数的增加而逐渐减小，意味着此次地震前历史强震发生所产生的断层面累积正应力变化为负值，断层区处于压性应力变化状态.而1950年察隅地震的发生，造成发震断层面较强的张性正应力变化.在此前提下，有效摩擦系数越大，有效正应力权重越大，库仑应力变化量增加的越多.

图10 不同有效摩擦系数下震源处累积库仑应力变化Fig.10 Cumulative Coulomb stress change at the hypocenter with different effective coefficients of friction

3.3 那曲地震发震机理探讨

在印度板块与欧亚板块的持续碰撞作用下，青藏高原内部发生强烈的变形，并伴随强震的发生(Li et al., 2018; Liu-Zeng et al., 2020; Wang and Shen, 2020).基于实测地应力数据的分析结果显示青藏高原内部最大水平主应力的优势方向为NE-SW向(杨树新等, 2012)，与历史强震震源机制解反演的结果(Bai et al., 2017)和本文平滑计算后的主应力方向结果基本一致(图3).现今GPS观测结果显示，相对于稳定的欧亚板块，那曲地震震源区周边速度场从NNE向逐渐向近EW向发生顺时针旋转(Gan et al., 2007; Zheng et al., 2017)，最大主压应变率的方向也同时从NNE向转变为近EW向(Allmendinger et al., 2007; Wang and Shen, 2020)，与本文利用1991—2016年较长时间尺度GPS数据解算的结果一致(图4).2015年尼泊尔M7.8地震的发生，造成震源区周边应力变化为-10～100 Pa(Zha and Dai, 2017)；2017米林M6.9的发生对震源区的影响可以忽略(尹凤玲等, 2018)，这种较小的应力变化与本文计算的结果基本一致(表3).从震源处应力演化的结果来看(图5)，1411年当雄南M72/3以来9次历史强震的发生造成那曲地震震源处累积库仑应力变化9.06×104Pa.

综合上述分析，初步认为2021年3月19日那曲M6.1地震是在印度板块与欧亚板块持续挤压的动力学背景下发生的，区域历史强震活动起到明显的促进作用.

4 结论

本文基于那曲震中周边的应力数据和GPS形变观测数据，分析区域构造应力场和应变率场特征.采用黏弹性Burgers流变模型计算1411年当雄地震以来9次历史强震引起的那曲地震震源处同震、震后库仑应力演化和应力变化.结果显示：

(1) 在印度板块与欧亚板块持续碰撞的动力学环境下，那曲地震震中区域最大主压应力轴的方位为北北东向，与发震断层近北东-南西的走向近乎平行.该构造应力环境控制着那曲地震以正断为主兼具左旋走滑性质的运动.

(2) 震中区最大主压应变率的方位为北北东向，和最大主压应力方位基本一致，均与震源机制解显示的发震节面近乎平行.同时，地震发生在最大剪应变率的梯度带上.

(3) 1411年当雄南M72/3以来9次历史强震的发生造成那曲地震震源处累积库仑应力变化9.06×104Pa.其中，1950年察隅M8.7地震的发生造成震源处同震和震后应力变化分别为1.91×104Pa和8.27×104Pa，与同震和震后黏弹性松弛效应引起的附加位移场所指示的震源处应力加载效应的结果一致.

(4) 综合分析区域构造应力应变场与应力变化的研究结果，初步认为2021年3月19日那曲M6.1地震是在印度板块对青藏高原北东向挤压的动力学背景下发生的，区域历史强震促进那曲地震的提前发生.该研究为认识那曲地震的发震背景和动力学机理提供理论依据.

致谢在成文过程中，作者与崔效锋研究员、孟令媛研究员、赵静旸博士、李永生博士的有益交流和讨论.图件使用GMT软件绘制(Wessel et al., 2013)，在此一并表示感谢.