0.   引言

地震灾害可能会对隧道结构造成严重破坏，严重威胁人们的生命安全和城市的正常运行。因此，为提高隧道在地震作用下的安全性和可靠性，隧道结构的灾变机制和抗震性能评价已成为岩土地震工程领域的重要研究方向^[1-2]。

目前，针对地下隧道结构数值模拟的研究通常选择有限数量的地震动作为输入^[3-4]。然而，地震动具有时间和频率非平稳特性的复杂随机特性，采用基于确定性分析方法的隧道抗震性能研究难以刻画不确定性条件下隧道结构的动力响应以及地震动的随机性。

针对目前常规方法的缺陷，概率密度演化方法能够通过概率密度函数（probability density function，PDF）充分考虑地震动随机性的特点^[5-6]，评估随机地震荷载下的结构抗震安全可靠性。在岩土工程抗震领域，概率密度演化方法的精确性和高效性已得到充分验证^[7-9]。然而，目前隧道地震响应研究主要集中在地震动的确定性分析，缺乏对地震动随机性特点的考虑。

本文依托典型工程场地中的隧道工程案例，将概率密度演化分析方法进一步拓展至隧道结构抗震性能评价。采用基于物理机制的随机地震动模型，生成符合抗震设计规范要求的一组非平稳地震动。建立考虑土-结构动力相互作用的隧道动力分析模型，通过非线性动力时程分析，结合概率密度演化方法得到隧道抗震性能评价指标（隧道倾斜角）的概率密度函数和累积分布函数。同时确定不同破坏状态下隧道倾斜角的超越概率，为隧道结构在不确定性条件下的动力响应分析提供依据。

1.   概率密度演化方法

概率密度演化方法结合概率守恒原则可以快速获得隧道结构地震响应概率密度演化特征。本文将简要介绍概率密度演化理论的概念及其数值求解过程，更多细节可以参考Li等^[5-6]的相关研究工作。

首先，将随机地震动模型生成的不同地震动加载到数值计算模型中，获得结构随机地震响应。在此基础上，采用偏微分方程代替传统的统计分析方法可以得到结构地震响应的平均值和标准差。随机地震作用下结构的动力响应控制方程可视为随机动力系统，关注的结构地震响应量Y可以用以下公式计算：

式中：H_y为地震响应函数，Θ为独立的随机变量，其概率密度函数为P_Θ(θ)。上述随机系统满足概率守恒条件，根据概率守恒原理的随机事件描述，可以获得非线性系统的广义概率密度演化方程，从而可求解结构地震响应量Y的概率密度函数：

其中，当t = ${t_0}$时，存在如下初始条件：

式中：${P_{Y\Theta }}(y, \theta , t)$为随机动力系统的联合概率密度函数；δ为狄拉克函数；${y_{\text{0}}}$为t=${t_{\text{0}}}$时Y的值。将上述所有的离散数值解进行累积求和可知，任意时刻Y的PDF可以用下式进行计算：

直接推导获取公式（4）的解析解相当困难，可采用Li等^[6]提出的方法对公式（4）求解。随机动力系统地震响应概率密度演化方程的求解流程如下。

（1）在概率空间中选取典型代表点并确定相应的赋得概率。在基本随机变量的概率空间中选取一组离散代表点如θ_q，同时确定了各个代表点对应的赋得概率，并针对初始条件进行了相应的初始化。针对每个选定的代表点，生成随机地震动。

（2）开展确定性动力系统分析。对土-结构系统开展大量确定性非线性动力时程分析，针对每个代表点$\Theta = {\theta _q}$得到所关心物理量对时间的偏导数${\dot Y_j} = ({\theta _q}, {t_n})$，其中j=1，2⋯，n，以获取不同随机地震动作用下结构的动力响应。

（3）概率密度演化方程求解。将每个代表点$\Theta = {\theta _q}$对应的结构响应量对时间的偏导数${\dot Y_j} = ({\theta _q}, {t_n})$引入广义概率密度演化公式（2），并根据初始条件和边界条件获得数值解${P_{y\Theta }}(y, \theta , t)$。

（4）对所有离散点求和得到最终解。将式（3）中获得的所有离散解相加求和或采用式（4）对${P_{y\Theta }}(y, \theta , t)$进行积分，可以最终获得广义概率密度演化方程的数值解${P_z}(y, t)$。

2.   数值分析模型

2.1   随机地震动模型

