邯郸地区主要断裂库仑破裂应力模拟与地震危险性分析

茅远哲 吕国军 孟立朋

1 河北省地震局，石家庄市槐中路262号，050021 2 河北省地震动力学重点实验室，河北省三河市学院街465号，065201

地震的孕育和发生是孕震体应力应变能不断积累、进入临界状态并最终失稳的力学过程，地壳岩石破裂过程中的应力分布、应力扰动与断裂构造明显相关[1]。因此，要探讨强震的迁移规律并开展地震危险性分析，最根本的途径是研究断裂带的应力状态及其动态演化过程[2]。

本文首先基于邯郸地区活动断裂地表、浅层、中层、深层探测和小震重新定位等地球物理探测结果以及前人研究成果，建立该区域三维粘弹性有限元模拟的断层模型；然后根据历史地震和古地震研究结果，加入运动学特征边界条件，模拟地震发生在目标区活动断裂上引起的应力应变场变化；最后计算1830年磁县地震引起的周围断裂面和滑动方向上产生的库仑破裂应力变化，给出研究区活动断裂地震危险性综合评价结果。

1 计算方法与模型

1.1 目标区内主要断裂分布

邯郸地区主要断裂呈北北东向和北西西向分布(图1，蓝色线段表示断裂，红色实心圆为磁县地震)，北北东向的断裂自西向东为紫山西断裂、太行山山前断裂和邯东断裂(邯郸隐伏断裂)，北西西向断裂自北向南为永年断裂和磁县断裂。这些断裂近似呈两两正交，控制着研究区内的构造活动。研究区主要活断层探测结果见表1。

图1 目标区主要构造背景

表1 邯郸地区主要断裂要素

1.2 几何模型的建立及介质参数确定

1.2.1 几何模型

根据地表和深部断裂展布资料，建立目标区几何模型(图2)。为防止边界出现突变，要求模型区域比目标区稍大，模型范围以紫山西断裂、太行山山前断裂、邯东断裂(邯郸隐伏断裂)、永年断裂和磁县断裂为界，每层划分为6个块体。因目标区内地壳厚度起伏不大，设计模型每层水平展布，模型各层厚度见表2。

图2 根据目标区断裂分布得到的平面几何模型

表2 模型中不同深度块体的横波波速

由浅中层地震勘探获得了主要断裂的倾角，但在地壳深部这些主要断裂的展布情况不明，因此在模型中将主要断裂均设为近直立断裂。由于目标区主要在模型中心部分，因此在模型中部对网格进行加密划分，分成10 563个节点、8 736个单元(图3，图中蓝色线段表示断裂)。

1.2.2 模型介质杨氏模量

利用有限元程序TEKTON进行模拟计算[3]，该程序使用的计算方法为根据岩石上施加的应力控制岩石破裂的机制，将莫尔圆实际半径(O1-O3)/2与破裂时半径之间的比值定义岩石极限破裂值Pf，可表示为：

式中，σ1与σ3分别为被施加的最大主应力与最小主应力。计算需要确定模型介质的杨氏模量和粘滞性系数等参数。

利用深部探测结果来计算模型内各不同深度块体的杨氏模量。由于宽频带地震台阵探测范围比模型区域要小，主要在东西向跨太行山山前断裂带和邯东断裂(邯郸隐伏断裂)、南北向跨永年断裂和磁县断裂的区域内布设；同时模型的主要关注区域在模型的中心部分。因此，用图3中A1、A2、A5、A6块体的波速代表整个块体的波速(有一定程度的简化)。根据宽频带台阵深部探测结果获取不同深度模型块体的横波波速，见表2。

图3 有限元模型结构和单元划分

地震波速与介质弹性常数的关系为：

(1)

E=2μ(1+ν)

(2)

式中，E为杨氏模量，VS为S波速度，ρ为岩石密度，μ为剪切模量，ν为泊松比。根据地壳岩石的常用量，ν设定为0.25；根据文献[4]，ρ取2.7×109g/m3。由式(1)可用VS和ρ求得μ，再根据式(2)求得E，由此得到有限元模型中不同深度块体的杨氏模量(表3)。

表3 有限元模型不同深度块体的杨氏模量

1.2.3 模型介质有效粘滞性系数

周永胜等[5]在总结近30 a来高温高压流变实验资料的基础上，应用流变数据结合地震震源深度分布，对华北地壳流变性质进行了研究，得到华北地区岩石圈强度剖面模型中岩石圈各层的流变强度数据。本文根据其研究成果得出30 km深度范围内的有效粘滞性系数(表4)。根据GPS观测的中国大陆地壳水平位移场计算出的最大剪应变速率值，目标区的剪切应变率在1.268×10-15/s左右，但对于模型中30～80 km深度范围，由于不确定其应变速率是否与地表一致，因此应变速率采用平均值1×10-15/s。

