地震波诱发的深部地层动应力机制及断裂位移研究−以汶川地震对龙门山前陆盆地的作用为例*

张子平 袁 青 熊齐欢

1)中交第二航务工程局有限公司，武汉 430040

2)长大桥梁建设施工技术交通行业重点实验室，武汉 430040

引言

2008年5月12日四川省汶川县发生了8.0级大地震。地震造成了巨大的经济损失和民众伤亡。伴随着地震的发生，周边地区地层产生了众多的形变构造，给后续的工程重建和人民财产安全带来了许多隐患。

通过对震后遥感记录和实地考察所得资料的调查研究，发现汶川地震后诱发了众多的滑坡和褶皱断层。曹俊兴等[1-2]发现在距离汶川地震破裂带30—40 km的龙门山前陆盆地地腹构造处曾存在断裂和褶皱变形。这些研究表明在汶川地震的影响下，近场至中场区域的地腹岩层出现了断裂与褶皱变形。断层作为一种分布极为广泛的构造形态，在许多构造环境中均有发育。它不仅是一个与构造地质学、地震学、地球动力学有关的重要地质现象，也同矿产资源的形成与分布、大型工程的基础稳定性、地震危险区划分等有十分密切的联系。

近年来，关于地震与断层间的应力应变研究已有发展，早在20世纪中后期，地震学家们就用库仑破裂应力来描述地震引起的地层应力应变效果[3-6]。1951年，Anderson[7]根据断层形成时3个主应力的作用方向与相对大小，提出了断层作用理论。Papadimitriou等[8]利用中国四川西部地震活动区域内地震数据研究了1893年以来库仑破裂函数的变化以及近110年该地区应力场的演变，计算出地震引发的弹性半空间模型中的静态位错。戴俊[9]基于爆炸应力波随距离衰减原理研究冲击波如何转变成地震波。Aiken等[10]引入流体的影响来研究大范围地震的动态库仑应力作用。Okada[11]利用半无限空间中任意倾角的倾向滑动剪切断层错动位移量计算出产生的应力场与位移场。Wang等[12]利用格林函数来计算地震区域内的水平分层粘弹性介质模型。前人对地震作用下岩石构造的破坏分析历经多年的研究，常采用的研究方法有Newmark分析法、动力有限元时程分析法[13]、动力有限元强度折减法[14]等，这些方法或理论分析，或构建模型，均综合考虑了地震影响下岩土体的破坏性问题。

这些研究方法虽然成果丰富，但少有对于深部地层破坏的研究，且没有过多考虑岩层间波的折射衰减效应。基于此，本文通过波动力学和爆破力学衍生出的力学方程和地震波的层间折射效应，提出岩层在地震波作用下的受力机制。依托汶川地震对龙门山前陆盆地深部地层的具体破坏效应，根据实际地层的物理力学参数，对深部地层进行应力平衡分析。建立地层与断层面模型，设置应力约束和阻尼边界，据Coulomb准则和摩尔应力圆判断是否发生剪切破坏与位移情况。根据最后的位移量结果与监测主应力情况，得出在汶川地震波的冲击作用下，龙门山前陆盆地深部地层的断层破坏形变情况。同时证实此地震波应力计算机制的可行性以及大地震对断层构造的诱发影响。

1 计算原理

1.1 一维情况下地震波应力效应

天然地震传播的地震波主要包含纵波和横波，而由此衍生的应力也是由横波和纵波产生的水平应力和垂直应力。其应力效应在一维情况下可由波动力学公式求解[15]：

式中，σ为地震波产生的动应力，纵波影响垂直应力，横波影响水平应力；ρ为岩体密度；c为岩体中的波速；µ为岩体中质点的振动速度。爆破地震学[16]中由动量守恒、质量守恒和能量守恒三大定律导出的岩石中冲击波参数计算公式：

