跨活断层隧道抗位错精细化分析方法研究

2022-09-15赵伯明赵天次周玉书

铁道学报 2022年8期

赵伯明，赵天次，周玉书

(1.北京交通大学城市地下工程教育部重点实验室, 北京 100044；2.北京交通大学土木建筑工程学院, 北京 100044)

活断层在我国分布广泛，突发性的断层错动是产生地震的主要根源，长期以来给人类的生命财产和结构工程的安全性带来了很大威胁。大量研究表明，活断层错动导致的同震地表破裂、土体变形、强地面震动及其他次生地质灾害均会给穿越活断层隧道带来危害[1-11]。然而在大型工程的选线和设计中不可避免的会穿越或临近活断层，因此跨活断层隧道的抗位错研究具有重要的理论意义和工程价值。

目前针对穿越活断层隧道抗位错分析国内外开展了诸多研究，研究方法包括现场调研、模型试验、解析计算和数值模拟等。

现场调研可以直接得到最真实的一手资料，通过震后短时间内的仪器监测可得到断层附近的地表同震位移[12]，运用回归分析手段将实地调研数据进行处理，建立地表同震位移与震级、上覆土层厚度等指标的定量关系[13-17]。但该方法受限于观测点的数量和密度，不同区域的地质构造和土层性质差异也会导致回归关系式的差异[18]，且大多数研究重点关注了地表的同震位移，对隧道等地下工程而言，土体内部的变形对结构受力影响更大，断层错动导致的围岩状态变化同样值得关注。

部分学者针对断层错动对隧道结构的影响进行了模型实验研究[19-24]，但在实验中由于模型尺寸的限制，对断层错动的模拟通常简化为单侧底板的抬升，并不能模拟真实的土体变形情况。

在基于弹性半空间的解析解法方面，学者们[25-27]应用简化的断层模型和位错理论，给出了半无限空间中矩形滑动断层的地震位移场的解析解。该方法推导清晰，操作便捷，其中应用最为广泛的是Okada提出的均匀半无限空间解析解法，但其缺点是对于多层介质模型或三维复杂情况并不适用。

随着地质资料和勘测手段的进步，Okada的解析解法已经略显粗糙，在均匀半无限空间解析解基础上进一步发展而来的层状半空间中的静态弹性变形计算方法考虑了地球的层状构造，位移矢量需要在柱坐标系下利用柱谐函数来表达，并通过传播矩阵的算法求解[28-29]，近年来已经逐渐成为主流趋势。

对于地球层状构造对同震变形的影响程度，文献[30-31]的研究表明：在地壳内发生的地震，在震源距不大于100 km的情况下，层状构造的影响可达到20%。

因此，随着计算理论的发展，人们对断层位错发生后土体的变形了解逐步加深，有限元模型中的输入也应尽量接近真实的构造环境，以提高计算结果的可靠性，更加精细化的位错输入方法有待提出。

本文首先对现有的跨断层隧道抗位错计算方法进行误差分析，宏观上分析了土体层状构造和断层破裂面上位错分布导致的同震变形场差异，微观上分析了不同输入方法导致的土体变形差异，揭示了现有隧道抗位错计算中同震变形输入方法的不足。然后基于正交化格林函数方法编制程序，对层状半无限空间中的静态弹性变形进行计算，结合Abaqus的input文件将断层错动导致的土体变形场以位移边界的形式施加到有限元模型，并采用Matlab建立了输入流程的GUI界面，使建模流程大为简化。最后将该方法应用于成兰铁路某隧道的三维抗位错分析中，实现了隧道工程抗位错分析中断层错动的精细化输入。

1 现有跨断层隧道抗位错机理研究

我国西部地区地质条件复杂，活断层分布密集，随着我国西部地区铁路隧道逐渐增多，部分隧道不可避免的需要穿越众多的活断层。当隧道工程处于断层附近时往往需要进行抗位错分析，以得到隧道结构在断层错动作用下变形效应和破坏机理，隧道工程与活断层之间的的空间关系见图1，黑色单箭头代表断层的错动。

图1 跨断层隧道空间位置关系

在隧道工程的抗位错分析中，通常截取有限区域为计算区域，并设置人工边界来模拟活断层的位错[10,32]，但现有的建模过程中往往存在着过度的简化，简化分析模型见图2。

图2 现有研究中的简化模型[32]

由图2可知，上述简化建模方法中将断层破裂面直接置于隧道中部，通过断层上盘底面的整体抬升模拟逆断层的错动。而实际地层环境中大部分断层的破裂面并未直接截断隧道，对于埋深不为零的隐伏断层或破裂面未直达地表的断层而言，该模拟方法显然是不适用的。此外，上述研究中将断层破裂面上的位错量设置为均匀分布，与实际断层面上的位错分布存在较大的差别。通过分析可知建模过程中的误差来源如下：①土体的层状结构；②破裂面上的位错分布；③断层错动的输入方法。本节以典型地震为例，通过对比分析得到各误差源的影响大小，说明了建立精细化位错输入方法的必要性。

