2017年11月18日西藏米林 M6.9地震对后续地震的静态库伦应力的影响
李振月3), 万永革1,2),*, 靳志同1,2), 杨帆4), 胡晓辉3), 李泽潇1)
1)防灾科技学院, 三河 065201
2)河北省地震动力学重点实验室, 三河 065201
3)中国科学技术大学, 合肥 230026
4)云南省地震局, 昆明 650224
*通讯作者: 万永革, 男, 1967年生, 研究员, 主要从事构造应力场、 地震应力触发等方面研究工作, E-mail: wanyg217217@vip.sina.com.cn

〔作者简介〕 李振月, 男, 1994年生, 2020年于中国地震局地球物理研究所获地球探测与信息技术专业硕士学位, 现为中国科学技术大学地球物理学专业在读博士生, 主要从事构造应力场、 应力触发方面的研究工作, E-mail: 952934956@qq.com

摘要

基于2017年11月18日西藏米林 M6.9地震的破裂模型, 分别计算了2种接收断层选择方案下, 不同深度的静态库伦破裂应力变化的空间分布。 一种方案认为接收断层的参数和发震断层一致; 另一种方案假设断层的方位随机分布, 选择在主震同震应力场作用下最容易发生错动的断层面作为接收断层面。 根据以上2种方案下静态库伦破裂应力变化的结果, 讨论主震对震后短期内余震以及2019年4月24日西藏墨脱 M6.3地震的静态库伦应力的影响。 结果表明: 1)当接收断层和发震断层一致时, 各深度处于静态库伦破裂应力变化正值区的余震比例较小, 正值区余震的震源机制和主震相似, 分析认为是主震继续破裂所导致。 2)各深度大部分余震处于静态库伦破裂应力变化的负值区, 考虑可能是由于这些余震和主震的震源机制差异较大所致。 通过选择最容易发生错动的断层面作为接收断层面计算静态库伦破裂应力变化的结果可发现, 不同深度的余震均处于静态库伦破裂应力变化大于触发阈值0.01MPa的范围内, 说明所有余震均有被触发的可能性。 推测处于静态库伦破裂应力变化负值区的余震可能发生在主震剧烈破裂导致的破碎带内, 从而导致其震源机制解与主震存在较大差异。 破碎带的存在会显著降低震源区的弹性常数值, 文中的结果也说明考虑震源区内外弹性常数的差异对准确估计震源区内地震之间的静态库伦应力影响是有益的。 3)根据不同机构和作者利用不同资料和方法给出的墨脱地震的震源机制解, 计算得到一个与这些震源机制差别最小的中心震源机制解, 进而据此中心震源机制解定量计算了米林地震对墨脱地震的静态库伦应力影响。 结果显示, 米林地震在墨脱地震中心震源机制解的2个节面上产生的静态库伦破裂应力变化的量值都很小, 说明墨脱地震是一次独立于米林地震的事件。

关键词: 米林 M6.9地震; 破裂模型; 接收断层; 静态库伦破裂应力变化; 墨脱 M6.3地震
中图分类号:P315.2 文献标志码:A 文章编号:0253-4967(2020)05-1091-18
THE STATIC COULOMB STRESS INFLUENCE OF THE MAINLING M6.9 EARTHQUAKE IN TIBET ON NOVEMBER 18, 2017 TO THE SUBSEQUENT EARTHQUAKES
LI Zhen-yue3), WAN Yong-ge1,2), JIN Zhi-tong1,2), YANG Fan4), HU Xiao-hui3), LI Ze-xiao1)
1)Institute of Disaster Prevention, Sanhe, Hebei 065201, China
2)Hebei Key Laboratory of Earthquake Dynamics, Sanhe, Hebei 065201, China
3)University of Science and Technology of China, Hefei 230026, China
4)Yunnan Earthquake Agency, Kunming 650224, China
Abstract

Based on the rupture model of Mainling M6.9 earthquake in Tibet on November 18, 2017, the spatial distribution of static Coulomb failure stress change at different depths are calculated respectively according to two different receiving fault selection schemes. The one scheme is that we set the parameters of receiving fault at different position to be consistent with the main shock; The other scheme is on the assumption that fault's orientation is randomly distributed under the ground, and we select the receiving fault which is most prone to slide under the influence of coseismic stress field produced by main shock. Therefore, the geometrical orientation of receiving fault will vary with space. According to the above two results of static Coulomb failure stress change, we discussed the static Coulomb stress influence produced by the main shock to short-term aftershocks and the Medog M6.3 earthquake in Tibet on April 24, 2019, respectively. The result shows that: 1)When the parameters of receiving fault are same with the main shock, the proportion of aftershocks in the positive zone of static Coulomb failure stress change is small at each depth. The focal mechanisms of aftershocks in the positive zone of static coulomb fracture stress are deemed similar to the main shock. We thought that they are motivated by the continuous rupture of the main shock. 2)Most of the aftershocks are in the negative zone of static Coulomb failure stress change at each depth. We inferred that this phenomenon which may be on account of the focal mechanisms of these aftershocks is quite different with the main shock. From the result of receiving fault to choose the most prone to slide under the coseismic stress field produced by main shock, we can clearly see that all the aftershocks are within the zone of static Coulomb failure stress change greater than the trigger threshold of 0.01MPa at different depths. It indicates that all the aftershocks are likely to be triggered. It was speculated that the aftershocks in the negative zone of static Coulomb failure stress change occurred in the crushed zone caused by violent rupture of the main shock. There are countless cracks in the crushed zone, and the orientation of these cracks is abundant. Perhaps, because most aftershocks occurred on these various cracks, their focal mechanisms are quite different from the main shock. The value of elastic constants will be reduced significantly in the crushed zone. All the results in this paper also indicate that considering the elastic constants difference between in and out of the source region is beneficial to accurately estimate the static Coulomb stress influence between earthquakes in the source region. 3)Different institutes and authors used different data and methods to get several different focal mechanisms of the Medog earthquake. According to these results, we calculated a central focal mechanism solution, which has a minimum difference with these focal mechanisms. On the basis of this central focal mechanism solution, the static Coulomb stress influence of the Mainling earthquake to the Medog earthquake is calculated quantitatively. Result indicates that the magnitude of static Coulomb failure stress change generated by the Mainling earthquake is quite small on both two nodal planes of the central focal mechanism solution of the Medog earthquake, this means that the Medog earthquake is independent of the Mainling earthquake.