表4 有限元模型不同深度的有效粘滞性系数

1.2.4 模型边界条件

影响目标区动力学特征的主要因素为岩石圈介质特性和运动学特征，因此模拟时将目标区的运动学特征作为边界条件代入有限元模型进行计算。

假设模型底部垂直于地表方向位移固定、地表自由，由于模型只到80 km深度，根据GPS位移场计算出模型4个边界节点上的速度矢量(图4)，其中目标区内部断裂标注较多，其GPS点位所示方向与区外一致，呈SE向。

图4 华北地区GPS速度场

2 粘弹性有限元模型的建立和地震危险性分析

利用得到的断裂活动模型、历史强震数据以及区内岩石圈介质分布，可以建立粘弹性有限元模型，模拟地震发生在目标区内活动断裂上引起的应力应变场变化，计算活动断裂的断裂面和滑动方向上产生的库仑破裂应力变化，以分析目标区未来强震危险性。

1830年磁县断裂西段发生7.5级强震(表5)，利用断层有限元计算中的劈节点技术[3]，在模拟的弹性步中按表5数据施加在磁县断裂上。

表5 1830年磁县地震主要参数

研究区内紫山西断裂和磁县断裂中、西段现今时有小震发生，据此设计模型：磁县强震发生，并且紫山西断裂和磁县断裂中、西段小震不断。

1830年磁县地震距今已有190 a，而华北地区强震复发周期约在千年尺度，据此模拟的时间过程设计为：1)发震；2)震后几天的调整；3)震后100 a的时间过程；4)震后1 000 a的时间过程，设计模拟的时间步为：1 d×2+2 a×2+10 a×10+100 a×3+1 000 a×3。

参照文献[6-7]的方法，取μ′= 0.4。数值实验表明，改变此值对计算得到的库仑破裂应力变化的空间分布影响不大，但对其大小有一定的影响。将剪切应力变化投影到假定的滑动方向，与假定的滑动方向一致取正，反之取负。假定的滑动方向取自后续破裂事件的震源机制。

3 模拟结果与讨论

根据以上模型设计，得到模型的模拟结果。首先根据时间尺度的不同，选取0 a、4 a、204 a、3 404 a的时间节点，提取时间步t0(弹性步)、t4(4 a)、t15(204 a)和t20(3 404 a)的位移场进行分析(图5、6)。可以看出，磁县地震发生时(图5)，由于存在水平和垂向位错，使得震中附近的水平位移场发生扭转，相对于远离震中的位置，目标区内震中的位移量明显大于其他地区。

图5 磁县地震时和震后4 a模型的位移场

震后4 a，目标区内不同位置位移量的差异减小，但区内位移方向分布没有明显变化，说明强震引起的局部位移异常虽然在慢慢减弱，但是还未消除。震后204 a和3 404 a的位移图像(图6)表明，整个区域内位移比较均匀，已经看不到局部位移方向和大小的异常，说明磁县地震的影响已经被完全吸收。

图6 磁县地震震后204 a和3 404 a模型的位移场

图7(单位MPa)为1830年磁县地震引起的周围断裂面和滑动方向上的库仑破裂应力变化。图中蓝色至绿色表示断裂上库仑应力降低，地震危险性降低；红色至黄色表示断裂上库仑应力增加，地震危险性增加。由于紫山西断裂和磁县断裂中、西段深部小震沿断裂分布[8]，确定为2条现今活动断裂，其库仑应力的变化对于地震危险性评价至关重要。模型显示，磁县强震发生后，磁县断裂中、西段库仑应力加强，在震中附近位置库仑应力降低。对于紫山西断裂，在断裂南端，特别是与磁县断裂交接处库仑应力加强；在断裂中段，库仑应力减弱；在断裂中偏北的一小段断裂上，库仑应力增强明显；在断裂北段，库仑应力增量值较小。说明在磁县断裂中、西段两端与2条北东向断裂交接的部位，地震危险性值得关注。

图7 磁县地震在周围活动断裂的断裂面和滑动方向上产生的库仑破裂应力变化

4 结 语

研究区内，磁县断裂中、西段两端与2条北东向断裂(紫山西断裂和太行山山前断裂)交接处库仑应力有所加强，未来地震危险性值得关注。

致谢：感谢中国地震局地质研究所陶玮副研究员为本研究提供巨大帮助。