泸定地震诱发磨西断裂局部应力场偏转的离散元模拟分析

2023-01-31曾扬农苏占东孙进忠张建勇吴雪丽

防灾科技学院学报 2022年4期

曾扬农，苏占东，2，孙进忠，张建勇，唐磊，周剑，张涛，吴雪丽

（1.防灾科技学院地质工程学院，河北三河 065201；2.河北省地震灾害防御与风险评价重点实验室，河北三河 065201；3.中国地质大学（北京）工程技术学院，北京 100083；4.中国地震台网中心，北京 100045；5.北京工业大学城市与工程安全减灾教育部重点实验室，北京 100124）

0 引言

地应力是赋存于原位地质体内的应力，主要由构造应力和上覆岩层的自重应力构成［1，2］，它不仅是决定区域稳定性的重要因素，而且是各种地下或地面开挖岩土工程变形和破坏的主要驱动力，也是岩体行为预测、隧道围岩稳定评价和支护设计的基础数据。除此以外，地应力是固体地壳最重要的物理量之一，是影响和控制地震孕育和发生、断裂失稳活动的直接力源［2-5］。因此了解强震区和活动断裂附近的应力状态具有重要的科学意义，同时也是深入理解板块运动、构造变形以及地震活动等地球力学行为的重要基础［3，6］。

构造应力场的变化控制地震的孕育和发生，地震的发生势必造成地应力场的变化［7，8］。大量研究表明，地震的发生会使活动断裂附近应力场相对于区域构造应力场产生偏转。现场应力测试和分析发现地震前后断裂局部应力场会发生变化。Liao等［9］利用压磁应力传感器测量得到地震前后断裂附近岩体中同一地点的应力，结果表明昆仑地震后断裂附近地应力最大水平主应力方向发生了约20°的偏转；李方全等［10］对华北及郯庐断裂进行了原位地应力测量，发现经过大地震后，震中区域的主应力方向与区域构造应力场的方向有较大的偏转；针对震源处应力场的变化，何金等［11］采用同一个地震多个震源机制中心解的方法筛去重复地震事件的震源机制解对汶川地震进行剖析，发现西南到东北主压应力轴方向有所变化；孙业君等［12］以郯庐断裂带中南段及周边为研究区，采用研究区中小地震震源机制解求出日本9.0级地震使其应力结构及应力方向发生变化。此外，还有研究人员采用数值模拟的方法，对应力场偏转特征进行分析，Su等［13］运用PFC2D4.0软件模拟发现玉树地震引起甘孜-玉树断裂附近应力场偏转，且不同位置的局部应力场偏转规律存在显著差异。以上研究表明，无论是地壳浅表层的原位地应力测量和震源处的震源机制解反演结果，还是数值模拟和分析方法，都反映出了地震发生时断裂局部应力场会发生偏转。

由于断裂构造组成复杂，对于应力场偏转的机理认识还不够深入，还需要参照现有的同震监测数据对其进行多角度分析。为此，本文以2022年泸定地震时姑咱和天全两个地震台站的YRY-4型四分量钻孔应变仪监测数据为基础，运用离散元软件（PFC2D5.0）模拟泸定地震前后磨西断裂局部应力场的时空演化特征，为断裂局部应力场偏转机理的认识提供支撑。

1 研究区概况

本文选择2022年泸定地震发震断裂—磨西断裂为研究对象。磨西断裂为鲜水河断裂带南东段的一条主干断裂，位于鲜水河断裂带向安宁河与大凉山两条断裂转换的构造部位，具有变形强烈、性质复杂、多期活动的特点，一直被视为7级以上强震危险区域［14-16］。

鲜水河断裂带北起甘孜东谷附近，大体呈北西-南东向展布，全长约350km，是南北活动构造带的重要组成部分，具有雁列组合和断裂弯曲的几何特性，由炉霍断裂、道孚断裂、乾宁断裂、中谷断裂、雅拉河断裂、色拉哈-康定断裂、折多塘断裂、木格措南断裂和磨西断裂9条分支断裂组成，均为全新世活动断裂，并在石棉附近与龙门山断裂带和安宁河断裂带交汇构成了川西地区著名的“Y 字形”断裂带（图 1蓝色透明带）［17，18］。鲜水河断裂带是中国西南山区一条现今活动强烈的大型左旋走滑断裂带，具有规模大、活动性强、地震频度高等特点，自1700年以来已发生6级以上地震23次，其中7级以上地震8次，距今最近的一次强震为康定6.3级地震（图1红色实心圆）［14，19-21］。