Keyword: Mainling M6.9 earthquake; rupture model; receiving fault; static Coulomb failure stress change; Medog M6.3 earthquake
0 引言

2017年11月18日西藏米林县M6.9地震(29.75° N, 95.02° E, 简称米林地震)发生后不久, 2019年4月24日在距离较近的墨脱县发生了M6.3地震(28.4° N, 94.61° E, 简称墨脱地震)(中国地震台网), 这2次地震均发生在青藏高原SE部旋转变形最为强烈的南迦巴瓦构造结(图 2)。 米林地震发生后, 中国地震台网中心发布的《2017年11月18日西藏米林6.9级地震图集和一些学者(白玲等, 2017; 吴宝峰, 2017; 刘云华等, 2018)对其发震方式的研究结果均表明, 米林地震是一次以逆冲为主的事件。 印度板块对欧亚板块NE向的碰撞挤压, 导致东喜马拉雅构造结的逆冲推覆和青藏高原物质的侧向扩张(Shen et al., 2017; 李煜航, 2018)是米林地震的主要发震背景(白玲等, 2017)。 墨脱地震后, 中国地震台网中心发布的《2019年4月24日西藏墨脱6.3级地震图集》和国外网站发布的该地震的震源机制结果(表1)显示, 墨脱地震也是一次以逆冲为主的事件。 同时, 米林地震与墨脱地震的震中位置较近, 墨脱地震是否为米林地震的余震?米林地震对墨脱地震的发生是否有促进或延缓作用?为了解答以上问题, 需要研究米林地震对墨脱地震断层面静态库伦应力的影响。 在米林地震与墨脱地震的关系未揭晓之前, 下文所述的米林地震的余震序列中均不包含墨脱地震。

表1 西藏墨脱M6.3地震的不同震源机制解及将其作为初始解得到的中心震源机制解 Table1 Focal mechanisms of the Medog M6.3 earthquake in Tibet and the corresponding central focal mechanisms derived as initial solution

使用描述米林地震破裂分布的模型可以更准确地评估该地震对周围断层的静态库伦应力的影响。 基于米林地震的破裂模型, 本文首先将分析米林地震对余震静态库伦应力的影响, 之后再研究米林地震对墨脱地震的影响。 米林地震发生后伴随有大量余震(韩佳东, 2018; 韦伟等, 2018; 尹昕忠等, 2018), 但目前没有针对米林地震余震序列震源机制解研究的报道, 故对余震静态库伦应力的作用仅限于定性讨论。 一般在研究主震对余震静态库伦应力的作用时, 均假定余震与主震的震源机制一致, 统计处于静态库伦破裂应力变化正值区(下文称之为促发区, 主震对促发区内的余震的作用称为促发作用)的余震比例(李瑶等, 2017; 靳志同等, 2019)。 这是在未知余震震源机制解情况下的一种妥协方式, 如果认为余震是主震破裂的延续, 则这种做法是合理的。 往往也会有一些余震处于静态库伦破裂应力变化的负值区(下文称之为抑制区, 主震对抑制区内的余震的作用称为抑制作用), 这部分余震是否确实被抑制?既然被抑制, 为什么还会发生?若这部分余震的震源机制与主震不同, 是否将得到相反的结果?对于这些问题, 目前还未见有相关的详细论述。 本文将基于米林地震的余震数据深入研究此问题, 这也是此文立意之初衷。

1 数据

为了精确地研究米林地震对后续地震静态库伦应力的影响, 需要利用米林地震的破裂模型、 后续地震的精定位数据及其震源机制解。 然而目前无法得到米林地震余震序列的震源机制解, 主要因为: 1)米林地震序列中绝大部分余震的震级较小(图3b), 无论使用哪种方法求解余震的震源机制解, 结果中均包含很大的误差(这是主要原因); 2)余震的数量太多, 求解震源机制解工作量太大。 而墨脱地震震级较大, 可以得到可信的震源机制解。

1.1 米林地震的破裂模型

《2017年11月18日西藏米林6.9级地震图集》中给出了根据远场体波资料反演得到的米林地震震源断层走向(109° )、 倾角(27° )不变, 滑动角可变的破裂模型。 之后, Li等(2019)基于远场体波和InSAR资料也给出了米林地震震源断层走向(109° )、 倾角(29° )不变, 滑动角可变的破裂分布。 由于Li等(2019)的结果另外使用了InSAR资料作为约束, 具有更高的准确度, 因此本研究将采用Li等(2019)给出的破裂模型, 如图 1 所示。 破裂面沿走向长51km, 沿倾向长27km, 图中每个小块对应1个3km× 3km的子断层, 第一行子断层的中心深度为1.5km, 破裂模型给出了各子断层的滑动方向和滑动量。 小箭头所指方向表示子断层的滑动方向, 箭杆长短和子断层颜色表示滑动量的大小。 由图 1 可见, 破裂主要集中在中部, 中部左、 右两侧的破裂方式差异较大, 左侧以逆冲为主, 右侧以走滑为主。

图 1 米林地震的破裂模型(数据来源于Li et al., 2019)Fig. 1 Rupture model of the Mainling earthquake(The data comes from Li et al., 2019).

下文的研究需要用到米林地震的滑动角, 那么如何根据破裂模型给出一个恰当的平均滑动角呢?本文的基本思路是: 对于滑动量大的子断层, 其滑动角应占有较大权重。 具体做法是: 首先, 将各子断层的滑动方向变换为北东下(NED)坐标系下表达的滑动矢量, 将子断层的滑动量作为滑动矢量的模长; 然后, 将各子断层的滑动矢量叠加得到和矢量, 以和矢量的方向作为米林地震破裂模型的平均滑动方向, 平均滑动方向和走向的夹角即为平均滑动角。 用上述方法计算得到米林地震的平均滑动角为42.33° , 说明整体上米林地震破裂的走滑分量与逆冲分量相当。 震源破裂过程(Li et al., 2019)显示破裂自左侧逆冲部分开始向右侧发展, 故《2017年11月18日西藏米林6.9级地震图集》和其他学者(白玲等, 2017; 吴宝峰, 2017; 刘云华等, 2018)给出的米林地震震源机制以逆冲为主。

1.2 米林地震的余震序列

