2021年5月22日2时4分，青海果洛州玛多县(34.59°N，98.34°E)发生M_S7.4地震(图 1)，震源深度17 km(中国地震台网正式测定)^[1]。该地震极震区的烈度最大值为Ⅹ度，Ⅵ度及以上区域面积约为53 704 km²，宏观震中和微观震中分别为玛多县玛查理镇和黄河乡^[2]。从发震构造上看，青海玛多地震发生在青藏高原中部的Ⅱ级地块巴颜喀拉块体北部，距其北边界东昆仑断裂带约70 km。近年来，巴颜喀拉块体周缘发生一系列强震，如2008年汶川8.0级地震、2011年玉树7.1级地震、2013年芦山7.0级地震、2017年九寨沟7.0级地震，而玛多地震却是唯一发生在巴颜喀拉块体内部的强震^[3]，因此受到科研人员的广泛关注。人们从震源机制、震源破裂过程、发震断裂、余震序列、同震滑动分布、震区构造应力场、库仑应力变化和地震危险性等方面对玛多地震进行了研究。

中国地震台网、美国地质调查局等机构对此次地震的震源机制和震源破裂过程反演后认为该地震为一次左旋走滑型地震事件，兼有正断分量^[4]，为双侧破裂方式，震区构造应力场最大主应力轴方向为NE-SW，最小主应力轴方向为NNW-SSE^[1]。李志才等^[5]对GNSS数据进行反演确定的地震矩、断裂的几何参数及同震滑动分布，也揭示了发震断裂左旋走滑特性，同时得到同震滑动量的数值及分布。余震双差重定位研究发现，余震带与昆仑山口—江错断裂距离最近且两者走向有一定的重合度，说明玛多地震的发震断裂为昆仑山口—江错断裂^{[3, 6]}，与野外地震考察结果相一致^[7-10]。前人对于玛多地震同震及震后变形研究可以分为大地测量数据反演及同震、震后变形及相关应力变化模拟多个方面。研究者基于InSAR、GPS等大地测量数据与图像，揭示出最大滑动量为6.0 m^[11-12]，并确定了同震变形、静态滑动分布以及断层破裂过程^[13]，量化了地震间应变率、断层几何以及同震滑动分布^[14]。基于InSAR数据反演的断层滑动模型和PREM地球模型，研究者计算了此次地震的理论同震位移、大地水准面、重力及应变的变化，得到地壳变形的空间分布，为研究地震断层及解释大地测量数据提供了参考^[15]。在同震、震后变形模拟及库仑应力计算方面，杨光远等^[16]使用Coulomb3.3程序^[17-18]计算玛多地震引起的地震破裂面及周边主要断层上的库仑应力变化；谈洪波等^[19]和Li等^[20]使用均匀半空间解析解PSGRN/PSCRM程序^[21]对玛多地震同震、震后变形、重力变化及发震断层及周围区域的库仑应力变化进行了计算模拟。

前人使用解析法或半解析法对同震、震后及库仑应力变化进行理论计算，模型比较简单，计算速度快，因此，在大地震发生不久，通常会有相关的同震、震后及库仑应力变化成果发布，对大地震周围区域的地震危险性做出第一时间的评估和判断。发震断层周围介质的非均匀性、几何的复杂性，对大地震的同震变形影响很大，简单的解析解和半解析解很难考虑这些复杂因素，而数值模拟可以发挥独特的优势。在数值模拟中，有限元方法可以根据实际地质体的几何、材料及边界条件灵活地进行网格划分，同时进行并行计算，实现对复杂地质体高分辨率模拟。分裂节点法通过将位错加入到载荷向量而不改变刚度矩阵，可以简洁有效地将断层位错引入有限元计算。为此，本文采用C语言编写了分裂节点法的三维并行弹性有限元程序，用该程序计算了青海玛多地震的同震变形和应力变化，并对研究区域主要断层地震危险性进行探讨。

1 有限元方法 1.1 使用分裂节点法的有限元模型构建

用有限元方法模拟地震变形时，首先面临的问题是如何将断层的同震位错引入有限元数值计算中。分裂节点法可以简单有效地解决该问题^[22]。如图 2所示，单元1、2分别位于断层的左右两侧，U为位移，上标表示单元号，下标表示单元内的节点号；ΔU表示位错量^[22]。

则有

$ U_2^1=\overline{U}_2^1+\Delta U_2^1, $

(1)

$ U_1^2=\overline{U}_1^2+\Delta U_1^2 .$