1.1 土体层状结构对地表同震变形的影响

土体的层状构造在同震变形场的计算中存在影响，为量化这一误差，本节选取典型地震2013年Mw6.7级芦山地震作为算例，对成层土体模型和均匀土体模型下的地表同震变形进行对比分析，研究土体层状结构对同震变形场的影响。

由王卫民等[33]给出的芦山地震反演结果可知，发震断层为一条埋深2 km的隐伏逆断层，断层破裂面并未延伸到地表，断层参数见表1。断层参数及坐标系如图3所示。其中断层长度为L，宽度为W，倾角为β，断层面上位错滑动角为α，滑动量为u。

图3 断层空间坐标系

表1 芦山地震震源参数

成层土体模型的土层划分参考参考芦山地震震源区地壳地震波速度模型[34]及前人的研究成果[35-41]，综合给出了芦山震区的成层土体模型，自地表向下共划分了4个不均匀波速层，土层参数见表2。均匀土体模型与成层土体模型的最上层土采用相同的力学参数。

表2 芦山地震震区地壳结构模型

首先计算均匀土体模型中芦山地震导致的地表同震变形，竖向位移等值线和水平位移矢量分别见图4，其中虚线矩形框表示断层的地表投影。

图4 均匀位错模型地表竖向、水平位移矢量图

地表位移的分布受断层产状影响，竖向位移的数值总体上大于水平位移。由图4(a)可知，地表隆起主要分布在断层上盘区域即Y轴正半轴，沉降集中下断层下盘区域，同震变形在断层上断线的地表投影附近变化最为剧烈。竖向位移的隆起峰值为0.802 m，沉降峰值为0.062 m，这是由于断层倾角较小，断层面上的位错量经过土层传导出现了一定的偏移。

由图4(b)可知，上盘和下盘的土体均向断层上断线的地表投影挤压，随着与上断线地表投影线距离的减小，水平位移值逐渐增大，峰值为0.338 m。

基于分层土体模型的芦山地震地表同震变形计算结果见图5，从结果来看，地表竖向位移和水平位移的分布和数值都与GPS监测数据[42-43]给出的结果相近。

图5 分层位错模型地表竖向、水平位移矢量图

对比图5两种模型的地表同震变形计算结果，考虑地壳层状构造之后地表的竖向位移和水平位移在分布规律上绝大部分都保持一致，竖向位移隆起峰值为0.775 m，沉降峰值为0.046 m，水平位移峰值为0.395 m，但是在峰值位置上有所差异。均匀土体模型的隆起峰值和沉降峰值都要大于分层土体模型，但水平位移峰值小于分层土体模型，这一差异是由均匀土体模型较分层土体模型更软的土体力学参数导致的。为清晰的展示两者差异，将两种模型的计算结果取差值，见图6。

图6 分层模型和均匀模型的地表竖向、水平位移差值

由图6(a)可知，竖向位移的差异主要集中在上断线投影附近，地表隆起存在一条沿断层走向的0.1 m以上分布带，地表沉降则对称分布在上断线投影，量值上仅0.05 m左右。这说明均匀模型具有更大的隆起和沉降峰值，竖向位移最大差值为0.308 m，约占分层土体模型竖向位移最大值的39.7%。

由图6(b)可知，水平位移的差异同样体现在上断线投影处，差异带宽度10 km左右，最大差值0.236 m，达到分层土体模型水平位移最大值的59.8%。

由芦山地震震例的计算结果可见，土体的分层结构对地表同震位移的影响是不可忽略的，位移差异达40%～60%，尤其是在靠近断层上断线地表投影区域，需要考虑地震发生地区地壳层状构造的影响。

1.2 破裂面上位错量分布对地表同震变形的影响

在真实的地震中，断层破裂面的几何形状与位错量分布都很复杂。在以往的跨断层隧道计算中，往往假设断层面位错是均匀分布的，但实际断层面上的位错分布是非一致的。吴忠良[44]通过对多个地震断层面上的滑动量调查统计发现断层面上的位错分布是高度不均匀的；Banerjee等[45]利用41个连续观测的远场GPS资料对Sumatra-Andaman地震断层的滑动量进行计算，发现断层上的最大滑动量出现在断层中部附近。

在经典半无限空间的位错理论中，Okada均匀土体模型的有限矩形位错源公式是由点源公式在断层面上积分所得，分层土体模型中有限断层源是由离散的点源格林函数线性叠加而得，两种模型均可对断层进行离散化，使破裂面上的位错量接近真实的位错分布，再求得对应的同震变形场。