2016年中国地震局地球物理研究所在南迦巴瓦地区布设了流动地震台网, 运行期间发生了米林地震。 韩佳东(2018)基于流动地震台网观测资料, 利用双差定位方法得到了2017年11月18日— 2017年12月17日米林地震及其余震共1 508个事件的精定位结果, 震中分布见图 2。 由图 2 可见, 余震震中分布的走向与破裂模型给出破裂面的走向一致。 图 3 统计了米林地震余震的震源深度和震级。 由图 3 可见, 米林地震绝大部分余震发生在0~24km深度范围内, 大部分余震震级为ML1~4。

图 2 米林地震及其后续地震的震中分布
红色五角星表示精定位后的米林地震震中, 黑色小圆圈表示米林地震余震的震中, 与五角星相连的图形为破裂模型给出的米林地震震源断层的下半球施密特投影; 墨脱地震的震源机制解来自于中国地震台网中心
Fig. 2 The epicenter distribution of the Mainling earthquake and its subsequent earthquakes.

图 3 米林地震余震的深度(a)和震级(b)统计Fig. 3 Depth and magnitude statistics of aftershocks of the Mainling earthquake.

1.3 墨脱地震的震源机制解

中国地震台网中心发布的《2019年4月24日西藏墨脱6.3级地震图集》和一些其它机构发布的墨脱M6.3地震的震源机制解见表1和图 4。 从图 4 可见, 不同震源机制解的2个节面的方位分布比较集中。 为了将这些震源机制解之间的差别定量化, 把这些震源机制解两两组合, 对所有组合方式的最小空间旋转角(Kagan, 1991; 万永革, 2019)求均值, 得到平均最小空间旋转角为11.29° 。 图 4 和平均最小空间旋转角表明这些由不同资料和方法得到的震源机制解之间的差别较小, 互相印证了这些震源机制解的可靠性。 尽管不同机构或研究者根据不同资料和方法给出的墨脱M6.3地震的震源机制解差别较小, 但下文定量计算米林地震对墨脱地震的静态库伦应力影响需要使用墨脱地震确切的断层面参数, 由此将面临着选择数据的困难。

图 4 墨脱地震不同震源机制解及其中心震源机制解的下半球施密特投影
图中紫色弧线表示墨脱地震不同的震源机制节面; 黑色弧线表示中心震源机制解的2个节面, 绿色弧线覆盖区域为其不确定范围; 红色、 蓝色和黄色的点表示中心震源机制解的P轴、 T轴和B轴, 其周围对应颜色的封闭曲线表示其不确定性范围; 绿点、 黑点和青色的点表示各个机构给出震源机制的P轴、 T轴和B轴的投影
Fig. 4 Lower hemisphere Schmidt projection of different focal mechanisms of the Medog earthquake and its central focal mechanism.

本节利用万永革(2019)根据同一地震多个震源机制解确定一个与所有震源机制解差别最小的 “ 中心震源机制解” 的方法求取墨脱地震的中心震源机制解。 该方法将不同震源机制解之间的最小空间旋转角作为它们的差别, 以某个震源机制解作为初始解, 在初始解的邻域内进行一阶泰勒展开, 采用Levenberg-Marquardt方法(Levenberg, 1944; Marquardt, 1963)迭代求解一个与所有震源机制解最小空间旋转角的平方和最小的中心震源机制解。 分别以各震源机制解作为初始解得到对应的中心震源机制解, 发现不同中心震源机制解与所有震源机制解最小空间旋转角的标准差及其差别都很小(表1第四列), 说明不同初始解得到的结果基本一致。 以IPGP给出的震源机制解作为初始解得到的中心震源机制解的标准差最小, 故将该中心震源机制解作为最终结果, 即: 节面I和节面Ⅱ 的走向、 倾角和滑动角分别为91.50° 、 85.40° 、 96.62° 和216.15° 、 8.05° 、 34.91° ; P轴的走向和倾伏角分别为175.39° 和40.05° , 不确定范围分别为167.50° ~183.50° 和31.99° ~47.89° ; T轴的走向和倾伏角分别为8.66° 和49.18° , 不确定范围分别为0.12° ~16.96° 和41.37° ~57.24° ; B轴的走向和倾伏角分别为270.97° 和6.60° , 不确定范围分别为263.08° ~279.08° 和-1.49° ~14.46° 。 墨脱地震的中心震源机制解及其不确定范围见图 4, 中心震源机制解与不同机构或研究者给出的震源机制解的最小空间旋转角最小为1.71° , 最大为15.93° (表1第五列)。

2 方法

Okada(1992)给出了均匀弹性半空间模型下有限尺度断层产生的应变场解析表达式, 通过应力、 应变之间的本构关系(徐芝纶, 2006)可得到弹性半空间中的应力场。 将1.1节所述的米林地震破裂模型中每个子断层在半空间产生的应力场叠加, 就得到米林地震产生的应力场。 判定米林地震对后续地震影响的方法是计算米林地震在后续地震断层面上产生的静态库伦破裂应力变化, 相关研究利用式(1)计算静态库伦破裂应力变化, 它描述了主震释放的同震应力场对断层错动的贡献:

ΔCFS=Δτ+μ'Δσn(1)

式中, Δ CFS表示接收断层上产生的静态库伦破裂应力变化; Δ τ 为滑动方向上的剪应力变化; Δ σ n为正应力变化, 拉张为正; μ '为视摩擦系数, 一般取0.2~0.8(Stein, 1999)。 静态库伦破裂应力变化为正, 表示对断层错动有促进作用; 静态库伦破裂应力变化为负, 表示对断层错动有抑制作用。