(2)

单元1的节点2和单元2的节点1上的位移可以分为两部分，一部分是平均位移，即$ \overline{U}_2^1$、$ \overline{U}_1^2$这一部分位移对于2个单元来说是相同的，即

$\overline{U}_2^1=\overline{U}_1^2 .$

(3)

另一部分是位错量，即ΔU₂¹、ΔU₁²，满足

$\Delta U_2^1=-\Delta U_1^2 .$

(4)

单元1的单元刚度矩阵方程为

$\left[\begin{array}{ll} K_{11}^1 & K_{12}^1 \\ K_{21}^1 & K_{22}^1 \end{array}\right]\left[\begin{array}{l} U_1^1 \\ \overline{U}_2^1 \end{array}\right]=\left[\begin{array}{l} F_1^1-K_{12}^1 \Delta U_2^1 \\ F_2^1-K_{22}^1 \Delta U_2^1 \end{array}\right] .$

(5)

同样，单元2的单元刚度矩阵方程为

$\left[\begin{array}{ll} K_{11}^2 & K_{12}^2 \\ K_{21}^2 & K_{22}^2 \end{array}\right]\left[\begin{array}{l} \overline{U}_1^2 \\ U_2^2 \end{array}\right]=\left[\begin{array}{l} F_1^2-K_{12}^2 \Delta U_2^2 \\ F_2^2-K_{22}^2 \Delta U_2^2 \end{array}\right] .$

(6)

总体刚度矩阵方程为

$\begin{array}{c} {\left[\begin{array}{ccc} K_{11}^1 & K_{12}^1 & 0 \\ K_{21}^1 & K_{22}^1+K_{11}^2 & K_{12}^2 \\ 0 & K_{21}^2 & K_{22}^2 \end{array}\right]\left[\begin{array}{l} U_1 \\ U_2 \\ U_3 \end{array}\right]=} \\ {\left[\begin{array}{c} F_1-K_{12}^1 \Delta U_2^1 \\ F_2-K_{22}^1 \Delta U_2^1-K_{11}^2 \Delta U_1^2 \\ F_3-K_{21}^2 \Delta U_1^2 \end{array}\right] .} \end{array}$

(7)

由此可以看出，在分裂节点法中，断层位错的引入体现在单元、总体的载荷向量上，而对于单元、总体的刚度矩阵是没有影响的，大大简化了刚度矩阵和载荷向量的组装，且对于计算的准确性和稳定性没有影响^[19]，体现了该方法的简洁性和有效性。

1.2 有限元与理论解结果对比

为验证我们自主研发的采用分裂节点技术的三维并行弹性有限元程序的有效性，我们建立了一个三维走滑断层地震模型(模型1)和一个三维逆断层地震模型(模型2)，分别利用有限元程序和基于半无限空间位错理论^[23]的Coulomb 3.3^[17-18]及edgrn/edcmp^[21]计算2个断层模型的地表同震位移场，并对比两者的结果。2个断层模型的参数如表 1所示。

图 3分别展示了理论解和有限元计算得到的走滑断层和逆断层的同震地表位移的等值线图和剖面结果图。图 3(a)~3(f)左侧为Coulomb3.3^[17-18]解析解结果，右侧为有限元模拟结果，表明有限元得到的走滑断层地震和逆断层地震引起的地表同震位移的分布形态与理论解有很好的一致性。图 3(g)~3(i)分别给出了有限元和理论解的过断层中心点垂直断层走向的剖面上走滑断层地震的南北向同震位移、逆断层地震的东西向位移和垂直向位移的对比图，FEM表示有限元模拟结果，Theory1表示Coulomb3.3^[17-18]计算结果，Theory2表示edgrn/edcmp^[21]计算结果，表明有限元与理论解有很好的一致性。

2 青海玛多地震的同震变形和应力变化数值模拟与比较 2.1 有限元模型的设置

采用前人反演得到的非均匀位错模型以及我们基于分裂节点技术^[22]的有限元方法，本文计算了玛多地震的同震位移和同震应力场。有限元模型东西长1 100 km，南北长1 000 km，深度100 km，共有64个并行计算分区，有限元网格数量为2 171 715，可以实现较高分辨率的有限元并行计算。本文采用中国地震局地质研究所的非均匀位错有限断层模型^[15]，如图 4所示。