式中，g为重力加速度；U为波阵面上质点的运动速度；Ub为冲击波传播速度（地震波波速）；γ0为比重（震前岩石密度）。可以看出，两者的计算公式具有高度的一致性，因为天然地震波是弹性波[17]，而爆炸发出的冲击波和应力波向外传播后也会衰减为弹性波。

1.2 地震波在分界面的折射与应力衰减

岩层密度和岩体中的波速均可由实验或当地地质情况得到，故岩体中质点的振动速度就是所需要测算的关键因子。对均质地层分层后进行横波纵波透射后的位移量计算：

假设层面一共有n+1个层面（图1），下部介质为第1层，上部介质为第n+1层，中间n−1层。P波和S波的入射波为Pr、Sr；反射波为Pf、Sf；透射波为Pt、St。其中的入射角为。其中上部介质，第n+1层的介质的总位移量（包括入射波与反射波）为[18]：

图1 地震波在厚地层中的反射与透射Fig. 1 Reflection and transmission of seismic waves in thick formations

P波：

S波：

式中，A，B是P波和S波的振幅；j为常数，ω是波形数据。且，这里的αn+1和βn+1对应为相应层面的波速，与其对应的入射角有关，，C即为射线参数。将式（3），（4）对时间t求导后分别代入式（2），在考虑透射系数的情况下可以得到P波、S波在通过地层分界层面（两层岩石组成不同、物理性质不同）后对目的地层的应力影响：

式中，j为常数。ρ1为所要求解层面的密度，cP1为其P波波速，cS1为其S波波速；ρ3则为入射进此层面的初始层面的密度，cP3为其P波波速，cS3为其S波波速。iP，iS为P波和S波的入射角。x是目的层位距离震源的水平距离，z是目的层位距离震源的竖直距离。A、B分别为P波、S波对应的振幅。ω为波形函数。由于地震波速多为不规律的，因而采取某一时刻的数据进行计算，t为对应的时间。

2 模型构建与应力计算监测

2.1 模型构建

龙门山前陆盆地位于龙门山冲断带与龙泉山前陆隆起之间，北起安县秀水，南抵名山、彭山一线，西部已卷入龙门山冲断带[19]。龙门山前陆盆地自古以来就充填了大量的海陆相沉积物[20]，以时间可分为三叠系、侏罗系、白垩系、古近系、新近系和第四系6个层位。根据本文的研究内容，基于前人勘探资料，对岩层的物理力学参数进行平均化计算。共计建立8个层位，表1为其物理力学参数[21-23]。

表1 各地层的物理力学参数与围岩应力Table 1 Physical mechanical parameters and surrounding rock stress of each layer

基于研究地震波引发动力对深部地层造成的破坏效应，将地层看作一个整体构造，建立15 km×5 km×10 km的立方体模型（图2），内部设置斜面模型观测地震波动力效果下断层面上的应力应变情况。由于采集数据的地震台站所处区域的断层带多为30°倾角的逆断层，故设计此结构对其进行分析，通过多次研究发现，20°—45°倾角范围内整个地层结构体的受力和位移情况相差不大，其他倾角情况下，倾角越小，越容易产生错断位移情况。

图2 主体地层模型Fig. 2 The model of main stratigraphic

基于莫尔—库仑失效准则对于断层失稳的判断分析，对断层前后层面周边的地层情况进行数值模拟，通过监测岩层的应力与位移形变信息来具体分析。设置在斜面上和周围区域共13个监测点（图3）。

图3 监测点在XOZ平面上的投影Fig. 3 Projection of the monitoring point on the XOZ plane

2.2 应力计算与监测

震源距离所取数据的地震台站水平距离36 km，如图4所示分为9等份，每4 km取一个质点。z对应每层的厚度（图4）。采集地震台站（江油市含增镇编号51JYH）记录的加速度时程信息作为初始条件，模拟地震波影响下的动力阻尼。

图4 地下分层情况与质点分布Fig. 4 Underground stratification and particle distribution

