许 鑫,万永革,2,李振月,崔华伟,冯 淦,李 枭

(1.防灾科技学院,河北 三河 065201;2.河北省地震动力学重点实验室,河北 三河 065201;3.中国科学技术大学,安徽 合肥 230026;4.山东省地震局,山东 济南 250014)

0 引言

据中国台网中心(http://www.ceic.ac.cn/)测定,2021年5月21日当地时间21时48分云南漾濞县发生MS6.4地震(25.67°N,99.87°E,以下简称漾濞地震),震中距大理市区仅30 km。此次地震事件损失严重,已导致3人死亡多人受伤,多间房屋倒塌。漾濞地震发生前27分,在其东南向不远处发生了5.6级地震(25.63°N,99.92°E,根据地震成核理论可将此地震定义为漾濞地震的前震)。漾濞地震发生后,短时间内该地区余震不下上千次之多,截至2021年5月31日,记录到的余震最大震级达5.2级(25.59°N,99.97°E)。

前人研究经验表明,通过计算地震引起的库仑破裂应力变化可以有效分析前震-主震-余震序列间的相互联系,为后续地震的发震时间及地点提供参考价值[1-6]。万永革等计算唐山地震等4个大地震第一次子破裂事件的静态库仑应力变化的结果均表明对后续破裂存在触发作用[2];又通过计算汶川MW7.9地震对周围大断裂带的库仑应力变化,揭示了部分断裂上的大震发震时间受其影响被触发或抑制[3]。沈正康等[4]通过黏弹性介质模型计算出东昆仑断裂带上地震间的触发作用。Robinson等[5]在已知断层滑动性质条件下,可以有效预测余震的发震位置盛书中等[6]通过计算尼泊尔强震产生的库仑应力变化,分析了其对中国大陆的影响,发现仅新疆和西藏的小部分地区影响相对明显。本次漾濞地震对周围主要断裂带影响如何？该序列中前震、主震和余震间是否存在触发关系？这是人们关心的问题。

针对以上问题,本研究首先根据Wang等[7]给出的漾濞地震的破裂模型参数,结合Okada[8]解析式计算此次漾濞地震对周围断层产生的静态库仑破裂应力变化以及对周围地区产生的水平面上的应力影响。然后计算了前震在主震破裂面和滑动方向上产生的库仑破裂应力变化,最后计算了前震和主震二者产生的库仑破裂应力变化,分析其对余震的影响,并着重讨论了震源破裂深度附近处于库仑破裂应力变化负值区内余震是否有被触发的可能。

1 漾濞地震对邻区的影响

1.1 漾濞地震的破裂模型

本文收集了Wang等[7]、张克亮等[9]和于书媛等[10]提供的破裂模型,考虑到后两者为大地测量方法反演的破裂模型均含有前震序列效应和余震序列效应,会影响到前震和主震对余震的触发结果,因此本文统一使用Wang等[7]采用主震的地震波形反演的破裂模型。

此破裂模型的走向为137.7°,与余震呈NW-SE向展布一致(图1)。他们将断层面分为3 km×3 km的子断层进行反演破裂分布,破裂面沿走向长33 km,沿倾向长21 km,图1中箭头方向和长短分别对应子断层的滑动方向及滑动量的大小,滑动主要集中在中部,以走滑为主。

图1 漾濞地震的破裂模型[7]Fig.1 The rupture model of Yangbi earthquake[7]

1.2 漾濞地震产生的库仑破裂应力变化对周围主断裂带的影响

本研究首先根据破裂模型得到位移u,结合Okada[8]解析表达式,计算出位移u对空间的偏导数,得到应变张量e,从而通过本构关系(视介质为Poisson体,两个拉梅常数均取3.0×1010Pa)得到应力变化张量σ,最后结合有限断层面参数,计算出正应力σn和剪应力τs,由于文献[6]已给出具体推导过程,所以在此不做赘述。由此可见在已知地震破裂模型的前提下,可以严格求解出弹性体内部产生的应变场[6,11-12]。考虑到活动断层上库仑破裂应力的变化受剪应力、正应力和摩擦系数共同影响,所以静态库仑应力变化可以表示为:

Δσf=τs+μ(Δσn-ΔP)

(1)

式中:σf为在接收断层上产生的库仑应力变化,如果该值不为0则表明此处受到前震影响,从而导致发震时间被提前或推迟。前人研究经验表明σf为正并超过0.01 MPa(阈值)时,对后续地震触发作用明显[1,13];τs为沿滑动方向断层面的剪切应力变化;σn为正应力变化,式(1)表明不论τs取值为正或负,当σn>0时,即正应力呈拉张状态,均利于断层滑动,反之则抑制;μ为摩擦系数;ΔP为孔隙压力变化,可用视摩擦系数μ′[μ′=μ(1-B),B为Skempton系数]表示孔隙压力减少摩擦系数的效应,μ′受地下介质影响,该系数常规取值范围为0.2～0.8[14],本文参考前人研究经验[15-18]取视摩擦系数μ′=0.4。

漾濞地震震中在维西—乔后断裂带附近,东侧发育有横向断裂较多的深大断裂——红河断裂带,维西—乔后断裂及红河断裂北端均为滇西北区最重要的活动断层之一[19-20]。为研究本次漾濞地震对周围断裂带的影响,本研究搜集了附近断层数据[21]以及主要断层的几何性质和滑动特征作为接收断层面,结合Wang等[7]给出的破裂模型以及各类型断层的解析式[8],计算出漾濞地震在震源深度(7 km)处产生的库仑破裂应力变化对周围断层的影响(图2)。

由图2可见:本次漾濞地震对维西—乔后断裂、澜沧江断裂带、红河断裂带等断裂带有一定影响。其中维西—乔后断裂中段库仑应力增加0.004 421 MPa,南端和北端减少0.017 66 MPa和0.002 736 MPa,故而该断裂南北端相对较稳定,并未有余震发生。澜沧江断裂南北两端库仑应力增加,中段减少,其中北端最大增加量为0.001 790 MPa,南端增加量较小,仅有0.000 219 MPa,均未达阈值,中段减少0.000 136 MPa。对邻近的深大断裂红河断裂及怒江断裂影响较小,对红河断裂库仑应力最大加载区也仅为0.001 762 MPa,对怒江断裂最大加载区仅为0.000 289 MPa,该结果与尹风玲等[19]对红河断裂带上库仑应力变化的时空演变进程相符。

通过对漾濞地震产生的库仑应力对周围断裂带的影响可知,本次MS6.4漾濞地震对周围断裂带影响较小,均未显示出明显的触发作用。值得关注的阈值只是前人经验给出的触发临界值,小的库仑应力变化也会影响该区的地震活动性[22],特别是库仑应力增加较多的维西—乔后断裂带的活动性一直很强[20]。所以应注意维西—乔后断裂及红河断裂带南端库仑应力变化较高区域的危险性。

1.3 漾濞地震产生的应力场在震源破裂处的水平面投影

为了能够更好地理解漾濞地震产生的应力场结果,本节首先计算了漾濞地震在地表产生的位移场。结果如图3(b),震中南北两侧物质向震中汇聚,且距震中越远位移量越小,说明区域内拉张作用明显,与垂直位移表现出的沉降较吻合。同理,东西两侧物质向外流出,距震中越远位移量越小,说明区域内挤压作用明显,符合东北侧垂直位移表现出的隆升。该位移场结果与张克亮等[9]通过高密度覆盖GNSS观测台阵连续观测给出的三维同震形变场以及于书媛等[10]利用Sentinel-1A卫星数据和D-InSAR技术获取的漾濞地震的同震形变场的结果趋于一致,这也能说明本研究使用Wang等[7]提供的破裂模型得出的结果是可靠的。

(将库仑应力变化表示于断层所在位置,填充色表示投影在对应地震破裂面上的库仑应力变化,红色对应于库仑应力加载区,蓝色表示库仑应力下降区,应力变化单位为MPa)图2 漾濞地震在7 km深度处产生的库仑破裂应力变化Fig.2 Coulomb failure stress change at 7-km depth caused by Yangbi earthquake

图3 漾濞地震产生的应力场在深度7 km处的水平面投影及地表位移场Fig.3 Horizontal projection of the stress field generated by the Yangbi earthquake at a depth of 7 km and the surface displacement field

计算的漾濞地震在震源破裂处产生的应力场在水平面的投影结果如图3(a)所示。本次漾濞地震震中附近面应力在0.1 MPa左右,震中东西侧产生较大面应力挤压,背离震中向东西两侧衰减,因紧挨维西—乔后断裂,至其南端时已减至0.02 MPa左右,到达红河断裂北端及怒江断裂南端时已锐减至数千帕,甚至更低;震中南北两侧产生较大水平面扩张应力,背离震中向南北两侧衰减,到达维西—乔后断裂中段和澜沧江北端时仅剩数千帕了,到达怒江断裂中段时仅有数百或数十帕。

在膨胀区内,各点以张应力(白色箭头)为主,压应力(黑色箭头)以震中为圆心的同心圆向外扩散;压缩区内,压应力为主,张应力以震中为圆心的同心圆向外扩散。整体上,主压应力和主张应力均呈条形磁铁磁场线分布状向外辐射。

从水平最大主应力(主张应力)和最小主应力(主压应力)方向来看,我们得到此次地震在震中区域附近变化的应力场结果(主压应力呈近EW向,主张应力呈近NS向)与Wan[23]得到的中国现代构造应力场的结果(NS向挤压,EW向拉张)恰好相反,这意味着此次地震产生的应力场在一定程度上抵消了该区域构造应力场,所以该地震是在区域构造应力场背景下的一次正常应变能释放。

2 漾濞地震序列间静态库仑应力影响分析

2.1 漾濞5.6级前震对主震的静态库仑应力影响

前震常被假定为发生在主震断层上的前兆滑动,因断层不规则加速导致主震发生,所以对主震的发生有着一定的预报意义[24-25]。此次漾濞地震为研究典型前震-主震序列提供了极佳的机会。本节重点讨论其间相互的静态库仑应力影响关系。首先采用震源机制中心解算法[26]分别计算了前震(节面Ⅰ走向306.83°,倾角81.27°,滑动角-166.98°;节面Ⅱ走向214.82°,倾角77.13°,滑动角-8.96°)和主震(节面Ⅰ走向45.45°,倾角84.25°,滑动角-2.85°;节面Ⅱ走向135.74°,倾角87.17°,滑动角-174.24°)的中心震源机制解,并根据余震展布将主震中心震源机制解的节面Ⅱ(走向135.74°,倾角87.17°,滑动角-174.24°)作为接收断层面,计算出前震在主震断层面和滑动方向上产生的剪切应力为0.12 MPa,正应力为0.039 MPa,根据式(1)可算出前震在主震断层面上产生的库仑破裂应力变化为0.13 MPa,表明前震引起的库仑破裂应力促进了主震的发生。

2.2 漾濞前震和主震对余震的静态库仑应力影响

除Wang等[7]提供的破裂模型外,还需要余震的精定位数据和对应的震源机制解。为此,本研究收集了漾濞地区5月18—28日的2 873个地震事件,使用Waldhauser等27]提出的双差定位法。参与反演计算的台站33个,P波权重1.0,S波权重0.5,8次迭代后垂直向定位误差缩小至0.837 km,水平东西向定位误差减至0.495 km,水平南北向误差减至0.671 km,得到2 115个精定位结果。其中最大前震MS5.6的定位结果为25.632 9°N,99.924 0°E;主震MS6.4的定位结果为25.689 1°N,99.873 8°E;最大余震MS5.2的定位结果为25.586 2°N,99.975 6°E,但考虑到MS≤3.0的中小余震有2 036个,占比96.26%,求解出的震源机制解缺乏可信度,加之余震持续发生,数据量庞大,所以本文无法一一给出漾濞地震余震序列的震源机制解。为此,本研究搜集段梦乔等[28]确定的漾濞地震震源机制序列,总体分析得到除主震之外地震序列的总体震源机制节面,发现与主震震源机制的最小三维空间旋转角仅为12°[29],所以可认定余震与主震的释放具有相同的形式,两者施密特等面积投影见图4。并结合前人[30-31]基于余震是主震破裂面延续的经验,假定余震接收断层面参数与主震破裂面一致的原则,给出余震的接收断层面参数。