接收断层上的静态库伦破裂应力变化是正应力变化作用和滑动方向上剪应力变化作用的综合效果。 式(1)表明无论Δ τ 为正或负, Δ σ n为正(拉张)均有助于断层错动, Δ σ n为负(挤压)则抑制断层错动。 实际上, 在Δ τ 为正或负2种情况下, Δ σ n的正负对断层错动的作用不同, 下面将结合图 5 进行具体分析。 由图 5 可见, Δ τ 是最大剪应力变化(Δ σ s)的分量(Δ τ σ scosα )。 当α ≤ 90° 时, Δ τ ≥ 0, Δ τ 的方向与滑动方向一致。 若Δ σ n> 0, 正应力表现为拉张作用, μ σ n有助于断层向滑动方向错动; 若Δ σ n< 0, 正应力表现为挤压作用, μ σ n抑制断层向滑动方向错动。 当α > 90° 时, Δ τ < 0, Δ τ 的方向与滑动方向相反。 若Δ σ n> 0, 正应力表现为拉张作用, μ σ n有利于断层向滑动方向的反方向错动, 也就抑制了断层向滑动方向错动; 若Δ σ n< 0, 正应力表现为挤压作用, μ σ n抑制断层向滑动方向的反方向错动, 但正应力的挤压作用只能相对减弱断层向滑动方向反方向错动的趋势, 对滑动方向不会有促进作用。 总之, 正应力为正, 有利于断层向滑动方向所在直线上剪应力变化为正的方向滑动。 万永革(2001)曾在研究中考虑到上述问题, 将静态库伦破裂应力变化计算公式定义为式(2), 本研究将利用式(2)研究米林地震对其余震和墨脱地震的静态应力影响, 仿照前人的研究(Stein et al., 1992; King et al., 1994; 万永革等, 2010), 取μ '=0.4。

ΔCFS=Δτ+μ'Δσn(Δτ0)Δτ-μ'Δσn(Δτ< 0, Δτ-μ'Δσn0)0(Δτ< 0, Δτ-μ'Δσn> 0)(2)

图 5 断层面上正应力变化(Δ σ n)、 最大剪应力变化(Δ σ s)和滑动方向上剪应力变化(Δ τ )的关系示意图Fig. 5 Diagram of relationship between normal stress change(Δ σ n), maximum shear stress change(Δ σ s) and shear stress change in the slide direction(Δ τ )of the fault.

要计算式(2)表达的静态库伦破裂应力变化, 必须明确掌握接收断层面的参数。 由于米林地震余震的震源机制解未知, 因此无法确切计算米林地震对余震的静态应力影响。 本研究首先仿照其他研究者(李瑶等, 2017; 靳志同等, 2019)的做法, 设置主震周围空间接收断层的参数与主震一致, 观察主震周围的Δ CFS分布以及处于促发区和抑制区余震的比例。 余震在主震断层面附近丛集, 若认为余震是主震破裂的延续, 这种做法是合理的。 然而, 并不是所有余震都是主震的继续破裂所导致的, 也有一些余震发生在主震断层面附近的破碎带内。 这些余震的断层方位是随机的, 和主震断层面的方位有较大的差异, 这也许是一部分余震处于抑制区的主要原因。 那么抑制区的余震有没有被促发的可能性?本文的另一种做法是计算在主震同震应力场的作用下空间中每个点最容易发生错动的断层方位, 该断层上的Δ CFS最大, 称为该点的最大静态库伦破裂应力变化(Δ CFSmax)。 计算每个点不同方位断层上的Δ CFS时, 把最大剪应力变化的方向作为滑动方向。 如图 5 所示, 滑动方向上的剪应力变化是最大剪应力变化的分量, 结合式(2)可知, 将最大剪应力变化的方向作为滑动方向可以保证各点的Δ CFSmax最大。 本研究感兴趣的是Δ CFSmax大于触发阈值0.01MPa的区域(本文称为触发区, 主震对触发区内地震的作用称为触发作用)。 主震对触发区外的地震可能存在抑制作用, 也可能存在促发作用。 但即便存在促发作用, 促发作用也达不到触发的量级, 或者说主震不是余震发生的直接原因。 触发区内的地震有被主震触发的可能性, 取决于地震断层的空间取向, Δ CFSmax越大, 满足被触发的断层方位越多, 被触发的可能性也就越大。 实际上, Δ CFSmax描述了主震能够触发余震的最大空间范围。 研究米林地震对墨脱地震的影响, 首先可以根据墨脱地震的震中位置和主震能够触发的范围开展定性分析, 若墨脱地震处于能被触发的范围之内, 说明墨脱地震有被触发的可能性; 若墨脱地震处于能被触发的范围之外, 则可以肯定米林地震不是墨脱地震发生的直接原因, 具体影响仍需结合墨脱地震的震源机制解和式(2)进行定量计算。

3 结果

根据Okada(1992)的公式和应力、 应变之间的本构关系计算有限断层面产生的应力场时, 取介质为泊松介质(泊松比为0.25), 介质的剪切模量取为地壳的剪切模量, 即3× 1010Pa。 本文首先计算接收断层参数与破裂模型(走向、 倾角和滑动角分别为109° 、 29° 和42.33° )一致时静态库伦破裂应力变化的空间分布; 其次, 把主震同震应力场作用下最容易错动的断层面作为接收断层, 得到最大静态库伦破裂应力变化的空间分布。 之后, 结合这2个结果分析米林地震对余震静态库伦应力的影响; 最后, 根据最大静态库伦破裂应力变化的空间分布情况和墨脱地震2个节面上产生的静态库伦破裂应力变化, 分析米林地震对墨脱地震的影响。

3.1 米林地震对其余震的静态库伦应力的影响

假定空间中接收断层的参数与主震一致, 根据图3a米林地震余震的震源深度分布情况, 分别计算了3km、 9km、 15km和21km 4个深度(Z)上的静态库伦破裂应力变化, 结果见图 6, 各子图中的灰色小圆圈代表对应深度范围[Z-3, Z+3]的余震。 从图 6 可见, 4个深度上的静态库伦破裂应力变化的空间分布不同, Z=21km时静态库伦破裂应力变化的高值区面积最小, 这是因为21km深度处于破裂面的下方, 距离破裂面较远。 在Z=3km深度上余震处于促发区的比例为36.16%; 在Z=9km深度上余震处于促发区的比例为21.57%; 在Z=15km深度上余震处于促发区的比例为39.92%; 在Z=21km深度上余震处于促发区的比例为18.80%。 这些数据说明不同深度上仅有少部分余震发生在促发区, 考虑到所用计算公式不同导致的结果差异, 作者还利用式(1)计算了4个深度上产生的静态库伦破裂应力变化的空间分布, 发现二者所得结果基本一致, 排除了计算方法对结果的影响。 若不考虑计算静态库伦破裂应力变化所用模型和数据的误差, 造成该现象的主要原因有: 1)静态库伦破裂应力变化是按一定深度计算的, 而对应余震是按深度范围给出的, 两者之间在对应上存在一定偏差; 2)抑制区余震的发震方式(指震源机制)和主震不一致; 3)未考虑时间上先发生的余震对后发生余震产生的应力影响。 对于第一种因素, 由于本文采用的是双差精定位结果, 这种偏差应该很小; 对于第3种因素, 万永革等(2005)已验证小余震对后续地球物理场的影响与主震相比可以忽略。 故探索余震为何处于抑制区, 仅需着重考虑第二个因素。

