山西南部区域应变场演化与地震关系的数值模拟

李宏伟，刘瑞春，陈永前

(1.山西省地震局，山西 太原 030021；2.太原大陆裂谷动力学国家野外科学观测研究站，山西 太原 030025)

0 引言

山西中南部至晋陕交界地区历史上曾发生过多次强震，2014年先后又出现了多项准同步前兆异常，部分异常测项于2015年下半年又同步转折变化，且与山西南部GPS区域应变参数转折变化较为一致。从区域应变特征来看，最小主应变于2014年5月出现偏离背景状态的加速负值变化，表明区域压性应变在持续增强，与基线的缩短状态也较为吻合。在此背景下，最大主应变和第一剪应变在2015年5月至12月同时出现短期增大-恢复的变化过程，之后发生了2016年3月12日运城盐湖ML4.8地震；另外，ML≥3.0地震视应力空间图像显示山西南部区域高值异常仍在持续[1-2]。因此需要开展基于地球物理场观测的具有物理意义的地震预测工作，为山西南部地震趋势判断提供科学依据。

地震的孕育是区域应力积累到逐渐加强的过程，当断裂积累的应变能超过断裂所能承受的临界阈值时，就会以地震的形式突然释放应变能。研究者在实践中也不断验证了这一结论：地震的发生与区域应变能的积累程度有关[3-6]。近年来，随着GPS观测精度以及测点密度的提高，逐渐被应用到数值模型的边界约束条件中。从1999年到现在，先后有刘杰等[7]、冯向东等[8]、刘峡等[9-11]对华北地区进行了三维有限元数值模拟。其中，2007年刘峡[9]利用研究区域的GPS速度场作为模型的边界约束条件，分析了模型内部的GPS速度场与模拟结果的一致性，反复调整模型材料参数，使两者尽可能误差最小，最终得到该区域的应力场特征。2013年刘峡等[11]采用“块体加载”的思想，模拟了汶川地震前后山西断陷带的应力场空间格局，结果显示，汶川地震后，整个山西断陷带呈现构造挤压增强态势。随后，山西断陷带于2010年相继发生河津4.8级地震、大同4.6级地震以及阳曲4.6级地震。2017年李平恩等[12]利用三维粘弹性模型模拟了太原盆地应力场演化与强震的关系，表明盆地的地震活动性受区域应力场控制。目前，山西断陷带的地震视应力计算结果显示，山西断陷带南部存在高视应力区[1-2]，因此，该文利用ANSYS有限元软件，结合研究区地质构造以及断层几何分布建立山西南部三维有限元模型，以2011-2016年GPS速度场作为边界加载条件，分析山西南部近期应变场空间分布特征。

1 山西断陷带南部地质构造背景及模型的构建

1.1 研究区地质构造背景

山西断陷带位于鄂尔多斯块体东缘，整体上呈“S”型展布，长约1 200 km，宽约几十公里，属于拉张型断陷。山西断陷带南部主要包括临汾盆地和运城盆地两个断陷盆地。其中临汾盆地由洪洞凹陷和临汾凹陷组成，受控断裂主要有苏堡断裂、罗云山断裂以及大阳断裂。1303年和1695年先后发生在洪洞凹陷和临汾凹陷的2个历史大地震，使临汾盆地划为高烈度区。运城盆地的主要活动断裂呈NEE向，西南长边走向NEE至近EW向，往东北逐渐变为NE、NNE向,盆地由两凹一垒组成，即侯马凹陷、运城凹陷和峨嵋台地。侯马凹陷受吕梁山山前龙门山断裂控制，单侧断陷，新生界厚约1 200 m；运城凹陷受中条山北麓断裂控制，单侧断陷，新生界达5 500 m[13-16]。研究区内主要活动断裂空间分布及断裂几何特征如图1和表1所示。

图1 研究区地质构造简图Fig.1 Geological structure of the study area

断裂编号断层名称走向倾向倾角(°)运动性质长度kmF1苏堡断裂NWWS高角度正断层30F2罗云山山前断裂NNESE80°以上右旋兼正断145F3韩城断裂NE30°NEE60°左右正断层100F4峨眉台地北缘断裂NE70°N70°～80°正断层156F5峨眉台地南缘断裂NESE70°～80°正断层50F6中条山断裂NE-NEENW58°～75°正断层137

1.2 模型的构建