为合理考虑地震动的随机性，本文基于随机地震动模型和概率空间的剖分选取输入地震动，采用Li等^[10]提出的基于物理机制的随机地震动模型，该随机地震动模型已广泛应用于工程结构动力可靠度分析中^[11-12]。通过合理确定随机地震动模型中基本随机变量的均值和变异系数，进一步划分概率空间来选择代表点，以生成一组相应的地震动加速度时程曲线。本文采用切球选点法^[6]共生成282个代表点。每个代表点都被分配了相应的概率，最终得到了282条输入地震动时程曲线。输入地震动的峰值加速度的平均值和标准差分别为0.80 m/s²和0.57，且人工合成地震动的加速度反应谱与规范设计反应谱较为吻合。

2.2   场地计算模型与参数

本文选取的工程场地覆盖层厚度为70 m，该场地的土层情况及其物理参数如表 1所示，土层等效剪切波速为v_se=273.97 m/s，参考《建筑抗震设计规范：GB50011—2010》中工程场地类别的划分，为中国规范中的Ⅱ类工程场地，该场地为场地卓越周期T_G=0.4 s。表 1分别给出了剪切波速${v_{\text{s}}}$、不排水剪切强度S_u、密度和泊松比随地层深度的变化规律。

表  1  土层剖面参数表

Table  1.  Parameters of soil strata

土层	埋深/m	剪切波速${v_{\text{s}}}$/(m·s-1)	剪切强度Su/kPa	密度/(kg·m-3)	泊松比
#1	0~10	250	100	1800	0.30
#2	10~25	300	280	1800	0.30
#3	25~35	400	400	1900	0.30
#4	35~70	600	500	1900	0.30

下载: 导出CSV 

| 显示表格

4个地层的剪切模量和阻尼比随应变衰减规律由Ishibashi等^[13]提供的典型G-γ-D曲线反映，如图 1所示。图 1中4组G-γ-D曲线将用于下文中校准土体动力非线性本构模型的参数。基岩密度和剪切波速分别设为2200 kg/m³和1000 m/s。

图  1  土层剪切模量与阻尼比曲线

Figure  1.  Curves of shear modulus and damping ratio of soils

下载: 全尺寸图片 幻灯片

为了模拟地震荷载作用下土体的非线性力学行为，本文采用Anastasopoulos等^[14-15]提出的非线性运动硬化模型（Non-linear kinematic hardening model, NKH），根据该参数校正方法，通过循环单剪数值分析，对四个黏土层的参数$\lambda$进行了校准，并确定了相应的NKH模型输入参数。校正结果显示，对于四个黏土层，$\lambda$的取值范围从1/4到1/25。以土层1和2为例，图 2（a）给出了NKH模型与“目标”G-γ-D曲线的校准对比图，可见两者吻合较好，验证了采用参数的合理性，而图 2（b）和（c）分别展示了通过循环单剪试验数值分析获得的土层1和土层2的剪应力-剪应变曲线，可以反映黏土在不同剪应变幅值作用下的非线性滞回特性。

图  2  土体非线性本构模型参数标定

Figure  2.  Parameter correction for nonlinear constitutive model of soils

下载: 全尺寸图片 幻灯片

2.3   有限元数值模型

本文以典型隧道结构为研究对象，隧道埋深为15 m，为典型中埋隧道，隧道的外径为6.2 m，厚度为0.35 m。采用通用有限元软件ABAQUS建立土-隧道系统二维非线性数值分析模型。通过选择不同宽度进行参数敏感性分析，最终将有限元模型宽度确定为200 m，该宽度可保证由侧边界处反射回结构区域的地震波足够小。

在数值模型中，隧道结构采用两节点线性梁单元，周围土体采用四节点二次减缩积分平面应变单元模拟，为确保地震波在网格单位中的有效传播，根据研究的地震波频率范围（即0.5~10 Hz），确定了合理的土体网格单元大小，同时，对隧道附近的土体网格进行了局部加密处理。为了提高数值分析模型的计算效率，隧道衬砌采用了线弹性材料。隧道衬砌和土体的接触界面设置为面面接触，接触面法线方向采用硬接触，接触面的切线方向设置为摩擦系数为0.4的库仑模型。所有分析土层的阻尼比均设置为5%，并按照Hashash等^[16]提出的瑞利阻尼参数校正方法合理确定了相应的阻尼参数。

模型的边界条件设置是影响隧道地震响应的重要因素。本文在模型两侧边界上相同高度节点处设置了绑定约束，使得两侧边界在剪切波作用下具有相同的变形模式。此外，模型底部边界设定为竖向固定，水平方向采用Lysmer-Kuhlemeyer黏性边界^[17]以模拟弹性半空间下卧层吸收人工边界的反射波，边界黏性阻尼系数由下卧基岩密度ρ、基岩剪切波速${v_{\text{s}}}$和基岩与土体接触面积S确定。