图 6 米林地震在不同深度上产生的静态库伦破裂应力变化
白色区域的静态库伦破裂应力变化为负, 黑色矩形是破裂面在水平面的投影, 灰色小圆圈表示对应深度的余震, 各子图中的余震深度为[Z-3, Z+3]
Fig. 6 Static Coulomb failure stress change at different depths caused by the Mainling earthquake.

如果促发区内余震的实际发震方式和主震相似, 可认为是主震继续破裂所导致的; 另外的一大部分余震发生在抑制区, 这些余震和促发区内的余震在时空分布上并没有明显的区分, 因此这些余震可能并非是被抑制的, 否则它们也不会在主震后短时间内大量发生。 如果发生在抑制区的余震和主震的发震方式不同是否将得到不同的结果?实际上, 由于地震断层面附近的破裂较为剧烈, 会在主震断层面附近产生破碎带, 破碎带内存在各种方位的薄弱面, 在米林地震同震应力场的作用下, 薄弱面有向同震应力场投影于薄弱面的剪应力方向错动的趋势, 在某些方位的薄弱面上静态库伦破裂应力变化达到使断层错动的阈值就促发了余震; 各种方位薄弱面的存在使得余震更容易发生, 震源区应力恢复过程中产生的应力扰动可能导致余震发生; 如果主震未将构造背景积累的应力充分释放, 剩余的背景应力作用于薄弱面也会导致余震发生。 上述3种情况中, 余震和主震发震方式的关系是相互独立的, 第一种情况主震对余震起促发作用, 后2种情况主震对余震的静态应力影响需结合余震的发震机制进行分析。 另一方面, 各种方位薄弱面的存在使得近距离余震的发震方式呈现出多样性(万永革, 2020)。

为了研究抑制区内的余震是否有被促发的可能性, 本文计算了3km、 9km、 15km和21km 4个深度(Z)上如第二节所述的最大静态库伦应力变化(Δ CFSmax), 各空间点断层的走向和倾角以5° 为间隔进行试算。 图 7 给出了4个深度上触发区的范围, 若地震位于该范围外, 即使米林地震对该地震具有促发作用, 也会由于作用太小而可以被忽略; 若地震位于该范围内, 则有米林地震触发该地震的可能性, 并且该触发作用是地震发生的直接原因。 不同深度上的触发区围绕破裂面分布, 距离破裂面越近, Δ CFSmax越大。 在深部(Z=21km)由于远离破裂面, 触发区的范围变小, Δ CFSmax的2个高值区对应于破裂面中部破裂最剧烈的逆冲部分和走滑部分(图7d), 这些结果与预期一致。 不同深度所有的余震都处于触发区内, 并且绝大部分余震处于最大静态库伦应力变化的高值区。 最大静态库伦应力变化越大, 满足被触发的断层方位越多, 被触发的可能性越大。 图 7 说明所有余震均有被触发的可能性, 佐证了根据图 6 结果猜测的抑制区大量余震也可能受到主震促发作用影响, 而处于抑制区是由于余震的发震方式和主震不同。 以上所做阐述都是可能性推测, 由于余震震级太小, 难以得知其震源机制解, 无法准确计算米林地震对其余震的真正作用。

图 7 米林地震在不同深度上产生的最大静态库伦破裂应力变化
白色区域的最大静态库伦破裂应力变化< 0.01MPa, 黑色矩形是破裂面在水平面的投影, 灰色小圆圈表示对应深度的余震, 各子图中的余震深度为[Z-3, Z+3]
Fig. 7 Maximum static Coulomb stress failure change at different depths caused by the Mainling earthquake.

米林地震释放的同震应力场在空间上的变化是渐变的, 空间上的最大静态库伦破裂应力变化所对应的接收断层面参数的变化也应该呈现渐变的形态。 为了检验图 7 结果中接收断层面选择的合理性, 以图7b对应的接收断层面参数(走向、 倾角、 滑动角)为例, 将其绘于图 8。 图 8 显示, 接收断层面参数对应的震源机制在空间的变化呈现出很强的规律性: NE-SW向为走滑型, NW-SE向为逆冲型, 其它位置的震源机制为这2种类型的过渡。 然而图 8 展示的接收断层的方位随空间的渐变效果并不尽如人意。 图 8 中洋红色沙滩球对应的接收断层方位与其周围接收断层的方位差别很大, 以此为例研究小区域内接收断层方位突变的原因。 提取出米林地震在图 8 中洋红色沙滩球的位置(29.3° N, 94.6° E)产生的应力张量, 走向和倾角均以5° 为间隔, 遍历计算不同方位断层的最大静态库伦破裂变化, 结果见图 9。 图 9 显示, 最大静态库伦破裂应力变化值与断层方位并非一一对应, 一个最大静态库伦破裂应力变化值可对应多个断层方位。 图 8 中洋红色沙滩球附近的接收断层面包括倾角接近90° 、 走向大致为280° 和350° 的2类, 图 9 中这2类断层上的最大静态库伦破裂变化基本一致, 很好地解释了图 8 中接收断层在空间上方位突变的原因。 图 8 与图 9 的结果说明, 计算最大静态库伦破裂应力变化选择的接收断层面是合理的。

图 8 与图7b对应的接收断层面的参数(走向、 倾角和滑动角)随空间的变化
图中沙滩球的投影方式为下半球施密特投影, 填充部分为压缩象限; 黑色弧线为最大静态库伦破裂应力变化对应的接收断层面
Fig. 8 Variation of parameters of receiving fault(strike, dip and rake)corresponding to Figure 7b.

图 9 根据米林地震在图 8 中洋红色沙滩球的位置产生的应力张量计算得到的不同方位断层上的最大静态库伦破裂应力变化
图中的洋红色沙滩球与图 8 相同, 其坐标为接收断层面的方位
Fig. 9 Maximal static Coulomb failure stress change on faults with different orientation calculated based on the stress tensor generated by the Mainling earthquake at the location of the magenta beach ball in Figure 8.