基于研究区断裂分布以及GPS站点位置，在经度109°E～113°E、纬度34°N～37°N区域内建立山西断陷南部三维有限元模型，垂向上包括上地壳、中地壳和下地壳，建模时采用高斯投影变换将研究区内的大地坐标转化为平面坐标。研究区断裂带由很多细小的断层系、断层束构成，空间分布繁杂，因此在构建过程中做一定的简化，即按照断裂带的主体走向，将属于同一断裂带的断层系处理为一条断层，将断层作为相对容易变形的软弱化带处理，同时适当增加断层及其附近区域的网格密度，确保之后的模拟计算具有比较好的精度和稳定度。根据研究区地壳速度结构的反演结果，综合设定上地壳为15 km、中地壳为12 km、下地壳为12 km，断层带宽度约5 km左右，断层深度延伸至中地壳底部。采用有限元静力学分析方法，假设研究区地壳为弹性结构，具体的材料参数如表2所示，利用ANSYS软件单元库中solid185划分单元，最终建立模型的节点数为375 143个，单元数为2 379 224个(见图2)。

表2 研究区弹性材料参数表Table 2 Elastic material parameter of the study area

图2 研究区三维有限元模型Fig.2 Three-dimensional finite element model of the study area

1.3 模型边界约束条件

一般认为，在几年时间尺度范围内，岩石的力学性质主要以弹性为主，区域应力的大小和方向也认为是基本稳定，因此可以将某一时段的GPS速度场作为模型边界约束条件，得到该时期研究区内应变场增量的空间分布图。目前由于GPS点位分布较稀疏，同时也存在一些速度不连续的点位，模型边界约束条件需要通过节点周围的GPS点进行插值计算，因此需要将误差较大、速度不连续的点位进行剔除，插值计算过程中尽量避免断层附近及断陷带内部的点，选择较为稳定的GPS站点。假定速度场在地壳的垂直方向上不随深度的变化而改变，将经插值计算得到的边界约束条件加载到模型周缘，同时将模型底部的垂直方向固定，水平方向设为自由面，模型表面不做任何约束，计算时将模型视为无自重弹性体。本次研究采用中国地震局第一监测中心提供的中国大陆GPS水平形变速度场，参考基准为区域无旋转基准。原始观测数据来源为“中国地壳运动观测网络”和“中国大陆构造环境监测网络”项目产出的连续和流动观测数据以及IGS提供的连续观测数据，第9页图3为研究区2011-2016年GPS水平形变速度场。由图3可知，2011-2016年研究区内变形速率较小，整体运动方向SEE向。该速度场所反应的变形性质与地质研究结果所给出的山西地震带NW-SE向拉张，NE-SW向压缩的应力场状态较为一致。

图3 研究区内水平形变速度场Fig.3 Horizontal deformation velocity field in the study area

2 计算结果分析

本文所述的应变场为该计算时段的应变场增量，不代表各断裂实际的应变积累程度，定义应变场结果的正值表示拉张变形，负值表示挤压变形，剪切应变为水平面方向的剪切应变。

通过GPS速度场给定位移约束，得到2011-2016年时段的主压形变场和剪切形变场图像(见图4)，结果显示，该时段整体形变量较弱，除断裂端部、断裂转折处应变量约为10-7量级，断层整体的变化量级约10-8，表明汶川地震对山西断陷带的影响逐渐减弱，但区域形变场仍处于弱挤压形变状态。在中条山断裂附近和罗云山山前断裂附近有相对较高的应变增量集中区，同时剪切形变也具有一致性，罗云山山前断裂附近的临汾跨断层水准也显示2013年之后断层张性活动明显减弱，逐渐转为压性变化(见图5)，与该时段数值模拟计算结果较为一致。2016年盐湖ML4.8地震发生在中条山断裂应变增量集中区附近，表明应力场的扰动对地震有促进作用，因此通过分析应力场的相对变化来研究形变场分布与地震的发生具有一定的意义。

图4 研究区应变场分布图Fig.4 Distribution of strain field in the study area

图5 临汾跨断层4-1水准测线观测值Fig.5 Observation of the 4-1 leveling line of Linfen cross-fault

3 结论

在收集研究区断层几何参数以及地质材料参数的基础上，运用ANSYS软件构建三维地质模型，将GPS观测结果作为模型边界约束条件，计算得到研究区的主压应变场以及水平面方向的剪切应变场空间分布图像，同时结合2016年盐湖ML4.8地震进行了验证分析，初步得到以下结论：

(1) 数值模拟结果显示，2011-2016年山西断陷南部区域以挤压形变为主，跨断层水准测量结果也在该时段出现张性活动减弱状态，分析可能受区域应力场调整的影响。

(2) 2016年盐湖ML4.8地震发生在形变场累积集中区附近，表明利用形变场的相对变化来研究、预测形变场的扰动与地震之间的关系有一定的意义。

(3) 断层的几何形态会对应变场空间分布有一定的影响，尤其是在断层端部、转折处以及断裂交汇部位。

(4) 模型计算的应变场表示的是区域应变场的应变增量，断裂危险性的大小主要取决于断裂绝对应变能的积累程度，因此将某一期的GPS速度场作为模型的约束条件仅能对本时段各断裂应变积累的相对量作一定的比较。

衷心感谢中国地震局第一监测中心和中国地震局地震研究所为本次研究提供GPS观测数据。