为分析断层破裂面上位错量分布对同震变形场的影响，选取典型地震1989年美国Mw6.95级Loma地震分别计算位错量均匀分布和真实分布两种情况的地表同震变形，Loma地震发震断层参数见表3，Loma地区的土体参数见表4[47]。

表3 Loma地震震源参数

表4 Loma地震震区地壳结构模型

沿走向和倾向将发震断层划分为2 km×2 km的子断层，每个子断层上的位错量反演结果见图7[46]。

图7 Loma地震断层破裂面位错量分布(单位：cm)

图7中每个子断层上位错量的精确反演结果反映了断层面上位错的真实分布情况，再求取各个子断层面上位错量的平均值作为均匀分布模型的位错量，采用分层土体模型分别计算Loma地震两种位错分布模式下的地表同震位移。

真实位错分布模式下地表三个方向的位移等值线分别见图8。

图8 真实位错分布下竖向、X向、Y向地表位移等值线(单位：m)

由图8可知，从结果来看地表位移分布符合斜滑断层的运动机制。在断层倾向抬升作用和走向挤压作用的叠加影响下，竖向位移有整体向左移动的趋势，竖向位移峰值在断层中部左侧约7 km的位置，峰值量约为0.7 m。X方向位移主要受走滑分量的控制，表现为上盘区域沿走向向左移动，下盘区域沿走向向右移动，X向位移峰值在断层中部左侧约12 km的位置，峰值量约为0.4 m。Y方向位移受到走滑分量和倾滑分量的共同作用，表现为大部分区域向Y正方向移动，只有断层右侧小部分区域向Y负方向移动，Y向位移峰值在断层中部左侧约8 km的位置，峰值量约为0.3 m。

为了量化位错分布模式对地表同震位移的影响，将均匀分布模式下的结果与真实分布模式的结果取差值进行分析，三个方向上的位移差值分布分别见图9。由图9可知，均匀分布下的三方向位移在整个断层区域内都与真实分布有较大的差别，误差最大的位置对应位移峰值所在位置。

图9 位错均匀分布与真实分布的竖向、X向、Y向地表位移差(单位：m)

考察断层走向中线剖面位置处的地表同震位移见图10。由图10可知，三个方向的位移值都随着远离断层而逐渐减小，随着与断层距离的减小，均匀分布模式会产生较大的误差，具体数值见表5。

图10 断层走向中线剖面的地表同震位移对比

表5 不同位错分布模式的位移峰值误差

由表5可见，均匀分布模式在水平和竖向都会给地表同震位移带来误差，水平同震位移峰值误差在X方向达到36.6%，在Y方向达到8.4%；竖向同震位移峰值误差在Z方向达到9.6%，这些误差放大了断层错动的破坏效应。

综上所述，将断层面上的位错量按均匀分布考虑会给计算结果带来误差，在计算土体的同震位移场时，采用反演得到的真实位错分布具有一定的必要性。

1.3 断层错动输入方法对地表同震变形的影响

在求得断层错动导致的土体同震变形场之后，将该同震变形场准确合理的输入截取的计算区域是保证跨断层隧道抗位错分析准确性的重要一环。现有研究中通常采用断层上盘底面的整体抬升模拟逆断层的错动，见图2，这是一种过度简化的手段，显然会带来一定的误差。本小节选取矩震级Mw7.0级典型地震对断层错动导致的模型底面同震变形进行计算，并将模型底面的真实变形与简化方法的输入值进行对比，量化了简化方法的误差，说明了选择合理断层错动输入方法的必要性。

为使对比结果更为清晰，忽略土体分层和位错分布的影响，按均匀土体模型且位错均匀分布计算，土层参数设置为表2中的第一层土，选取逆断层进行分析，由断层破裂参数和矩震级Mw之间的统计公式[46,48-49]得Mw7.0级典型地震的震源参数见表6。

表6 Mw7.0级典型地震震源参数

选取断层上断线中轴面为观察剖面，断层模型及观测剖面位置如图11，分别计算深度为0、50、100、150、200 m处观察剖面的土体水平和竖向同震变形值，将各深度真实的土体同震变形值和简化输入变形值对比，见图12。

图11 断层模型及观测剖面

由图12(a)可知，竖直方向土体的真实同震变形与简化方法中输入的土体变形值趋势一致，随着深度的增加简化值与真实值的差值逐渐减小，上盘区域差值约0.2 m.下盘区域差值约0.25 m，误差为真实变形值的25%左右；由图12(b)可知，而水平方向土体的真实同震变形与简化输入变形值差异显著，简化方法中认为断层上盘固定不动仅下盘抬升的输入方式显然是不准确的，上盘和下盘的输入差值都为0.4 m左右，水平变形输入误差达80%左右。

2 跨断层隧道抗位错精细化计算