本节采用2种接收断层面选择方案研究米林地震对余震静态库伦应力影响的结果说明, 余震很可能发生于主震破裂面附近剧烈破碎产生的破碎带内各种方位的薄弱面, 这种情况下, 不能将主震的断层参数作为优势断层(接收断层)的参数来研究主震对断层区(震源区)内余震的静态库伦应力影响。 主震断层面附近的剧烈破碎会导致断层区内的弹性常数大幅度降低。 万永革等(2008)曾在唐山地震序列应力触发的研究中考虑了断层区内与断层区外的弹性常数的区别; Li等(1990, 1994)利用围陷波给出了断层区内弹性常数的取值。 这些研究均支持本节结果, 预示着研究断层区内地震之间的静态库伦应力影响需要考虑断层区内与断层区外弹性常数的差异, 获得断层区内准确的弹性常数对正确估计断层区内地震之间的静态库伦应力影响是有益的。

3.2 米林地震对墨脱地震的静态库伦应力影响

中国地震台网中心给出墨脱地震震中为(28.4° N, 94.61° E), 深度为10km。 从图7b米林地震在9km深度上产生的触发区范围来看, 墨脱地震的位置距离触发区较远, 说明米林地震对墨脱地震不存在触发作用, 但不能确定是否存在较弱的促发作用或抑制作用。 在已知墨脱地震断层参数的情况下, 最直接的方法是计算米林地震在墨脱地震断层面上产生的静态库伦破裂应力变化。 由于震源机制解给出了2个等可能的发震断层面, 本节将分别对2个断层面计算, 根据计算结果再开展进一步研究。

根据米林地震破裂模型, 计算得到米林地震在墨脱地震震源处产生的同震应力张量为

σ=881.98143.14-95.09143.14-90.05-11.05-95.09-11.058.01(3)

在本文1.3节中, 已根据不同机构和作者利用不同资料和方法得到的墨脱地震震源机制解求取了其中心震源机制解, 下面将根据中心震源机制解给出的断层面参数计算米林地震对墨脱地震2个节面的静态库伦应力影响, 计算得到2个节面的正应力变化Δ σ n、 滑动方向上的剪应力变化Δ τ σ 作用在节面上的最大剪应力方向与震源机制给出滑动方向之间的夹角α 以及静态库伦破裂应力变化(Δ CFS)等参数见表2。 从表2所示的结果来看, 米林地震对2个节面都表现出抑制作用, 且量值都很小, 因此可忽略其对墨脱地震的影响, 即认为墨脱地震是一次独立于米林地震的事件。

表2 墨脱地震2个节面的Δ σ n、 Δ τ α 和Δ CFS Table2 Δ σ n, Δ τ , α and Δ CFS obtained from two nodal planes of the Medog earthquake
4 结论与讨论

本文利用米林地震的破裂模型、 米林地震余震的精定位数据以及墨脱地震的震源机制解, 计算了均匀弹性半空间模型下基于2种接收断层选择方案的米林地震在其附近空间产生的静态库伦破裂应力变化, 讨论了米林地震对其余震以及墨脱地震的影响, 主要得到以下几点认识:

(1)设置接收断层参数和主震一致得到不同深度的静态库伦破裂应力变化结果显示, 不同深度仅有少部分余震处于促发区, 大部分余震处于抑制区。 若促发区余震的实际震源机制和主震相似, 可认为是主震的继续破裂引发余震, 分析认为大量余震处于抑制区的主要原因可能是它们发生在主震剧烈破裂导致的破碎带内的各种方位的薄弱面上, 和主震的震源机制有较大差异。

(2)假设破碎带内断层的方位是随机的, 计算了在主震影响下, 主震附近不同深度最可能破裂的断层上产生的静态库伦破裂应力变化, 即最大静态库伦破裂应力变化。 结果发现, 不同深度上所有余震均位于最大静态库伦破裂应力变化大于触发阈值0.01MPa的范围内, 说明不确定震源机制的情况下所有余震均有被主震促发的可能性, 因此有理由认为在接收断层和主震一致的假设下抑制区的大量余震并非都是被抑制的, 否则无法解释它们在主震后短期内大量发生。

(3)墨脱地震震源位于米林地震最大静态库伦破裂应力变化< 0.01MPa的区域内, 说明米林地震对墨脱地震不存在触发作用。 根据墨脱地震的中心震源机制解, 定量计算了2个节面上产生的静态库伦破裂应力变化, 结果表明2个节面上的静态库伦破裂应力变化为负且量值很小, 说明米林地震没有影响到墨脱地震的孕育进程。 利用不同资料和方法得到的墨脱地震震源机制之间存在的差异无疑会给米林地震对墨脱地震静态库伦应力影响的计算结果带来不确定性。 1.3节的结果说明利用不同资料和方法所得墨脱地震震源机制之间的差别很小, 因而将它们的中心震源机制解作为代表来研究米林地震对墨脱地震的影响是可行的, 以其它震源机制解给出的断层参数来研究米林地震对墨脱地震的影响与根据中心震源机制解得到的结论一致。

本文以米林地震对后续地震的静态库伦应力影响为例, 重点分析了接收断层的不确定性为结果带来的影响。 在接收断层参数与主震一致的假设下, 本文的震例与一些震例(姜辉等, 2011; 靳志同等, 2019)的结果相似, 均发现很大比例的余震处于抑制区, 推测其原因是由于余震发生于主震附近剧烈破裂导致的破碎带内各种方位的薄弱面, 发震方式与主震不一致。 米林地震的余震全部处于最大静态库伦破裂应力变化的触发区, 说明余震均有被促发的可能性。 这进一步说明由于主震破裂面附近的破碎剧烈导致断层区的弹性常数降低, 没有将断层区内的弹性常数区别对待以及将主震的断层参数作为优势断层来研究主震对断层区内余震的静态库伦应力影响是不恰当的。 这也预示着在断层区内进行地震之间的静态库伦应力影响的研究需要准确地评估断层区内的弹性常数。 当然, 也有一些震例在接收断层参数与主震一致的假设下得到的静态库伦破裂应力变化与余震的空间分布对应较为理想, 如2016年12月8日新疆呼图壁MS6.2地震(王苏等, 2018)等, 这可能由于震源区的介质坚硬, 没有很好地发育破碎带, 绝大多数余震发生于主震破裂面, 也有可能是余震的发震机制与主震并不一致, 只是这些丛集的余震距离主震非常近, 余震既处于接收断层参数与主震一致假设下的促发区, 也处于最大静态库伦破裂的触发区。 例如, 王苏等(2018)在其研究中利用与计算2016年12月8日新疆呼图壁MS6.2地震产生的静态库伦应力变化完全相同的方法对2016年11月25日新疆阿克陶地震开展分析, 得到的结果却与本文结果类似, 即余震处于促发区的比例并不理想。 当然, 还有可能一部分余震由主震继续破裂导致, 发震方式与主震一致, 另一部分余震发生于破碎带内, 发震方式与主震不一定相同。