本文建立了2个有限元模型：均匀模型和分层模型，其材料参数如表 2所示，杨氏模量和泊松比可由表中的密度、P波速度和S波速度计算得到。由于有限元模型的尺寸大大超过主震断层，本文假定除模型上表面(地表)外，其余边界上的法向位移为0。

2.2 玛多地震地表同震变形和同震应力变化

根据上述三维弹性有限元模型，本文模拟得到的玛多地震地表同震位移场见图 5。东西向位移u大于零的部分位于断层南侧，小于零的部分主要集中于断层北侧，说明断层南侧以向东运动为主，断层北侧则以向西运动为主(图 5(a))；南北向位移v大于零的部分位于断层北侧，小于零的部分主要集中于断层东南部，说明断层南侧以向南运动为主，断层北侧以向北运动为主，呈现拉张的运动方式(图 5(b))。以上同震位移场分布特征与玛多地震具有明显的左旋走滑机制一致；垂向位移w大于零的部分位于断层北部，小于零的部分位于断层南部，断层北部为上升盘，南部为下降盘，上升盘为断层下盘，下降盘为上盘(图 5(c))，表明发震断层具有正断分量。

本文计算结果中玛多地震引起的研究区域内主要断层10 km深度的应力变化，如图 6所示。图中正应力变化大于零的区域表示拉应力增加，主要位于断层南西侧和北东侧，与拉应力减小的区域交错分布(图 6(a))；剪应力变化在发震断层两侧是小于零的，在断层尖端向外是大于零的(图 6(b))；库仑应力变化在断层两侧基本上是小于零的，库仑应力变化大于零的区域以断层两端为中心(图 6(c)和6(d)中的白色曲线分别是库仑应力变化的0和0.001 MPa的等值线)，在断层两侧呈蝴蝶状对称分布(图 6(c)、6(d))。9~11 km深度范围内除主震断层带两侧大于1.0级的余震基本位于10 km深度水平面上平行于主震的库仑应力变化大于零的区域(图 6(d))。

本文同时计算了研究区域内与玛多地震发震断层相邻的5条断层上的同震库仑应力变化，如图 7，计算使用的有效摩擦系数为0.4^[24]。结果表明，玛多地震引起鄂拉山断裂南段、昆中断裂东段、东昆仑断裂带的大部分(除中部靠西的小部分外)、玛多—甘德断裂北西段、南东段及主震震中附近以及达日断裂带北西段库仑应力增加，地震危险性增加，这些地方是需要密切关注、加强地震监测的区域。

3 讨论 3.1 玛多地震地表同震变形与大地测量观测数据及与前人计算结果的对比

为验证本文采用分裂节点有限元方法模拟实际大地震同震变形的有效性，本文设计了玛多地震的2种不同介质模型：均匀介质模型和分层介质模型，2种模型材料参数如表 2。本文分别计算了2种介质模型的同震变形场，并将其处理为卫星视线方向LOS位移量，与GPS观测的同震水平位移、InSAR升降轨LOS位移量观测结果进行对比，结果如图 8所示。图 8(a)、8(b)中的黑线表示玛多地震发震断层，蓝线表示有限元模型模拟结果，红线表示GPS观测结果(低频)，粉线表示GPS观测结果(高频)，左侧是均匀有限元模型模拟结果，右侧是分层有限元模型模拟结果。图 8(c)、8(e)、8(g)分别是InSAR升轨的观测数据、均匀和分层有限元模型模拟结果。图 8(d)、8(f)、8(h)分别是InSAR降轨的观测数据、均匀和分层有限元模型模拟结果，表明均匀有限元和分层有限元模型得到的升轨和降轨的卫星视线方向LOS位移量与InSAR观测得到的卫星视线方向LOS位移量有较好的一致性。

本文采用统计均方根(root mean square, RMS)方法对有限元模拟结果和观测变形量之间的差异进行定量说明。RMS^[15]定义如下

$\overset\frown{\sigma }=\sqrt{\sum\limits_{i=1}^n\left(u_{\mathrm{obs}}^i-u_{\mathrm{cal}}^i\right)^2 / n} .$

(8)

u_obsⁱ表示同震位移的GPS和InSAR观测数据；u_calⁱ表示同震位移的有限元模拟结果；n表示观测数据的数量。玛多地震同震位移的有限元模拟结果与观测位移量之间的RMS如表 3所示。

从表 3可以看出，有限元模拟同震位移的结果与大地测量观测结果之间的RMS在合理范围之内，与杨君妍等^[15]的计算结果一致，说明本文所采用的有限元方法对大地震同震变形的模拟是有效的。