图1 研究区构造背景Fig.1 Tectonic background of the study area

2022年9月5日四川省甘孜藏族自治州泸定县发生 MS6.8地震，震中位于 29.59°N，102.08°E，震源深度约为16km，最大烈度为Ⅸ（9）度［14，18，19，22，23］。

2 观测方法和数据

钻孔应变观测是研究地壳变形和地应力场变化的一种重要手段［24］。本次监测数据来源于YRY-4四分量钻孔应变仪，该应变仪分辨率可达10-10量级［25］，可以以每分钟或每小时一个数据点的采样率获得监测数据［26，27］。本文选取姑咱地震台和天全地震台在2022年泸定地震前后的钻孔应变监测数据，作为磨西断裂局部应力场变化分析的实际参考。

姑咱台站的应变测量传感器安装于40.69m深度处的基岩中，天全台站的应变测量传感器安装于39.5m深度处的基岩中，基岩岩性为花岗岩［26，28］。YRY-4四分量钻孔应变探头由4个依次呈45°角的元件组成（图2）［26］，以测量钻孔径向的应变变化。由于姑咱台站与天全台站分别自2006年11月1日、2020年1月1日开始运行，因此可以调查2022年9月5日泸定地震相关的应变（ε）变化（图3a、图3c），根据台站监测的应变数据，水平主应变方向θ由式（1）求出：

图2 四分量钻孔应变观测示意图Fig.2 Sketch diagram of four-component borehole strain observation

式中，εi（i＝1、2、3、4）是安装在钻孔中的四个元件的测量值；θ1为ε1对应传感器的方位角。

假设最大水平主应力（SH）方向与水平主应变方向一致（方向角均为θ），主应力方向的变化（偏转角变化值）记为Δθ。由图3b可以发现，从2006年11月到2022年9月姑咱台站Δθ值大约由-5°变为-30°，意味着姑咱台站钻孔中观察到的SH轴逆时针旋转了25°；由图3d可以发现，从2020年1月到2022年9月天全台站数据几乎没有变化。放大地震发生当天姑咱台站和天全台站的数据（分钟采样率），如图4所示，姑咱台站Δθ与原始状态相差约0.6°，天全台站Δθ变化很小，但是由于地震的发生，两个台站的数据都存在明显的阶跃现象，地震发生前后姑咱台站大约发生了1.5°偏转，天全台站大约发生了0.04°偏转。

图3 钻孔应变和偏转角变化值随时间的变化曲线Fig.3 Curves of drill hole str ain and deflection angle varied with monitoring date

图4 2022年9月5日偏转角变化值随时间的变化曲线Fig.4 Curves of the variation value of the deflection angle varied with time on Sep.5 in 2022

3 数值模型和方法

3.1 PFC 2D 5.0

PFC（Particle Flow Code）程序是Itasca公司开发的二维颗粒流离散单元程序，以牛顿第二定律及力-位移的物理理论为基础，通过离散单元方法来模拟圆形颗粒介质的运动和相互作用，分析大变形问题和颗粒细观结构的变化及其力学特性。模型中的每个颗粒表示为一个实体，颗粒组合在一起便可模拟变形多面体颗粒。由于只需定义颗粒与黏结性质，而不用定义整体的本构关系，其在岩体力学研究方面有很大的优势［29-32］。

3.2 磨西断裂模型

根据研究区地质模型，确定研究区几何模型（图5a）。如图5b所示，使用PFC2D5.0生成相应的磨西断裂数值模型，模型尺寸为120km×60km，区域空间由78444个半径120～200m的颗粒组成。为研究断裂应力场的空间变化特征，分别在断裂错动的前方（图5b，测量圆1和2）、断裂错动的后方（图5b，测量圆3和4）和断裂的端部（图5b，测量圆5和6）布置应力测量圆，并且在断裂的中部垂直剖面线上对称设置4个应力测量圆（图5b，测量圆7至10），在姑咱台站和天全台站对应位置分别布置应力测量圆（图5b，测量圆11和12）。