若主震将构造背景积蓄的能量充分释放, 在震源附近余震迅速发生的短时间内背景应力和应力恢复过程中产生的应力扰动的作用较小, 余震主要受主震同震应力的作用, 主震对大部分余震起到促发作用, 本文的数据和结果倾向于这种情况。 若主震未能将构造背景积蓄的能量充分释放, 能量的继续释放会导致较大震级的余震, 如果余震仍沿着主震断层面破裂, 余震的震源机制则与主震有较好的一致性。 这种情况下余震将受控于背景构造应力, 主震产生的同震应力对余震的发生可能起到促进或者延缓作用, 具体影响需要根据余震的震源机制来分析。 在余震的震源机制与主震相似的情况下, 假设接收断层参数与主震一致所得到的静态库伦破裂应力变化能较为真实地反映主震对余震的作用, 如2015年尼泊尔MW7.9地震(黄禄渊等, 2016)等。 一般由于余震震级小、 数量多, 无法准确给出震源机制解, 因此只能在假设余震震源机制的情况下定性讨论主震释放的同震应力场对余震的影响, 这也是本研究的主要工作。 不能把主震对小震级余震的作用定量化, 是进行主震对余震静态库伦应力影响研究的极大挑战, 也是最大的缺陷, 深入研究余震的发震机理或发展准确求解小地震震源机制的方法是弥补该缺陷的2个可行的方向。

致谢 北京大学张勇教授为本研究提供了米林地震的破裂模型; 刘敬光提供了《2019年4月24日西藏墨脱6.3级地震图集》; 审稿专家提出的建设性修改意见使本研究得以进一步完善。 在此一并表示感谢!