本文通过分裂节点法将位错数据加入有限元模型计算的同震位移结果表明，发震断层南侧以向东、向南及向下运动为主，断层北侧则以向西、向北及向上运动为主，总体表现为有正断分量的左旋走滑断层。该计算结果与杨君妍等^[15]使用球形分层地球模型位错理论、谈洪波等^[19]使用均匀半空间位错理论进行同震位错计算的解析法或半解析法所得结果基本一致。

3.2 玛多地震对巴颜喀拉地块及周边地震危险性的影响

本文玛多地震同震库仑应力计算结果表明，发震断层的两侧沿走向方向上库仑应力下降，而在断层两端及以两端为中心呈蝴蝶状发散的区域内，库仑应力增加。这个现象表明玛多地震对发震断层及断层周边地壳产生了显著的应力调整：玛多地震部分释放了发震断层所积累的应力(发震断层两侧的库仑应力变化大部分都是显著降低)，也造成断层两端及周边部分区域的地壳应力积累(图 6(c))。该结果与杨光远等^[16]利用Coulomb 3.3程序在同等条件(10 km深度处和主震震源机制一致的微元面)下的计算结果对比，同震库仑应力分布基本一致。我们的模拟结果还表明，9~11 km深度范围内除主震断层带两侧大于1.0级的余震基本位于10 km深度水平面上平行于主震的库仑应力变化大于零的区域(图 6(d))。研究区域内鄂拉山断裂、昆中断裂、东昆仑断裂、玛多—甘德断裂以及达日断裂这5条主要断裂的同震库仑应力结果表明: 玛多地震引起鄂拉山断裂南段，昆中断裂东段，东昆仑断裂带的大部分、玛多—甘德断裂西北段和东南段以及达日断裂带北西段库仑应力增加(图 7)。鄂拉山断裂南段与昆中断裂东段位于巴颜喀拉块体北部的东昆仑—柴达木块体，历史上该地区曾发生过M5~M6及M6~M7地震，此次玛多地震引起该地区库仑应力增加，具有较高的地震危险性，应对其密切关注。东昆仑断裂带历史上曾发生过如1937年M_S7.5地震(阿拉克湖—托索湖段和下大武段)、1963年M_S7.0地震(秀沟—阿拉克湖段)以及2001年M_S7.1地震(库赛湖段)，是一条活跃断裂带^[25]。玛多地震使玛沁—玛曲段东南段和托索湖—下大武段东南段库仑应力增加，这与Li等^[20]的结果一致。这2个断裂段在玛多地震时产生应力积累，加之玛沁—玛曲段近年来未发生大地震，导致该地区地震危险性大大增加，应当加强地震监测、做好防震防灾工作。玛多—甘德断裂中段库仑应力变化分布较为复杂，库仑应力显著增加与显著降低的区段交错分布，可能与其复杂的断层结构有关^[4]。玛多—甘德断裂部分区段库仑应力显著降低可能与玛多地震释放了其部分应力有关，但其两端的库仑应力显著增加，地震危险性有所提高。位于玛多地震发震断层南部的达日断裂带北段及中段库仑应力增加，而其余部分库仑应力水平下降，考虑到1947年达日M7.7地震所带来的达日断裂带西北段库仑应力增加及其较高的应力积累速率^[26]，其地震危险性仍需进一步研究与评估。

4 结论

本文采用分裂节点技术开发三维并行弹性有限元程序，并通过地震位错模型的理论解验证了有限元程序的正确性。该有限元程序能够考虑复杂几何、非均匀材料和位错非均匀分布，可以模拟计算大地震同震变形。利用该程序模拟计算了玛多M_S7.4地震的同震位移场、应力场及库仑应力变化，计算得到的地表同震位移结果与GPS和InSAR大地测量数据一致。模拟结果确认了玛多M_S7.4地震的发震断层为一条左旋走滑兼有正断分量的断层。计算得到的库仑应力增加的区域与玛多地震的余震分布有较好的对应性。库仑应力结果表明，玛多地震引起了鄂拉山断裂南段，昆中断裂东段，东昆仑断裂带的大部分区段、玛多—甘德断裂西北段和东南段以及达日断裂带北西段库仑应力增加，这些区域的地震危险性增加，是未来需重点关注的地区。

石耀霖院士和周永发博士对本研究自主研发的三维有限元程序编制给予了很多建设性的指导和帮助，周元泽教授对论文初稿提出了很多建设性的意见，在此一并表示感谢。