图5 磨西断裂模型示意图Fig.5 Diagram of the Moxi fault model

在PFC离散元模拟中，对于模拟岩土体材料，平行黏结模型和平直节理模型能更好地模拟其特性，但采用平直节理模型对岩体进行建模需要标定的细观参数数量增加，带来了更复杂的宏观-细观参数关系，导致标定困难［33］，因此本文选用平行黏结模型进行数值模拟。对于断裂（或节理）模型，一般是通过降低节理区域黏结强度或去除黏结实现，但是在节理滑动过程中颗粒相互阻挡，存在微观尺度的粗糙问题；光滑节理模型设置平行于节理的平面，允许颗粒相互重叠或通过［33］，故本文采用光滑节理模型模拟磨西断裂，模型参数如表1所示。

表1 数值模拟设置参数Tab.1 The setting of parameters of numerical simulation

假设模型剖面位于2km深处，应力大小根据杨树新等［39］提出的青藏地块水平构造应力进行计算：

式中，Z为深度（m）。

因此，最大水平主应力σT、最小水平主应力σt分别为45MPa和20MPa。最大主应力方向基于吴萍萍等［40］鲜水河断裂模型速度边界，并综合参考GPS速度场进行确定，最大主应力方向约为N49°W，在进行数值模拟时加载方向为水平方向，因此进行建模时对最大主应力方向与磨西断裂整体进行一定程度旋转，使最大主应力方向符合数值模拟中加载方向，如图5a所示。

4 结果与讨论

在PFC中，可以通过建立应力测量圆（图5b）来监测裂隙周围局部应力场的变化。应力测量圆记录的应力是测量圆中所有颗粒的应力均值。应力测量圆可以监测笛卡尔坐标系下所有的应力分量（σxx、σyy和σxy）［41］，运用式（4）对偏转角进行分析：

将测量圆11（图6a，模型上姑咱台站对应点位）与测量圆12（图6b，模型上天全台站对应点位）发生地震前后那一段的数据与图4中所示的现场观测相比，在地震发生前后，实测数据中姑咱台站的Δθ先增大后减小然后再增大，整体呈现减小趋势，应力测量圆11的Δθ在整体上同样呈现出减小的趋势，大约逆时针偏转了0.4°。在实测数据中姑咱台站Δθ突然增大，可能是因为断裂区域边界应力比增大，导致应力方向变化幅度增大［42］；而天全台站的Δθ在实测中变化范围较小，可能是该台站距离断裂较远，地震活动对其影响不大，实测数据所得的结果与测量圆12所得的结果都呈现上升趋势，但两者变化范围不在一个数量级，可能是因为模型模拟时间历程与实际观测时间历程之间的关系有所差别。两个测量圆变化情况与实际情况基本相似，表明数值模拟与实际情况基本吻合。泸定地震震后现场调查表明，断裂经过的地貌位置没有发生同震构造变形，除局部发育的小裂缝外，没有发现其他明显的地表变形现象［43］，表明本次地震未形成明显的地表错动，对浅表层的变形影响较小，因此整体而言最大主应力的方向偏转角度较小。而数值模拟与钻孔观测值不同，造成这种差异的原因可能是假设地震震源位于2km 的深度，而实际震源深度位于更深的位置，YRY-4四分量钻孔应变仪测量是在较浅的深度进行的。

图6 数值模拟台站应力测量圆的偏转角变化值随时间的变化曲线Fig.6 Curves of the variation value of the deflection angle of the stress measurement circle varied with time at the numerical simulation station