参考文献
[1] 白玲, 李国辉, 宋博文. 2017. 2017年西藏米林6. 9级地震震源参数及其构造意义[J]. 地球物理学报, 60(12): 49564963.
BAI Ling, LI Guo-hui, SONG Bo-wen. 2017. The source parameters of the M6. 9 Mainling, Tibet earthquake and its tectonic implications[J]. Chinese Journal of Geophysics, 60(12): 49564963(in Chinese). [本文引用:3]
[2] 韩佳东. 2018. 米林地震序列监测与南迦巴瓦地震活动性研究[D]. 北京: 中国地震局地球物理研究所.
HAN Jia-dong. 2018. Mainling M6. 9 earthquake sequences monitored and Namche Barwa seismicity analysis [D]. Institute of Geophysics, China Earthquake Administration, Beijing(in Chinese). [本文引用:2]
[3] 黄禄渊, 杨树新, 陈连旺, . 2016. 2015年尼泊尔 MW7. 9地震对强余震的触发作用和对周边断层的影响[J]. 地震研究, 39(2): 187195.
HUANG Lu-yuan, YANG Shu-xin, CHEN Lian-wang, et al. 2016. Strong aftershocks triggering and stress effect on nearby faults induced by Nepal MW7. 9 earthquake in 2015[J]. Journal of Seismological Research, 39(2): 187195(in Chinese). [本文引用:1]
[4] 姜辉, 邓志辉, 王海涛, . 2011. 从新疆地区主余震活动看地震静态应力触发模型[J]. 地震地质, 33(3): 586601. doi: 103969/j.issn.0253-4967.2011.03.009.
JIANG Hui, DENG Zhi-hui, WANG Hai-tao, et al. 2011. Study of static stress triggering model through the activities of main and aftershocks in Xinjiang[J]. Seismology and Geology, 33(3): 586601(in Chinese). [本文引用:1]
[5] 靳志同, 万永革, 刘兆才, . 2019. 2017年九寨沟 MS7. 0地震对周围地区的静态应力影响[J]. 地球物理学报, 62(4): 12821299.
JIN Zhi-tong, WAN Yong-ge, LIU Zhao-cai, et al. 2019. The static stress triggering influences of the 2017 MS7. 0 Jiuzhaigou earthquake on neighboring areas[J]. Chinese Journal of Geophysics, 62(4): 12821299(in Chinese). [本文引用:3]
[6] 李瑶, 万永革, 靳志同, . 2017. 新疆精河 MW6. 3地震产生的静态应力变化研究[J]. 中国地震, 33(4): 671681.
LI Yao, WAN Yong-ge, JIN Zhi-tong, et al. 2017. The study of static stress variation of the MW6. 3 Jinghe earthquake[J]. Earthquake Research in China, 33(4): 671681(in Chinese). [本文引用:2]
[7] 李煜航. 2018. 青藏高原东北横向扩展运动研究[J]. 国际地震动态, (5): 4547.
LI Yu-hang. 2018. Study on the lateral motion of northeastern Tibetan plateau[J]. Recent Developments in World Seismology, (5): 4547(in Chinese). [本文引用:1]
[8] 刘云华, 单新建, 张迎峰, . 2018. 基于地震波及InSAR数据的2017年11月18日西藏米林 MS6. 9地震发震构造[J]. 地震地质, 40(6): 12541275.
LIU Yun-hua, SHAN Xin-jian, ZHANG Ying-feng, et al. 2018. Use of seismic waveforms and InSAR data for determination of the seismotectonics of the Mainling MS6. 9 earthquake on Nov. 18, 2017[J]. Seismology and Geology, 40(6): 12541275(in Chinese). [本文引用:2]
[9] 万永革. 2001. “地震静态应力触发”问题的研究[D]. 北京: 中国地震局地球物理研究所.
WAN Yong-ge. 2001. Study on “static stress triggering” problem[D]. Institute of Geophysics, China Earthquake Administration, Beijing(in Chinese). [本文引用:1]
[10] 万永革. 2019. 同一地震多个震源机制中心解的确定[J]. 地球物理学报, 62(12): 47184728.
WAN Yong-ge. 2019. Determination of center of several focal mechanisms of the same earthquake[J]. Chinese Journal of Geophysics, 62(12): 47184728(in Chinese). [本文引用:2]
[11] 万永革. 2020. 震源机制与应力体系关系模拟研究[J]. 地球物理学报, 63(6): 22812296.
WAN Yong-ge. 2020. Simulation on relationship between stress regimes and focal mechanisms of earthquakes[J]. Chinese Journal of Geophysics, 63(6): 22812296(in Chinese). [本文引用:1]
[12] 万永革, 沈正康, 兰从欣. 2005. 兰德斯地震断层面及其附近余震产生的位移场研究[J]. 地震学报, 27(2): 139146.
WAN Yong-ge, SHEN Zheng-kang, LAN Cong-xin. 2005. Study on displacement field generated by aftershocks in Land ers seismic fault plane and its adjacent areas[J]. Acta Seismologica Sinica, 27(2): 139146(in Chinese). [本文引用:1]
[13] 万永革, 沈正康, 盛书中, . 2010. 2008年新疆于田7. 3级地震对周围断层的影响及其正断层机制的区域构造解释[J]. 地球物理学报, 53(2): 280289.
WAN Yong-ge, SHEN Zheng-kang, SHENG Shu-zhong, et al. 2010. The mechanical effects of the 2008 MS7. 3 Yutian, Xinjiang earthquake on the neighboring faults and its tectonic origin of normal faulting mechanism[J]. Chinese Journal of Geophysics, 53(2): 280289(in Chinese). [本文引用:1]
[14] 万永革, 沈正康, 曾跃华, . 2008. 唐山地震序列应力触发的黏弹性力学模型研究[J]. 地震学报, 30(6): 581593.
WAN Yong-ge, SHEN Zheng-kang, ZENG Yue-hua, et al. 2008. Study on visco-elastic stress triggering model of the 1976 Tangshan earthquake sequence[J]. Acta Seismologica Sinica, 30(6): 581593(in Chinese). [本文引用:1]
[15] 王苏, 李建有, 徐晓雅, . 2018. 2016年新疆阿克陶 MS6. 7地震和呼图壁 MS6. 2地震的余震触发研究[J]. 地震研究, 41(1): 98103.
WANG Su, LI Jian-you, XU Xiao-ya, et al. 2018. Study on aftershock triggering of the 2016 Aktao MS6. 7 and Hutubi MS6. 2 earthquakes in Xinjiang[J]. Journal of Seismological Research, 41(1): 98103(in Chinese). [本文引用:2]
[16] 韦伟, 谢超, 周本刚, . 2018. 西藏米林 M6. 9地震及其余震序列地震定位[J]. 科学通报, 63(15): 14931501.
WEI Wei, XIE Chao, ZHOU Ben-gang, et al. 2018. Location of the mainshock and aftershock sequences of the M6. 9 Mainling earthquake, Tibet[J]. Chinese Science Bulletin, 63(15): 14931501(in Chinese). [本文引用:1]
[17] 吴宝峰. 2017. 2017年11月18日西藏米林6. 9级地震震源机制[J]. 地震地磁观测与研究, 38(6): 2629.
WU Bao-feng. 2017. Determination of the focal mechanism of Milin M6. 9 earthquake on Nov. 18, 2017 in Tibet[J]. Seismological and Geomagnetic Observation and Research, 38(6): 2629(in Chinese). [本文引用:2]
[18] 徐芝纶. 2006. 弹性力学[M]. 北京: 高等教育出版社.
XU Zhi-lun. 2006. Elastomechanics[M]. Higher Education Press, Beijing(in Chinese). [本文引用:1]
[19] 尹昕忠, 周本刚, 陈九辉, . 2018. 西藏米林 M6. 9地震早期余震时空分布特征[J]. 地球物理学报, 61(6): 23222331.
YIN Xin-zhong, ZHOU Ben-gang, CHEN Jiu-hui, et al. 2018. Spatial-temporal distribution characteristics of early aftershocks following the M6. 9 Mainling earthquake in Tibet, China[J]. Chinese Journal of Geophysics, 61(6): 23222331(in Chinese). [本文引用:1]
[20] Kagan Y Y. 1991. 3-D rotation of double-couple earthquake sources[J]. Geophysical Journal International, 106(3): 709716. [本文引用:1]
[21] King G C P, Stein R S, Lin J. 1994. Static stress changes and the triggering of earthquakes[J]. Bulletin of the Seismological Society of America, 84(3): 935953. [本文引用:1]
[22] Levenberg K. 1944. A method for the solution of certain non-linear problems in least squares[J]. Quarterly of Applied Mathematics, 2(2): 164168. [本文引用:1]
[23] Li W, Xu C J, Yi L, et al. 2019. Source parameters and seismogenic structure of the 2017 MW6. 5 Mainling earthquake in the Eastern Himalayan Syntaxis(Tibet, China)[J]. Journal of Asian Earth Sciences, 169: 130138. [本文引用:4]
[24] Li Y G, Leary P, Aki K, et al. 1990. Seismic trapped modes in the Oroville and San Andreas fault zones[J]. Science, 249(4970): 763766. [本文引用:1]
[25] Li Y G, Vidale J E, Aki K, et al. 1994. Fine structure of the Land ers fault zone: Segmentation and the rupture process[J]. Science, 265(5170): 367370. [本文引用:1]
[26] Marquardt D W. 1963. An algorithm for least-squares estimation of nonlinear parameters[J]. Journal of the Society for Industrial and Applied Mathematics, 11(2): 431441. [本文引用:1]
[27] Okada Y. 1992. Internal deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 82(2): 10181040. [本文引用:2]
[28] Shen X Z, Liu M, Gao Y, et al. 2017. Lithospheric structure across the northeastern margin of the Tibetan plateau: Implications for the plateau's lateral growth[J]. Earth and Planetary Science Letters, 459(1): 8092. [本文引用:1]
[29] Stein R S. 1999. The role of stress transfer in earthquake occurrence[J]. Nature, 402(6762): 605609. [本文引用:1]
[30] Stein R S, King G C P, Lin J. 1992. Change in failure stress on the southern San Andreas fault system caused by the 1992 magnitude=7. 4 Land ers earthquake[J]. Science, 258(5086): 13281332. [本文引用:1]