对岩石地基，在底部施加加速度转化的应力荷载。采取以0.02 s为间隔的12500个样本时间数据（图5）。根据式（5）、（6）以及表1中的岩石物理参数计算对应的各地层质点P波、S波在通过地层分界层面（两层岩石组成不同、物理性质不同）后对目的地层的应力影响（以第7层质点2为例）（图6）。

图5 加速度时程信息Fig. 5 Acceleration time history information

图6 第7层质点2处的P波和S波应力影响Fig. 6 The influence of P-wave and S-wave stress at the particle #2 of the 7th layer

可以看出，P波和S波的应力在（−100，100） MPa区间内。S波在初始的10 s内峰值达到−36 MPa左右，而P波在初始的10 s内峰值达到−50 MPa。其中正数值代表拉应力，负数值代表压应力。前50 s内的应力数据具有代表性，应力伴随时间完成由0瞬间增加至压应力极值再逐渐递减的转变，也与地震导致岩石破裂的时间范围接近，后200 s的数据由于应力积累等原因呈现周期性反复。而由于地震波中横波应力效应一般大于纵波，也在图中体现为P波的峰值大于S波。

通过将式（5）、（6）计算好的各地层地震波P波和S波动应力影响施加于建立好的模型体上，发现随着地震波的震荡变化，地震应力和不平衡力变化（图7）也呈现上下波动的趋势。

图7 动力过程中的最大不平衡力Fig. 7 Maximum unbalanced force during dynamic process

随着地震波加速度数据的正负转化，与其成正相关关系的地震应力也随之发生变化，导致最大不平衡力一直在0与最大值之间波动，且由于大尺度模型的自由边界问题，未出现平衡收敛。从这也可以推测出，地层结构已出现了破坏。对静力平衡下及动力施加过程中上盘3个点的主应力情况（图8）进行统计。

图8 动力过程中的主应力变化Fig. 8 The change of principal stress during dynamic process

根据表2的数据绘制摩尔应力圆和大小主应力的变化曲线，蓝色摩尔圆由静力下的最大最小主应力峰值所绘制，选取此情况为极限平衡条件的摩尔圆。从动力施加下的主应力情况来看，110000动态时间步后主应力趋于稳定（表2）。各点上最大主应力的变化情况均大于最小主应力的变化，判断红色主应力圆会相交于同监测点蓝色稳定情况的摩尔圆（图9），即此位置的层面已发现了断裂破坏。初步判断在汶川地震影响下龙门山前陆盆地深部地层出现了断裂构造。

图9 监测点的主应力变化Fig. 9 The change of principal stress at monitoring points

表2 监测点的主应力变化Table 2 The change of principal stress at monitoring points

3 监测点与模型体的位移分析

3.1 监测点位移变化

主应力可初步判断是否发生断层错位，监测各点各层位的位移变化可以验证每层的错断与否并统计其位移变化情况。图10为13个监测点的Z向位移和中间斜面监测点的X向位移记录。

图10 动力过程中的监测点位移Fig. 10 Displacement of monitoring point during dynamic process

在整个动力过程中每个监测点的Z向位移变化趋势基本相同，初始都为0，随时间不断呈反向递增趋势，即整个模型体在纵向上不断发生负向位移。整体上，上部地层的Z向位移比下部地层要少，因为地震的加速度是直接施加在模型体的底面上。而监测点9至监测点10、监测点11至监测点12两个层面间的Z向位移量相差无几，此处的纵向位移不太明显。