将同侧应力测量圆进行作图分析，根据万永革等［44］的初步研究结果，将应力测量圆进行分区，测量圆1和2位于压缩应力区，测量圆3和4位于拉伸应力区。将分居断裂中部测线两侧的应力测量圆的偏转角随时间变化过程分别成图，如图7所示。由图7可见，I区与Ⅱ区中三种不同的区域变化情况一致，压缩应力区与拉伸应力区以及断裂端部的应力场偏转角呈现不同的变化趋势：压缩应力区的Δθ整体下降，SH轴逆时针旋转；而拉伸应力区与断裂端部的Δθ则是整体上升，SH轴顺时针旋转。由于地震的发生，大约在70s后压缩应力区中测得的偏转角大约逆时针偏转了0.3°（图7，应力测量圆1和2）；而拉伸应力区中的最大主应力顺时针偏转了约0.2°（图7，应力测量圆3和4）。两个压缩应力区中应力测量圆的Δθ先是缓慢下降，随后急剧下降；在不同拉伸应力区中应力测量圆的Δθ先是几乎不变，然后呈现上升的趋势，Δθ随着时间的变化整体是增加的。Hardebeck等［45］指出，1992年兰德斯地震导致圣安德烈亚斯断层SH大幅旋转，其在拉伸应力区SH方位角顺时针旋转，这与本文研究结果一致。在断裂端部测量圆整体呈现上升趋势，说明断裂端部测量圆位于拉伸应力区，顺时针偏转了约0.3°。3个区域Δθ在70s后随时间的变化显示出一些波动性，主要是由于断裂错动导致的损伤扩展所致［46］。

图7 不同区域测量圆的偏转角变化值随时间的变化曲线Fig.7 Curves of the variation value of the deflection angle of the measured cir cles in different areas varied with time

图8为断裂中部垂直于断裂走向测线上各测量圆处的应力场偏转角随时间的变化过程。由图8可见，垂直于断裂剖面的近场测量圆（图7，测量圆7和8）和远场测量圆（图7，测量圆9和10）的Δθ变化趋势几乎相似，且整体都呈现上升趋势；断裂NE盘测量圆（图7，测量圆7和9）的Δθ变化比断裂SW 盘测量圆（图7，测量圆8和10）的Δθ变化大，这可能与断裂两盘岩体中局部应力场的强度差异有关。图9展示了磨西断裂中部垂直断裂走向测线上最大偏转角变化值（Δθmax）与测量圆到断裂距离之间的关系。由图9可见，同侧靠近断裂测量圆的Δθmax比远离断裂测量圆的Δθmax大；从监测数据来看，本次地震造成姑咱、天全台站偏转角偏转较小，但是依然呈现了靠近断层偏转大，远离断层偏转小的特点，这主要是断层的错动在断层周围引起了一个较大的变形的积累。Δθmax出现在断裂NE盘靠近断裂的测量圆7处，且NE盘远离断裂测量圆的Δθmax比SW 盘的大，说明NE盘局部应力场受到的扰动比SW 盘大。NE盘测量圆7到测量圆9，由近到远Δθmax由大到小；SW 盘测量圆8到测量圆10，由近到远Δθmax也是由大到小。两盘Δθmax由近到远，由大到小的空间变化规律相同，只是NE盘的Δθmax大于SW 盘的Δθmax。不论是模拟数据还是监测数据，其结论与Su等［13］对于Δθmax随着距离断裂的测量距离增加而减小的结论一致。

图8 断裂中部垂直于断裂走向上各测量圆的偏转角变化值随时间的变化曲线Fig.8 Curves of the variation value of the deflection angle of each measured circle in the middle of the fracture perpendicular to the its strike varied with time

图9 断裂中部垂直于断裂走向上最大偏转角变化值与测点到断裂距离的关系图Fig.9 The relationship of the variation value of the maximum deflection angle in the middle of the fault perpendicular to its strike with the distance from the measurement point to the fault

在本次模拟研究中，假设磨西断裂附近的应力变化主要与泸定地震有关，周围的其余断裂可能受泸定地震的影响而产生局部滑移，并可能影响钻孔周围的应力变化。因此模拟时仅局限于磨西断裂，该断裂局部应力场的变化代表了研究区域的主要局部应力场的变化，在后期的研究中将会引入多个断裂，考虑多个断裂的影响，并调查其余断裂对姑咱台站和天全台站同震应力分布和应力方向的影响。此外，此次模拟没有采用严格的边界条件进行实质性的破裂模型的探索，而是基于模拟结果与姑咱、天全台站监测数据吻合程度的基础上来分析断裂周围应力场的时空演化规律。同时由于应力场会受到研究区的工程地质条件、流体、温度等因素的影响，所以目前的分析仅限于定性的趋势预测，并没有深入的定量化分析。