由上节分析可知，土体层状结构、破裂面上的位错量分布及错动输入方法均会导致地表同震变形的计算出现误差，其中影响最为显著的是错动输入方法，误差达25%～80%；土体层状结构影响次之，误差为40%～60%；破裂面上的位错量分布影响最小，仅10%左右，该误差在位错量反演不明确的条件下可忽略。

因此，本节提出了一套精细化的位错输入方法，首先由断层位错出发，求得断层错动导致的模型边界上的同震变形场，再将该同震变形场以位移人工边界的形式施加至有限元模型的边界上，以实现断层错动的精确输入。以成兰铁路跨越龙门山破裂带某隧道为例，探究隧道在断层错动作用下的破坏机理及力学特性，为隧道工程抗位错措施的建立提供指导。

2.1 精细化位错输入方法的实现过程

早期研究中将断层上盘或下盘的移动简化为模型单侧底面的抬升，这是较不合理的一种简化输入方法。本文提出的精细化断层错动输入方法建模思路见图13。

图13 精细化输入方法建模思路

计算流程如下：首先输入断层参数和土层参数，求得断层错动导致的有限元模型区域的同震变形场，再对计算结果进行提取，获得模型边界处的变形值并写入input文件，最后通过input文件将同震变形场施加到Abaqus有限元模型，以实现断层错动的精确输入，自编程序计算流程见图14。

图14 计算流程图

2.2 跨活断层隧道数值计算模型

以我国西部成兰铁路穿越北川—映秀活断层某隧道为例进行抗位错分析，北川—映秀是龙门山中央断裂带的主断裂，也是汶川地震的发震构造之一，本文将其简化为逆冲断层，断层参数见表6。

由图12可知，断层位错在土体产生的位错分布主要集中在上断线投影前后500 m范围内，在这个范围外的土体位移值虽然很大，但是位移分布一致，位移差值很小，隧道结构的破坏主要是由土体位移差值造成的，所以有限元模型的规模取纵向长为500 m。考虑围岩影响范围为3～5倍洞径，模型宽和高分别取为100 m、100 m，Abaqus模型见图15，坐标系原点为图中模型左下角点，模型整体位于第一象限。

图15 Abaqus有限元模型

为了验证数值模型的准确性，首先对不加隧道结构只有土体的数值模型进行计算，对该模型施加断层位错作用产生的土体同震变形场，提取地表沿纵向中线的水平和竖向位移值，与Okada解析解进行对比，结果见图16。由图16可知，从计算结果可以看出，本文精细化输入方法的数值解与使用Okada位错理论计算得到的解析解相差极小，验证了本文方法的准确性。

图16 有限元模型计算结果验证

下一步加入隧道结构形成土体-隧道有限元模型，隧道顶板覆土厚度为20 m，隧道断面选取圆形截面，尺寸参数见图17，模型中各材料的参数见表7。由于本文重点为为计算方法展示，为节约计算时间隧道衬砌结构采用弹性本构，若有更高的计算需求可替换为各类弹塑性本构，本文方法同样适用。土体采用Abaqus中的实体单元C3D20R模拟，隧道衬砌采用壳单元S4R模拟，土体-隧道有限元模型网格见图18，模型X轴为隧道横向宽度方向，Y轴为隧道纵向长度方向，Z轴为土体的竖向高度方向。

图17 隧道结构断面尺寸参数(单位：cm)

表7 有限元模型材料参数

图18 有限元模型网格划分

为计算断层错动而引发的隧道衬砌结构的受力与变形效应，将上述基于分层土体位错模型计算得到的土体同震变形场作为模型的位移边界施加在模型上来模拟断层的错动作用，得到隧道衬砌结构的破坏发展过程，断层错动前后的隧道变形绝对值云图见图19。

图19 断层错动前后的隧道变形绝对值

由图19可知，变形绝对值沿隧道纵向先减小后增大，隧道中部变形值的变化速率最大，为进一步分析隧道衬砌的受力分布情况，选取典型观测位置见图20，绘制隧道各处主应力随隧道纵向距离变化曲线见图21。

图20 隧道横截面观测位置分布

图21 各观测位置主应力随隧道纵向距离变化

由图21可知，隧道的拱腰和拱脚各处整体呈受压状态，最大压应力出现在拱腰附近，拱脚处也具有较大的压应力峰值；而拱顶和仰拱在断层上盘靠近上断线投影处呈受拉状态，这与断层上盘在该处的大幅抬升密切相关，最大拉应力出现在仰拱附近。由汶川地震的震害调研可知[24,50-51]，隧道的仰拱和拱顶主要出现抗拉破坏，拱腰和拱脚主要出现抗压破坏，拱肩损伤较小，这与计算结果相符。因此，本文基于分层土体模型提出的精细化断层错动输入方法可为工程中隧道的抗位错分析提供一定的指导作用。