沈千贺, 万永革,2
(1. 防灾科技学院, 河北 三河 065201; 2. 河北省地震动力学重点实验室, 河北 三河 065201)
2021年5月22日北京时间02:04:12青海省果洛州玛多县发生7.4级地震,地震位于青藏高原北部、巴颜喀拉块体内部,是青藏高原近期发生的一次较大地震。巴颜喀拉块体一直以来地震都比较活跃,以往地震主要发生在巴颜喀拉块体周缘,而玛多7.4级地震(34.56°N,98.34°E)发生在巴颜喀拉块体内部[1]。王未来等[2]认为该地震发震断裂位于玛多-甘德断裂与达日断裂之间,江错断裂附近。潘家伟等[3]进行野外实地考察,发现此次地震形成了一系列由张裂隙、张剪裂隙、剪切裂隙、挤压鼓包和裂陷等多类型破裂雁行状组合而成的复杂同震地表变形带,总体表现为左行走滑运动性质。该破裂带主要沿东昆仑断裂带南部的江错断裂分布,整体呈N105°E走向,全长151 km。这次地震的发震断裂为江错断裂,该断裂向西延伸可与2001年东昆仑8.1级地震的发震断层昆仑山口断裂相接。为了进一步了解玛多地震的发震背景,本研究采用玛多地震的震源机制聚类分析断层面。
2021年玛多地震发生后,各机构测定到大量的余震震源机制数据,本研究对玛多7.4级地震以及其余震序列进行聚类分析,得到玛多地震发震断层的基本产状,再通过应力场反演计算断层的滑动角,最终确定断层的几何信息。
利用震源机制解,通过聚类方法与应力场反演方法来确定断裂带产状的优势在于相较于直接测定断层面产状(钻孔、探槽地质资料测定产状)或重力资料测定产状,本研究可以降低地质环境等因素造成的干扰,并且可以得到无法直接测量的地壳深部断层形态。一般大型断裂带上地震频发,大震后余震较多,且有较多手段可以测定震源机制解,如P波初动求解[4]、CAP方法[5]、Snoke方法[6]、FOCMEC方法[7]等。因此利用震源机制聚类分析发震断层的几何形态,在数据获取方面十分方便。
玛多7.4级地震发生在周缘地震频发的巴颜喀拉块体内部,远离地震多发的边界区域[8],所以此次地震对地震动力学和地震预测有重要的实际意义。本研究旨在通过玛多地震序列震源机制聚类分析断层的几何形态,通过应力场分析其滑动性质,并分析该断层发生的应力背景,为研究玛多地震的发震构造和孕育过程研究提供资料。
不同机构以及作者给出的震源机制解不同,会对后续的聚类分析和应力场计算中数据的选择产生影响,所以在进行聚类前需要对同一地震不同震源机制解进行中心解的求解,以获得更为精准的数据。本文采用万永革[9]的方法来求解中心解,首先选择一个震源机制数据作为初始解,经过一阶泰勒展开后通过Levenberg-Marquardt 方法进行迭代求解,当解的改变量非常小或达到最大迭代次数时,终止迭代。所求中心解需要使所有的震源机制解与中心解的空间旋转角平方和最小。
下面以主震为例展示震源机制中心解的求解过程。我们搜集到8个机构和作者的震源机制解(赵韬等[1]、中国地震台网中心[10]、中国地震局地球物理研究所[11]、美国地质调查局[12]、德国波茨坦地学研究中心[13]、全球矩心矩张量数据中心[14]、张喆等[15]、张建勇等[16]),分别以各个机构震源机制为初始解,最终得到的最优中心解分别为:走向102.35°,倾向88.34°,滑动角-6,32°,标准差为19.74;中心震源机制的第一个节面走向、倾角和滑动角分别为:102.35°,88.34°,-6.32°;第二个节面走向、倾角和滑动角分别为:192.53°,83.69°,-178.33°。按照震源机制类型的划分[17],属于走滑型地震。震源机制与中心震源机制的最小三维空间旋转角为5.03°;中心震源机制的P轴走向和倾伏角分别为:57.28°,5.64°;B轴走向和倾伏角分别为:267.66°,83.47°;T轴走向和倾伏角分别为:147.60°,3.28°。同时得到该地震的海滩球图[图1(a)]和震源机制辐射图[图1(b)]。图1(a)中黑色弧线为中心震源机制的两个节面,绿色弧线区域为其不确定范围;红色的点为震源机制解的P轴,黑色的点为震源机制解的T轴,黄色的点表示震源机制解的B轴,其周围对应颜色的封闭曲线表示其不确定范围;绿色、黑色和蓝绿色的点为各个机构得到的震源机制解P轴、T轴和B轴;紫色弧线为各个机构和作者得到的震源机制节面。图1(b)中压缩象限为蓝色,膨胀象限为红色。
图1 玛多MS7.4地震的震源机制中心解及震源机制的辐射图型
按照对主震震源机制处理的相同策略,对所有搜集到多个震源机制解(包含此次地震的余震以及该地区的历史地震)的地震进行震源机制中心解求解,共得到含主震在内的137个震源机制结果,其中有4个地震远离余震破裂带,明显不属于该地震序列,删除这4个震源机制,将其余震源机制解绘制于图2,并将除主震外的132个震源机制数据列于附录1中。玛多地震序列震源机制有走滑型65个、逆断型20个、逆走滑型15个,其余为不确定型。可以看到这些地震沿NWW和SEE方向分布,余震分布在巴颜喀拉块体内部、玛多-甘德断裂和达日断裂中的江错断裂上,验证了王未来等[2]研究推断的发震断层为昆仑山口-江错断裂。
图2 本研究所用震源机制解(左下角的图片展示了断层所处位置,图中红线为发震断层所处的位置[2],黑线为其他断层)
本小节确定了发震断层的大概位置,得出了更为精确的震源机制解,为得到具体的发震断层几何数据,需要采用震源机制节面聚类方法对这些没有争议的震源机制节面进行聚类。
本节对之前筛选整理得到的更为精准的震源机制解节面进行基于密度的聚类分析,其主要原理为只要某一区域对象的密度大于某一阈值时,就将其加到相邻的聚类中,对于簇中每个对象,在给定半径ε的邻域中至少要包含最小数目k个对象。相比其他的聚类方法,基于密度的聚类方法(DBSCAN算法)可以在有噪音的数据中发现各种形状和各种大小的簇[18]。该算法基于一组邻域来描述样本集的紧密程度,它将簇定义为密度相连点的最大集合,能够把具有足够高密度的区域划分为簇,并可在有“噪声”的空间数据库中发现任意形状的聚类[19-21]。基于密度聚类的优势在于可以发现任意形状的簇,并且对噪声不敏感。基于密度的聚类方法需要2个参数,分别为密度邻域半径ε和ε邻域包含的对象最小数目k[22]。根据万永革[23]提出的断裂带震源机制节面聚类确定断裂带产状方法及其在2021年漾濞地震序列中的应用,最小数目k经验设定为:
k=int(m/25)
(1)
式中:m为需要分类的数据总个数。
聚类对象的邻域半径ε表达为:
(2)
通过对玛多地震的震源机制节面进行聚类分析可以得到两簇聚类中心,根据结果绘制出玛多地震震源机制聚类结果图(图3)。其中震源机制节面采用绿色弧线表示,聚类中心节面采用红色弧线表示;黑点表示震源机制节面的极点位置,红点表示聚类中心的极点位置。聚类中心极点周围的蓝绿色椭圆为聚类中心的置信区间。第一簇[图3(a)]节面数为65,得到的中心节面法向的方位角为 23.49°,倾伏角为 1.79°,标准差为20.00°,其中心节面的断层面走向为113.49°,不确定范围是100.58°~126.41°,倾角为88.21°,不确定范围为82.24°~85.82°;第二簇[图3(b)]节面数为56,得到的中心节面法向的方位角为104.84°、倾伏角为 7.14°、标准差为17.95°,其中心节面的断层面走向为 194.84°,不确定范围为 186.68°~203.01°,倾角为82.86°,不确定范围为74.73°~89.00°。第一类和第二类的中心节面夹角为81.19°,由此可以得到第一类中心节面与第二类中心节面接近垂直。噪声节面[图3(c)]数据个数为145个。
图3 2021年玛多地震震源机制解聚类结果
将得到的两个聚类节面与主震求解的中心解(走向102.35°,倾向88.34°,滑动角-6,32°)和潘家伟等[3]野外实地考察的结果(走向为N105°E,全长151 km)进行对比,发现第一类中心节面与两者具有较高的一致性,因此本研究选择第一类节面作为求解玛多地震序列发震断层的结果。聚类得出的结果与王未来等[2]和潘佳伟等[3]的研究结果相差无几,证明了本研究的准确性。
聚类方法只能确定断层的走向和倾角,无法得出断层的滑动角,震后的地质考察[3]也仅给出了断裂的走向。欲估计滑动角,需要求解该断裂带的应力场,根据应力场在该断层上投影的剪应力方向确定该断裂的滑动角。
断层的滑动特征是发震断层的重要研究内容。上一节通过对震源机制进行聚类得到了两个节面的走向和倾角,而想要得到断层的滑动性质,则需要先求解发震断层的应力场。本节采用多个地震的震源机制解确定应力场的方向和大小的方法[23]来求解应力场。本研究采用万永革[24]设计的程序,通过网格搜索法对断层进行应力场分析。
应力因子R值代表三个主应力的比例关系,应力因子R为[25]:
(3)
式中:σ1代表主压应力大小;σ2代表中间应力大小;σ3代表主张应力大小。
应力场分析采用的准则为滑动方向与剪切应力的方向最为一致,求解结束后采用F检验得到参数的估计空间[24]。经求解运算得到玛多地区断层的应力因子R=0.9,最优应力状态的剪切力最大的两个节面的走向、倾角、滑动角分别为 280.0°、75.9°、15.7°和186.0°、74.8°、165.3°。最优的压轴走向、倾伏角分别为52.87°、0.72°;中间轴走向、倾伏角分别为321.00°、69.00°;张轴走向、倾伏角分别为143.14°、20.99°;应力比值为0.90;滑动角差别的平方和为294 986.29(图4)。
[图(a)中绿色曲面为90%置信度下应力场的最大剪应力节面;大蓝箭头为张轴方向,大红箭头为压轴方向,相同颜色的封闭曲线表示其不确定性范围;小红箭头为理论滑动方向,小蓝箭头为观测滑动方向;应力因子R=0.9。图(b)红色代表压缩轴,蓝色代表拉伸轴]
图4中可以看出断层受到NE-SW向的挤压和NW-SE向拉张的作用,且应力因子R=0.9接近1,表明挤压应力占突出优势,而拉张呈现出各个方向的弥散情况。由此推测玛多地震是由于南部的印度板块向北挤压亚欧板块并在东北受到阻挡,产生的应力经过长期积累在巴颜喀拉块体内部发生破裂造成的。与张建勇等[16]研究得出的玛多地震可能是东北向挤压的青藏高原遇到四川盆地、华北地块等坚硬块体的阻挡,使得板块内部形成与板块边界相似的断裂,从而发生强震活动的结论相一致。
为得到断层的滑动性质,本研究根据万永革[23]提供的方法和上节反演得到的应力张量参数计算得出,在玛多主震断层面(走向102.35°,倾角88.34°)上的相对剪应力和相对正应力分别为0.97和-0.43,聚类得出的发震断层面(走向113.49°,倾角88.21°)上的相对剪应力和相对正应力分别为0.84和-0.79。图5展示了区域应力场下研究区域各个走向和倾角的断层面上的相对剪应力和相对正应力分布。图中可以看出该断层面参数空间存在四个剪应力大的区域(走向10°,倾角90°与走向360°,倾角90°是同一区域),分别为走向10°,倾角90°;走向100°,倾角90°;走向190°,倾角90°;走向280°,倾角90°的集中区。上述四组节面极易在剪应力作用下产生破裂,相对正应力存在两个最大的区域分别为走向150°,倾角90°和走向330°,倾角90°的两个集中区。玛多7.4级地震及其发震断层正处于剪应力大的区域上,剪应力强度非常高,易形成剪切作用下走滑性质的运动。
(底图颜色表示相对应力值的大小;震源机制分类按照万永革 [17]进行:NS为正走滑型,SS为走滑型,NF为正断型,TS为逆走滑型,TF为逆断型)
通过应力场信息拟合断层面,计算得到聚类断层的滑动角为-0.72°[26]。综上所述,最终得出江错断层走向为113.49°,倾角为88.21°,滑动角为-0.72°,表明江错断层是在剪切作用下的左旋走滑断层。该结论与玛多地震主震性质一致,而且玛多地震序列震源机制以走滑型为主,与断层运动性质相符,证明了断层滑动性质结论的准确性。但需要说明的是,主震的震源机制跟地震序列的整体破裂机制通常是不一样的.地震序列整体震源机制所反映的是整条断裂的破裂性质,而主震只是反映1次地震的情况。由于主震震级比较大,所以主震的震源机制和整体破裂所表现的震源机制相当接近(空间旋转角为13.3°),本次研究就说明了这一点。
本研究的理论根据是发生在某一断裂带上地震的震源机制节面可能与该断裂带几何形状一致[22]。本研究主要通过对玛多震源机制聚类分析来确定地震的断层面。
研究通过对玛多地震序列震源机制解的中心解求解和筛选,得到更为精准的震源机制解,并对这些震源机制解进行基于密度的DBSCAN法的聚类,确定断层的走向和倾角后进行应力场的反演,利用得到的应力场数据计算获得滑动角。本研究得到以下结论:
(1) 通过对多震源机制进行聚类方法的研究,得出与王未来等[2]、潘家伟等[3]的研究具有极高一致性的断层参数,即走向为113.5°,倾向为88.2°。同时证明了采用震源机制节面聚类求取断层面形状的准确性。
(2) 通过对发震断层进行应力场分析,得出玛多地震发震断层处于NE-SW挤压和NW-SE拉张的应力状态,且断层的滑动角为-0.72°,与玛多地震序列震源机制多为走滑型地震相符。表明玛多地震发震断层是受到NW-SE方向的挤压和NW-SE方向的拉张,在江错断裂上形成了较大的剪应力,致使江错断裂表现为左旋走滑断层错动。研究认为玛多地震是由于印度板块向北挤压亚欧板块,在巴颜喀拉块体内部受到阻挡,致使近东西走向的江错断裂产生破裂所致。
本研究统计了包括主震在内的133个震源机制(132个余震),其中有走滑型65个、逆断型20个、逆走滑型15个,其余为不确定型,与赵韬等[1]和张建勇等[16]研究得出的玛多主震为左旋走滑型地震的结论一致。通过聚类和应力场分析求解的断层信息为走向113.49°,倾角88.21°,滑动角-0.72°;潘家伟等[3]实地考察的结果为发震断裂沿N105° E走向的断裂分布,主要表现为左旋走滑断层;主震的中心解与求解的发震断层结果相近,差距较小,同时与实地考察结果相近,证明了本研究的准确性。
本研究采用对震源机制解进行聚类分析方法和应力场反演方法来确定断层面的几何数据,与实地考察和前人研究结果有很高的一致性,证明了聚类方法的可行性,同时为玛多地震发震构造研究提供了资料。由于一般发生在断裂带上的地震与断裂带本身几何形态相似,所以对大量的震源机制进行聚类可以得到一个一致性较高的集合,该集合节面中心的形状近似于断裂带的形状。此方法适用于拥有较多震源机制数据的地震。
致谢:本文采用GMT软件[27]绘制而成,特此致谢。