本文数值分析基于排水条件下的总应力原则，数值模型计算主要分为两步：初始地应力分析和动力分析。在初始地应力分析中，模型底部采用固定边界；在随后的动力分析步中，底部边界水平方向的自由度释放，同时通过阻尼器的水平方向施加竖向传播的剪切地震波。在动力分析步计算中采用隐式动力分析步，时间增量步长选用自动步长。

3.   概率密度演化分析

3.1   隧道抗震性能评价指标

对于圆形隧道，本文采用了Koizumi^[18]提出的地震作用下的倾斜角作为隧道结构抗震性能评价指标，用于下文的隧道抗震性能分析中，系统揭示随机地震作用下的隧道结构抗震性能演化规律。该性能评价指标合理考虑了地震作用下圆形隧道实际破坏情况与隧道倾斜角之间的高度相关性。隧道倾斜角为

式中：$\eta$为隧道拱顶和拱底之间的水平相对位移；${D_0}$为地震前隧道初始直径。

王继栋^[19]通过非线性静力推覆方法对圆形隧道结构的抗震性能进行了系统分析，并合理划分了倾斜角在不同破坏状态下的阈值及其范围。根据不同的结构地震响应特征，确定了4种破坏状态，并提出了相应的破坏阈值：无破坏（φ < 1/1400），轻微破坏(1/1400≤φ＜1/600），中等破坏(1/600≤φ＜1/250），以及严重破坏(1/250≤φ＜1/150)。

3.2   隧道抗震性能概率密度演化特征

根据概率密度演化理论的隧道随机动力分析的基本步骤，将所有输入地震动下的隧道倾斜角的时程曲线带入概率密度演化公式（4），并根据概率密度演化方程组的建立与求解流程进行求解，可以快速获得隧道结构基于倾斜角指标的概率密度演化函数。图 3给出了t = 5~10 s对应的隧道倾斜角三维概率密度演化曲面。可以看出，隧道倾斜角的概率密度演化函数随时间呈现先增大后减少的趋势。通过该图可以获得任意时刻隧道倾斜角的概率密度函数，和这种概率分布随时间的变化趋势，更全面地了解地震动的不确定性对隧道抗震性能的影响机制。

图  3  典型时段内的隧道倾斜角概率密度曲面

Figure  3.  Probability density surface of tunnel inclination angle in a typical time interval

下载: 全尺寸图片 幻灯片

基于隧道倾斜角的三维概率密度演化曲面以及不同地震动作用下的隧道倾斜角最大值，可进一步通过虚拟随机过程和等效极值事件^[6]分析隧道结构在不同破坏状态下的超越概率。图 4给出了隧道倾斜角极值的概率密度函数和累计分布函数。基于图 4可以快速获得不同破坏状态下的失效概率，即无破坏、轻微破坏、中等破坏和严重破坏时倾斜角的超越概率分别为98.88%，86.10%，30.90%，3.70%。

图  4  累计分布函数曲线

Figure  4.  Curve of cumulative distribution function

下载: 全尺寸图片 幻灯片

上述分析表明，通过采用概率密度演化方法、虚拟随机过程和等效极值事件等理论，可以快速获得随机地震作用下隧道结构地震动力响应的概率信息、等价极值事件的概率密度函数（PDF）和概率分布函数（CDF），以及隧道结构在不同破坏状态下的超越概率，更加准确地评估隧道结构的抗震能力。

4.   结论

本文将概率密度演化方法引入到隧道结构的抗震性能评价中，通过基于物理机制的随机地震动模型合理考虑了地震动的随机性，进行大量土-结构非线性动力时程分析，系统揭示了随机地震作用下隧道抗震性能演化规律，得到隧道倾斜角三维概率密度演化曲面，定量给出了不同倾斜角对应的隧道超越概率。根据上述分析，可得出以下3点结论。

（1）采用基于物理机制的随机地震动模型生成了一系列非平稳地震动时程曲线。所生成的地震动平均反应谱与规范设计谱较为一致，验证了该随机地震动模型用于隧道抗震性能分析的适用性，并能合理考虑地震动的随机特性。

（2）通过提出的概率密度演化方法可以准确获取隧道结构抗震性能的量化指标（即：隧道倾斜角）的瞬时概率信息，建立的相应概率密度函数可有效体现隧道倾斜角随时间变化规律。通过隧道倾斜角瞬时概率密度函数演化，进一步揭示了土-隧道系统动力分析的不确定性传播机制。

（3）基于隧道不同破坏状态划分及其阈值体系，结合隧道倾斜角概率密度函数，可快速获得不同破坏状态下隧道的超越概率，即无破坏、轻微破坏、中等破坏和严重破坏时隧道倾角的超越概率分别为98.88%，86.10%，30.90%，3.70%。

未来的研究可以进一步完善该方法，提高模型准确性和可靠性，同时更好地考虑地震荷载的多样性和不确定性。