2022年1月8日青海门源M6.9地震余震精确定位与断层面拟合

赵 莉, 刘白云, 范 兵, 王文才, 张卫东, 刘艳云, 杜建清

(1. 甘肃省敦煌文物保护研究中心, 甘肃 敦煌 736200; 2. 中国铁塔股份有限公司甘肃省分公司, 甘肃 兰州 730000;3. 中国地震局地质研究所地震动力学国家重点实验室, 北京 100029; 4. 甘肃省地震局, 甘肃 兰州 730000;5. 新疆众志盛合建设工程有限公司, 新疆 乌鲁木齐 830011)

0 引言

据中国地震台网正式测定,2022年1月8日1时45分,在青海省海北州门源县(37.77°N,101.26°E)发生M6.9地震,震源深度10 km(https://www.cenc.ac.cn/cenc/dzxx/396391/index.html)。此次地震震中距门源县城54 km,距西宁市138 km。

震区位于青藏高原东北缘,是青藏高原往大陆俯冲的前缘部位,它的东北部连接阿拉善块体,东部毗邻鄂尔多斯块体,新构造较为强烈。新生代以来,印度板块持续北向挤压,造成青藏高原内部地壳缩短增厚,物质向周边挤出[1]。青藏高原东北缘由于受到具有坚硬中上地壳的阿拉善块体阻挡[2],发育了一系列地壳尺度的大型断裂带。祁连—海原断裂带是调节青藏高原和阿拉善地块汇聚的边界,托莱山断裂与冷龙岭断裂分列于该大型断裂之西、东,相互为右阶斜列衔接,是断裂带中西段的重要组成部分[3]。GPS及其相关研究表明,托莱山断裂与冷龙岭断裂滑动速率差异较大,在主震震中附近闭锁程度高,应力积累显著[4]。因此,在这两条断裂交汇处,变形不均性增大、弹性应变能累积,从而导致在该区域内地震活动较为活跃。此次地震是继2016年1月21日门源M6.4地震后该块体内发生的又一次强震。震源机制解显示为走滑型,距离最近的断层是托勒山北缘断裂。

震后兰州大学与甘肃省科技厅会同中国地震局兰州地震研究所、中山大学和青海省地震局等10余名科研人员当天赶赴震区,对本次地震地表破裂带、发震构造和震害等进行了现场考察。经考察后发现,本次M6.9地震地表破裂带长度超过22 km,主要由2条破裂带组成。其中,北侧主破裂带沿广义海原断裂带中段的冷龙岭断裂西段分布,从东向西大致经过硫磺沟脑分水岭、兰新高铁大梁隧道、道沟至下大圈沟西支沟止,整体走向NWW,长度超过17 km,向两端逐渐衰减;南侧的次级破裂带分布在广义海原断裂带中西段的托勒山断裂东段局部段上,大致沿大西沟、狮子崖一线分布,走向近EW,长度约5 km,构成与主破裂带西段左阶斜列的次级分支破裂,应为地震向西破裂的端部效应(1)资料来源于2022年1月18日科技日报(记者颉满斌).。

为进一步探明此次地震的发震构造及其几何形态,本文利用双差地震定位方法对2022年1月8日青海门源M6.9地震震后几天的余震进行相对精确的重新定位,然后根据震区小震线性分布的疏密程度并考虑地表破裂带的展布,选择其中余震相对密集分布的线性矩形区域,根据这些矩形区域内的余震震源信息,利用模拟退火及高斯-牛顿的计算机算法拟合本次地震断层面的参数。从而基于地震学的研究方法推断本次地震的发震构造及在深部的延展形态。

1 理论方法

1.1 双差精确定位方法

地震的活动性可揭示断裂带的几何形态结构[5],地震活动性与构造活动结合是推敲地震活动构造较为易于入手的方法,但这种方式的使用却依靠于地震定位精确度。双差算法原理是通过各台站计取的2个相邻震事件的观测和理论走时差的残差来判定它的相对位置,因此可抵消共同误差产生的不确定因素,尤其是与速度结构特征在横向变化及与接收有关的影响,不必进行对台站校正和精确估计走时就可大幅度改善其相对定位精确程度。目前这种方法已在地震构造研究中得到比较广泛的应用[6-8]。

在3D模型中,相对同样一个观测站点k,震源i和震源j两者之间的观测到时与理论到时之差的方程如下[9]:

(1)

在双差地震精确定位法中,在震源i和震源j距离较近时,可以认为震源到两个地震台站之间的传播路径是相同的,此时若走时速度结构是已知的,则式(1)可简化为:

(2)

(3)

WGm=Wd

(4)

式中:G为M×4N阶的偏导数矩阵式(M为观测双差数目,N为事件数量);m为4N阶震源参数矢量;W为加权的对角矩阵表达式;d为双差数据矢量表达式。

在反演计算过程中,对所有余震事件经再次定位得到的所有震源参数(即3个方向上坐标与发震时刻)加上了使其平均移动为0的约束条件:

(5)

在计算时,一般采用共轭梯度数学算法求解上述方程,可以得到带阻尼的最小二乘解,并利用奇异值分解方法得到有关模型参数的分辨度、残差等。

1.2 断层面与滑动角计算方法

一般来说中强震发生后小地震经常在断层面及其附近发生,这时利用一些小震位置参数就可以拟合得到断层面参数。其原理为:找寻一个平面,使得所有小震的震源位置到这个平面长度的平方和为最小。建立小震到所求断层面竖向距离与观测误差两者比值平方和为最小的目标函数:

(6)

式中:E为跟ρ、φ、δ三个变量相关的三元非线性函数表达式;n表示参与拟合小震个数;δ为拟合断层面的倾角值;σi是第i个小震定位的误差。反复搜索三个参数ρ、φ、δ的值,E最小时对应的φ值就为拟合断层面的走向。

尽管基于余震疏密情况可计算得到断层面参数,但滑动角对断层面之间作用的判断也是重要的。基于假定——断层方向和构造应力场在拟合断层面上相互作用时剪切应力的方向是一致的。万永革等[10]提供了断层面滑动角的一般计算方法:

设应力场P轴走向是φp、俯角是δp,T轴的走向是φT、俯角是δT,则在NED的地理坐标系当中方向矢量表达式为:

T=(cosφTcosδT,sinφTcosδT,sinδT)

(7)

P=(cosφpcosδp,sinφpcosδp,sinδp)

(8)

中间轴(B轴)的方向表达式如下:

B=P×T

(9)

主应力相对R大小计算过程如下:

(10)

式中:S1、S2、S3分别为最大主压应力、中间主压应力与最小主压应力。

在前述得到的断层面走向φ、倾角δ基础上,断层面滑移方向及法向在地理坐标系下可表达为:

l=(cosφ,sinφ,0)

(11)

m=(-sinφcosδ,cosφcosδ,sinδ)

(12)

n=(sinφsinδ,-cosφsinδ,cosδ)

(13)

式中:l、m、n分别为走向、向下及法向方向。在断层坐标系l、m和n中,走向和倾向的应力量表达式如下:

T1=cos(T,n)cos(T,l)S1+cos(B,n)cos(B,l)S2+
cos(P,n)cos(P,l)S3

(14)

T2=cos(T,n)cos(T,m)S1+cos(B,n)cos(B,m)S2+
cos(P,n)cos(P,m)S3

(15)

此时,滑动角为:

(16)

通常情况下,上述问题的解法有局部和全局算法,但都依赖于初始解,若初始解的选取不恰当,实际解会与正确解相差较远。万永革等[10]发展了模拟退火的全局搜索与高斯-牛顿局部搜索相结合的算法解决了上述问题,避免了对初始解的过度依赖,既可在全局搜索出最优解,还可估计出断层面参数的误差值。另外,在求得断层面参数的基础上进一步利用局部构造应力场参数可得到相关断层面的滑动角。因而被广泛应用到断层面参数的确定研究当中[11-17]。

2 数据

本文收集了青海及周边地震台网87个地震台站(其中甘肃地区7个测震台、青海地区80个测震台)于1月8—13日记录的余震到时数据。地震台站的分布见图1。在收集到的地震事件中,原始震中位于研究区(37.3°～38.0°N,100.7°～102.0°E)余震事件共867次。剔除单台记录及震级小于0.5的地震187次,最终选择的680次余震事件震级分布范围为ML0.5～5.3(图2),其中ML0.5～0.9地震1次,ML1.0～1.9地震427次,ML2.0～2.9地震201次,ML3.0～3.9地震41次,ML4.0～4.9地震8次,ML5.0～5.9地震2次,最大余震为1月8日ML5.3地震。被选地震事件总的震相到时包含了Pg震相到时为5 667条(图3),Sg震相到时为3 604条。

F1:阿尔金断裂带;F2:祁连山北缘断裂;F3:托勒山断裂;F4:冷龙岭—天景山断裂;F5:毛毛山—老虎山—海原断裂;F6:青海南山—循化南山断裂;F7:西秦岭北缘断裂;F8:东昆仑断裂;F9:玛多—甘德断裂;F10:巴彦喀喇山主峰断裂;F11:甘孜—玉树—风火山断裂;F12:黄河—灵武断裂图1 震区主要断裂及地震台站分布Fig.1 Spatial distribution of main faults and stations in the earthquake area

图2 研究区重新定位前的小地震震中分布图Fig.2 Epicenter distribution map of small earthquakes in the study area before relocation

图3 P波走时曲线图Fig.3 Travel time curve of P-wave

定位结果既涉及定位方法的选择,还与采用的初始速度模型有关。本次在双差地震定位时所使用的速度模型(图4)综合考虑了青藏高原东北缘地区人工激发地震和层析成像等研究[18-19],其每层详细信息列于表1,任意点的P波速度值用插值法提供。

表1 地震重新定位所用的速度模型Table 1 Velocity model for seismic relocation in this paper

图4 双差定位时使用的初始1D P波速度结构模型 (据文献[18-19]确定)Fig.4 1-D initial P-wave velocity model used for double diff- erence relocation (After reference[18-19])

3 结果分析

3.1 重定位结果

利用双差地震定位法对上述收集并整理好的地震观测资料进行再次定位。双差定位程序在读取震相到时时,由于不同震相分辨难易程度不同,因而读取精度也有差异,需对不同类型震相资料作加权。一般来说,初始定位时Pg到时读取质量比Sg高,因此对这两种类型震相各自赋予1与0.5的比重。考虑到计算范围内余震数量较多、地震台站分布较密的特点,台站和地震间距最大值选取为200 km、地震对之间间距最大值选取为10 km,每个地震最多同时可以与10个地震组成地震对。重新定位后丢弃了47个独立地震,得到了633个地震事件的精确震源位置。重新定位前后地震走时均方根残差平均值从0.389 s下降到0.016 9 s(图5)。余震水平方向的定位误差值在43～485 m间(图6)。对比重新定位前(图2)、后(图7)余震的震源位置分布,可以看出重新定位后的地震事件相对更加集中地分布于托勒山断裂东段与冷龙岭断裂西段的两侧。

图5 定位前后地震走时残差均方根直方图Fig.5 Histogram of RMS of seismic travel time residual before and after location

图6 重新定位后水平方向误差分布图Fig.6 Horizontal error distribution after relocation

F1:托勒山断裂;F2:昌马—俄博断裂;F3:祁连山北缘断裂;F4:冷龙岭断裂图7 研究区重新定位后的小地震震中分布图(粗青色矩形框为不同的反演区域段)Fig.7 Epicenter distribution map of small earthquakes in the study area after relocation (Thick cyan rectangular frame indicates different inversion region)

对震源深度深入研究可为确定地震孕育与发生的构造背景环境、区域变形及构造动力学性质等方面提供基础参数[20]。对比本次研究区地震事件在重新定位前后震源深度的频度分布特征(图8)可以看出,重新定位前本次地震的余震初始震源深度集中分布在5～15 km,重新定位后余震震源在0～20 km深度范围内偏正态分布,这与Fan等[21]对此次地震序列的震源深度分布特征研究结论是一致的。

图8 研究区小地震重新定位前、后的震源深度统计图Fig.8 Statistics of focal depths of small earthquakes in the research area before and after location

3.2 断层面参数反演及分析

根据震区活动断裂以及本次地震地表破裂带展布情况,并结合重新定位后余震沿主干断裂线性分布特征,本文沿托勒山断裂东段与冷龙岭断裂西段选定了两个矩形区域。利用上述两个矩形区域内的余震事件信息(图7)进行断层面参数拟合计算。综合模拟退火与高斯-牛顿的算法,求得的托勒山断裂东段与冷龙岭断裂西段断层面走向值、倾角值及拟合断层面4个顶点位置列于表2。

表2 拟合2个小段的断层走向、倾角、滑动角和位置Table 2 Strike, dip angle, slip angle, and position of fault by fitting two small segments

从拟合得到的断层面参数值及深部形态来看,两个断层面均比较陡立[图9(c),图10(c)],托勒山断裂东段断层面走向为95°,即近EW向展布,总长度约为15 km。冷龙岭断裂西段断层面走向为315°,即呈NNW向,总长度约为12 km。它们共同组成了本次地震的发震构造。

图9 托勒山断裂东段重新定位后小震的投影分布及距离断层面情况(圆圈表示再次定位后余震,粗线为拟合断层面边界;AA′为边界端点;SD是走向,DD是倾向,DF是余震与断层面的之间距离)Fig.9 Projection distribution of small earthquakes in the eastern part of Toleshan fault after relocation and their distance from the fault plane

图10 冷龙岭断层西段再次定位后的小地震分布情况,其余说明和图9相同Fig.10 Distribution of small earthquakes in the west section of Lenglongling fault after relocation, and other descriptions are the same as those in Fig.9

由于选择的两个矩形区域内余震分布相对比较均匀,因此拟合得到的断层面走向误差也比较小。另外被选择的地震事件基本沿断层走向方向对称分布,因而得到的余震与拟合的断层面距离分布频度[图9(d)、图10(d)]也大致对称,说明绝大多数余震是分布在所拟合断层面附近的,且基本以断层面为几何中心分别向左右两个方向相对对称分布。

3.3 滑动角计算结果与分析

在研究构造应力场与震源机制相互间关系时,不仅需要P轴走向φP和俯角δP,T轴走向φT和俯角δT,此外还需表示主应力相对大小的量R。卜玉菲等[22]将甘肃省及附近地区分割为1°×1°子区,利用这些地区M3.5以上地震震源机制方面资料反演平均应力场参数,得到了上述地区的应力场方向及相对大小的信息。本文从他们的研究结果中选取能代表本次门源地震震源区构造应力场参数和R值的资料(表3)。

表3 门源地区的局部应力场参数Table 3 Regional stress field of Menyuan zone

假设倾伏角误差值为10°,T轴和P轴方位误差值为5°,在前文拟合两段断层面走向、倾角及其标准差的基础上,采用表3中的局部应力场参数值及其标准差求解它们的滑动角,得到了每个分段断层面滑动角(表2)。从结果来看,托勒山断裂东段断层面的滑动角为2.2°,冷龙岭断裂西段断层面的滑动角为9.6°,表明它们的活动性质均为左旋走滑。

将拟合得到的断层面参数与震后考察的地表破裂带展布特征做对比,可以看出本次拟合的断层面与地表破裂带分布的位置基本一致。从拟合的每段断层面长度与其相应的地表破裂带长度来看,拟合的冷龙岭断裂西段断层面长度与对应的地表破裂带长度比较一致,拟合的托勒山断裂东段断层面长度比相应的地表破裂带长度长,对于两者相比延伸长度相差较大的原因还需要进一步分析确定。此外,本次拟合的两段断层面与现场地质考察的发震断层滑动性质相比也是一致的。

4 讨论

大震发震构造研究除传统方法外,还可以采用余震分布独立给出发震构造的几何参数与滑动性质。近些年来,随着地震测站的密度增加及观测技术的进步,地震观测震相到时拾取越来越精确。基于小震空间分布特征分析或者在此基础上通过某种反演方式来推断相关的地质构造等越来越引起地球科学家的关注。基于此次门源M6.9震区小震资料信息,综合高斯-牛顿局部搜索法与模拟退火全局搜索法拟合的断层面与地质考察吻合比较好,说明基于余震分布推测大震发震构造是可行的。

从余震重新定位前、后震源深度变化可以看出,重新定位前研究区内余震的震源深度集中分布于5～15 km,重新定位后10～15 km深度范围内的地震频度明显减少,总体上表现在0～20 km深度范围内偏正态分布,且具有震源深度较浅的特征。同时,拟合的断层面整体埋深深度也在上述深度范围内。这或许是由于青藏块体的快速抬升和在本区向NE推挤的过程中遇到阿拉善块体和鄂尔多斯块体的阻挡,在这一地区形成了弧形断裂组成的断裂束,而这些断裂束的新活动性强、埋深浅,故在活动过程中形成了一系列震源深度较浅的地震。

对于拟合的托勒山断裂东段断层面长度与对应的地表破裂带长度相比差距较大的原因,可能是因为拟合该断层面时地震事件的选取考虑了此段整体的余震分布,而通过地质考察的地表破裂带在该断裂段两端末端反映并不清晰。除文中开展的此类研究外,还有诸如震源区地壳速度结构与发震构造之间的关系等科学问题有待深入研究。

5 结论

文中利用双差定位法对本次门源M6.9地震的余震进行了重新定位,并基于重新定位后的震源信息,利用模拟退火与高斯-牛顿相结合的算法进行断层面拟合计算,获得的主要结论如下:

(1) 重新定位后大部分余震基本沿托勒山断裂东段与冷龙岭断裂西段分布,这可能意味着本次地震是上述两条断裂末端相互挤压共同破裂形成的。其挤压力源可能是在NEE向区域主压应力作用下,祁连块体在运动过程中受到了东侧鄂尔多斯和北侧阿拉善刚性块体的阻挡[23],其内部的冷龙岭断裂西段与托勒山断裂东段相连的尾端区域在长期的挤压下形成应力集中区。当应力积累到一定程度后集中释放时就会造成大震的发生。

重新定位后的余震分布与地表破裂带与范莉苹等[21]对本次余震重新定位等方面的研究结果是一致的。

(2) 重新定位前余震震源深度主要集中分布于5～15 km范围内,重新定位后变化为在0～20 km的深度范围偏正态分布。

(3) 拟合得到的托勒山断裂东段断层面走向为95°,即近EW向,其长度约为15 km;冷龙岭断裂西段断层面走向为315°,即呈NWW向,其长度约为12 km;从倾角来看,上述两断层面均为陡立状态,符合纯走滑断层的形态。

(4) 滑动角反演表明,两个断层面的运动性质均为左旋走滑。

致谢：本研究使用了万永革研究员提供的断层面计算程序,在此表示衷心感谢!