本节并非求解单个地震事件产生的库仑破裂应力变化,所以先计算最大前震和主震的同震应力张量并叠加,再由同震应力场求解出库仑破裂应力变化。结合此次余震的平均深度在10 km附近,为能详细地计算前震主震对余震的触发作用,本研究计算4个深度上的库仑应力展布。结合Wang等[7]给出的破裂模型在深度为7 km附近滑动量最大,故将0～10 km处分为0～4 km和4～10 km两个深度范围区间在2 km及7 km处的库仑应力分布(图5)。最大余震的精定位深度为13.249 km,因此本研究将10～21 km的深度范围分为10～16 km和16～21 km分别计算出13 km处以及18.5 km处前震和主震的库仑应力分布对之后余震的静态库仑应力影响(图5)。为能更详细地分析余震展布情况,将分布于库仑破裂应力正值区或负值区的余震个数统计于表1,并给出库仑破裂应力正值区内余震的占比。

从图5和表1中发现,有54.52%的余震位于库仑应力变化正值区,说明一半以上的余震与主震的震源机制相似,余震破裂面与主震优势破裂面较为一致。在深度为2 km、13 km和18.5 km处,前震和主震产生的库仑应力变化为正区域内余震占比较高,在90%附近。在主震震源破裂处附近(7 km),破裂模型显示该深度滑动量虽然最大,但该深度处余震多分布在库仑应力变化负值区。本研究使用的破裂模型(图1)显示,滑动量主要集中在3～12 km深度范围内,其中最大滑动量在6～9 km的深度上,并且主震震源位于最大滑动量附近,主震的发生对周围的相同震源机制的节面上存在卸载作用,故而该深度抑制区余震事件明显多于加载区,在第3节将对其进行进一步讨论。