而对于中间斜面监测点的X向位移来说，点7：初始位移为0，108 000 step（step为内置时间步数，代表程序运行时间，基本400 000 step对应250 s）之后向X正向递增到7 m，再负向位移，158 000 step回到初始位置，继续负向位移至−28 m左右。点8：初始位移为0，104 000 step之后向X正向递增到7 m，再负向位移，140 000 step回到初始位置，继续负向位移至−47 m左右。点9：初始位移为0，100 000 step之后向X正向递增到6 m，再负向位移，124 000 step回到初始位置，继续负向位移至−70 m左右。点10：初始位移为0，100 000 step之后向X正向递增到4 m，再负向位移，117 000 step回到初始位置，继续负向位移至−84 m左右。点11：初始位移为0，100 000 step之后向X正向递增到3 m，再负向位移，105 000 step回到初始位置，继续负向位移至−118 m左右。点12：初始位移为0，基本一直负向递增位移，最后达到−136 m左右。点13：初始位移为0，基本一直负向递增位移，最后达到−168 m左右。

总体上说，各监测点的X向位移均有向负向递增的趋势。对于上盘3个点来说，深度越大，变化幅度越小，峰值越小。而对于下盘的3个点来说，深度越大，峰值越小。在斜面上的监测点呈现出与其Z向位移相当的特点，从整体上来说，深度越大的层面，其发生的X向位移量越小。且在110000 step这个应力激发时刻7、8、9、10，4个监测点出现了正向位移。之后4个点又向负向位移，最后稳定为负向位移。综合3个图，发现每个监测点的X向位移均有所不同，即每一层位都存在X向的位移偏差。

3.2 模型体位移变化

对每100 step的时间顺序依次截取X-Z面的剪切应变变化进行分析（图11）。

图11 X-Z面监测的剪切应变变化Fig. 11 Shear strain change monitored by X-Z plane

可以看出，因为施加的地震数据是置于底面的，整个模型的错位是从底层开始的，红色区域的层位首先发生断裂，然后往上，黄色，绿色，蓝色层位随着时间的推移相继发生断裂，最后形成断层构造。在整个过程中，下部岩层的剪切形变量均大于上部岩层。2 500 step后全部层位已基本完成错断，与1、2、3监测点的主应力信息相符，表明这些位置已发生断层偏移现象。地震加速度释放的350000 step时刻后，全部层位基本发生错位。而与剪切应变随时间变化的特征相类似，作为能代表整体形变特征的体积应变变化也是从底面开始，随时间推移，上部岩层慢慢完成层位的错断。伴随着分界斜面，也就是各图中的斜线，上面的岩层相对斜线以下岩层位移向上偏移。将斜线看作假想断层面的话，断层上盘向上位移，下盘向下位移，即发生了逆断层。而整个模型体的X向位移本来应该是连续而平缓的，由于上盘岩层向上位移，下盘岩层向下，导致整体X向位移出现不连续的断裂。

4 讨论与结论

本文提出了一种考虑地震波分层折射影响的地震应力计算机制，并对龙门山前陆地腹构造受汶川地震影响的应力与形变进行数值分析，探究大地震对深部岩层断层构造的诱发机制。得到以下结论：

（1）利用分层介质传播衰减规律计算出的汶川地震波应力影响能对应各时间点地震波计算出同时期的应力影响，并考虑了波在地层分界面之间的折射（地震波能量衰减或增强）效应，使结果具有时间上的连续性与动态效应，模拟出岩层破坏的变化趋势，反映地震波对深部地层的动应力影响及位移—时间变化。

（2）通过记录上盘监测点的最大、最小主应力变化情况，初步判断龙门山前陆盆地地腹岩石出现了剪切破坏情况，得出汶川地震确实引发了深部地层的断裂构造。跟踪这些点的X与Z向位移，发现各层位位移量出现差距，也佐证了这一点。在本次的模拟中，各层的Z向位移差别较小，而X向位移偏差较大，且随着深度的变化具有较为明显的差异性，深度越大，位移量越小。

（3）通过对模型X—Z面的剪切形变与体积形变进行分析，可以很明显的发现各层位都有错断结构，且根据体积应变的变化趋势和最后结果，可以判断断层为逆断层。对整个模型体在地震影响下的X向位移结果进行分析也可佐证此点。