(红色到蓝色表示压缩到膨胀的过渡,P,T,B分别表示压轴、张轴和中间轴)图4 2021漾濞地震序列主震震源机制与前震和余震序列整体震源机制)施密特等面积投影图Fig.4 Schmidt equal area projection of the focal mechanism of main shock and the overall focal mechanism of the foreshock and aftershock sequences of the 2021 Yangbi earthquake

(图中三个沙滩球分别对应5.6级前震,6.4级主震,5.2级余震,绿色圆点对应余震事件,黑色框对应于破裂模型的断层面投影,红色区域为库仑破裂应力加载区,蓝色为库仑破裂应力抑制区)图5 前震主震分别在四个深度上产生的库仑破裂应力分布图Fig.5 Coulomb failure stress distribution generated by foreshock and main shock at four depths

表1 漾濞地震前震及主震对余震序列影响统计Table 1 Statistics of the influence of foreshock and main shock on the aftershock sequence of Yangbi earthquake

3 结论与讨论

本研究根据2021年漾濞地震的破裂模型[7]计算了均匀弹性半空间下漾濞地震产生的库仑破裂应力变化对周围断层的影响,以及对周围地区产生的应力变化在水平面上的投影,并结合2 115个余震精定位结果分析了“前震-主震-余震”序列间的静态库仑应力影响,得出以下几点认识:

(1) 漾濞地震对维西—乔后断裂南北端、红河北段以及澜沧江断裂中段有明显的应力卸载作用,对维西—乔后断裂中段、澜沧江北端和邻近的红河断裂带中段库仑破裂应力增加在千帕的数量级。考虑本次漾濞地震为一次震级仅为MS6.4的右旋走滑型地震,故而产生的库仑破裂应力对周围断层影响较小。尽管此次地震产生的库仑破裂应力对周围断层的影响未达阈值,但考虑到阈值也只是前人经验给出的参考触发临界值,因此应注意维西—乔后断裂中段库仑破裂应力变化较高区域的危险性。

(2) 主震震中附近水平面应力在0.1 MPa左右,东西侧产生较大水平面挤压应力,背离震中向东西两侧衰减,至维西—乔后断裂带南端时已减至-0.02 MPa左右,到达澜沧江断裂中段、红河断裂北端及怒江断裂南端时,已锐减至数千帕,到达怒江断裂中段时仅有数百或数十帕。漾濞地震产生的应力场在水平面上的结果显示,震中区域附近呈NS向拉张,EW向挤压的应力状态。与区域构造应力场方向相反,减弱了区域构造应力场,表明此次地震是构造应力场积累而导致的地震应力释放,是该区域应力积累后的一次正常释放。

(3) 以主震的节面II(走向135.74°,倾角87.17°,滑动角-174.24°)作为接收断层面,计算出前震沿主震断层面滑动方向上的剪切应力为0.12 MPa,正应力为0.039 MPa,前震在主震断层面上产生的库仑破裂应力变化为0.13 MPa,表明前震所引起的库仑破裂应力促进了主震的发生。

(4) 通过计算不同深度上前震和主震对余震的静态库仑破裂应力变化,在深度为2 km、13 km和18.5 km处,前震和主震产生的库仑应力变化为正区域内余震占比较高,可见前震及主震破裂引起的库仑破裂应力变化是该区余震持续不断的原因之一。

值得关注的是在破裂模型滑动量最大处(7 km)的触发比例极低。这些被抑制的余震均是在主震之后大量发生的,并且与触发区的余震在时空分布上并无差异。那究竟是什么孕育了这些余震,这些余震的发震构造又是什么样的,是否为主震破裂面的延续？推测在主震发生时,附近会产生大量的破碎带,破碎带中发育方位各异的薄弱面,在漾濞地震同震应力场作用下,部分薄弱面上静态库仑应力达到阈值便会触发余震,这样被触发的余震的发震方式会呈多样性分布[32]。为讨论在主震震源破裂处附近处于库仑破裂应力变化负值区内的大部分余震是否有被触发的可能,本研究将最大剪应力变化的方向作为断层滑动方向,将空间上的点最易错动的面作为接收断层面,并以此接收断层面计算的库仑破裂应力变化 即被称作该点的最大静态库仑破裂应力变化(ΔCFSmax)。

类比于式(1),万永革[33]曾考虑到Δτrake<0时剪应力变化的分量与滑动方向相反情况下,Δσn>0便会利于断层向滑动方向反方向滑动,也可以表述为抑制断层沿滑动方向滑动,Δσn<0抑制断层向反方向滑动,对断层沿滑动方向滑动并无促进作用,具体表达为式(2)。本研究将利用此公式去分析漾濞地震对其余震的最大静态库仑应力(ΔCFSmax)影响。

(2)

通过求解出的最大静态库仑应力变化(ΔCFSmax),发现余震均处于ΔCFSmax较大区域(图6),这就意味着该区域内满足被触发的易错动断层居多,在主震同震应力场作用下极大可能发生错动诱发余震。同时也表明在震源破裂附近余震发震断层面受破碎带影响严重,可见以主震破裂面作为接收断层面判定主震对余震的影响这一做法并不总是奏效。

图中沙滩球表示同震应力影响下,最大静态库仑破裂应力变化对应接收断层面(黑色弧线)的理论震源机制图6 震源破裂附近(7 km)处,最大库仑破裂应力变化及接收断层面随空间的变化Fig.6 Variation of the maximum Coulomb failure stress and the spatial variation of receiving fault plane near the source rupture(7 km)

本文在对周围断层的影响和地震序列间的触发关系均为一级近似估算。因此采用Wang等[7]提供的相对简单的弹性半空间破裂模型,实际的不均匀地球介质一定程度上会影响计算结果的准确性;震后地球的黏弹性松弛效应会导致应变扩散[34-35],但是这种黏弹性松弛效应对发震间隔小于30年地震的影响可以忽略不计[5]。

本研究在计算漾濞地震对周围断层的影响时,由于缺少后续余震的震源破裂信息,未曾考虑后续大量余震对周围断层的影响,但美国Landers地震余震能量释放仅为主震的4%[36],因此本研究认为余震相对于主震对周围断层的影响可以忽略。

尽管有上述不确定性,但本文在一阶近似下给出了该地震破裂对周围断层的影响,并估计了前震对主震以及前震和主震对余震的静态库仑应力影响。这些对于该地区的地震动力学研究是有意义的。

致谢:青藏高原研究所王卫民研究员为本研究提供了漾濞地震的破裂模型,两位匿名审稿专家和中国科学技术大学胡晓辉博士为本研究提供建议,本研究绘图采用GMT软件[37]绘制,特此致谢！