富水地层矿山法隧道坍塌机理与安全系数研究

Study on Collapse Mechanism and Safety Factor of Tunneling by Mining Method in Water-rich Strata

摘要

广州地铁十一号线沙河站1^#横通道在富水地层采用矿山法施工时突发坍塌，造成巨大的损失与影响。为揭示其致灾机理，以该坍塌隧道工程为背景，通过建立反映地层变化特征的三维数值模型，考虑流固耦合作用，探究了富水地层中隧道坍塌的致灾因素，并系统分析了各因素对隧道稳定性的影响。此外，基于数值模拟所揭示的隧道坍塌体特征，构建了富水条件下矿山法隧道掌子面的失稳破坏机制，结合强度折减法与极限分析上限定理，提出可以考虑地下水作用和隔水层厚度的矿山法隧道掌子面安全系数上限解的计算方法，并通过与数值模拟结果和已有文献算例对比，验证了该方法的可靠性。结果表明：拱顶的高水头差、拱顶围岩强度过低，支护不及时和作为保护层（隔水层）的微风化含砾粉砂岩厚度不足是沙河站1^#横通道拱顶失稳的主要诱因。在现场施工条件下，拱顶水头超过13 m或拱顶围岩的黏聚力低于400 kPa时，围岩的变形将超出安全许可范围；未及时支护时，拱顶保护层厚度不应小于2 m。掌子面安全系数随拱顶水头高度和开挖直径增大而降低，随保护层厚度及围岩黏聚力的增大而升高。研究成果可为类似地质条件下矿山法隧道的稳定性分析与设计优化提供参考。

Abstract

Objective The collapse of shallow-buried tunnels in water-rich strata remains a persistent challenge for underground construction in densely populated urban areas. The sudden failure of the cross-passage at Shahe Station on Guangzhou Metro Line 11, which occurred during excavation using the mining method, highlights the complex interaction among groundwater, surrounding rock strength, and excavation support conditions. This study aims to clarify the mechanism responsible for the collapse and identify the dominant factors controlling face stability in fractured, water-bearing strata. Special emphasis is placed on establishing a rational evaluation method for the face safety factor under seepage conditions by integrating numerical simulation results with upper-bound limit analysis, providing design guidance for tunnel construction in similar hydrogeological environments. Methods A three-dimensional numerical model was developed to reproduce the excavation process of the cross-passage under fluid-solid coupling conditions. The computational framework combined finite-difference-based stress-deformation analysis with groundwater seepage modeling, enabling the simulation of the coupled interaction among pore water pressure, stress redistribution, and deformation evolution. Excavation was simulated in stages based on the practical excavation sequence. Key factors were parametrically varied, including hydraulic head, strength parameters of the weathered silty sandstone, support installation timing, and the thickness of the protective layer above the tunnel crown. The strength reduction method was applied to determine the reduction factor at which the numerical model reached failure to quantify tunnel face stability. In addition, an analytical formulation of the tunnel face safety factor was proposed by integrating the strength reduction approach with the upper bound theorem of limit analysis while incorporating groundwater seepage effects. The reliability of the proposed upper bound solution was validated through comparisons with published case studies and numerical simulation benchmarks. The collapse morphology revealed by the numerical model was further utilized to establish a mechanistic model describing the instability process of the tunnel face under water-rich conditions. Results and Discussions The simulations successfully reproduced the collapse characteristics of the Shahe Station cross-passage and captured the principal failure features, including crown subsidence, progressive deformation of the surrounding rock, and hydraulic gradient-induced instability. The results revealed four dominant destabilizing factors: 1) a high hydraulic head difference at the crown, 2) excessively low strength of the crown surrounding rock, 3) delayed installation of support structures, and 4) insufficient thickness of the weakly weathered silty sandstone serving as a protective layer. Quantitative threshold analyses indicated that when the crown hydraulic head exceeded approximately 13 m or the cohesion of the crown rock mass fell below 400 kPa, face deformation exceeded the allowable safety limit. Under conditions in which timely support was not provided, the protective layer thickness was required to be at least 2 m to prevent collapse. The proposed upper-bound method for calculating the tunnel face safety factor under groundwater influence showed strong agreement with published analytical solutions and independent numerical simulation results. Sensitivity analyses demonstrated that the tunnel face safety factor decreased as the hydraulic head and tunnel excavation diameter increased, whereas it increased with greater protective layer thickness and higher rock cohesion. The comparison between the numerical and analytical results confirmed that incorporating groundwater seepage into the upper bound limit analysis was essential for accurately assessing tunnel face stability. In addition, the collapse mechanism model established from the numerical observations revealed a progressive failure process consisting of initial pore pressure accumulation at the excavation crown, localized plastic deformation within the weak sandstone, rapid loss of load-bearing capacity due to hydraulic gradient amplification, and eventual roof collapse accompanied by large-scale deformation. Conclusions The study clarified the coupled hydro-mechanical factors governing the collapse of shallow-buried tunnels in water-rich strata. The main triggers of the Shahe Station collapse were identified as high crown water head, low crown rock strength, insufficient protective layer thickness, and delayed support installation. The integration of the strength reduction method and upper bound limit analysis produced a reliable upper bound solution for evaluating tunnel face safety factors under groundwater influence. The validated calculation method revealed consistent trends in safety factor variation with hydraulic head, excavation diameter, protective layer thickness, and rock cohesion. The findings provide a scientific basis for risk evaluation and design optimization of mining-method tunnels under similar hydrogeological conditions.

本刊网刊

浅埋隧道施工不可避免地会造成围岩扰动，从而引发地层变形，严重时可能发展为地表塌陷，造成不可挽回的损失^[1‒3]。2019年12月，广州地铁十一号线沙河站1^#横通道突然坍塌，事故造成3人死亡，直接经济损失约2 004.7万元，是近年来城市地铁隧道施工过程中发生的典型重大工程事故之一。截至目前，关于该事故的分析仅有官方发布的事故调查报告，还缺乏对事故隧道失稳力学机理的系统性研究。

目前，针对隧道稳定性的研究主要有数值模拟^[4‒7]与理论分析^[8‒11]两类。在富水地层隧道掌子面稳定性与失稳机制方面，已有研究多基于考虑流固耦合的数值模拟开展。如：Li^[12]和Chen^[13]等分别分析了青海某盾构深埋隧道和温州某断裂带隧道坍塌的破坏特征及机理；Xian等^[14]揭示了地下水侵蚀作用诱发隧道突涌水的力学机制。但是，目前研究多针对单一、均质地层，而沙河站横通道位于高水压区段，其拱顶上方为裂隙发育的强风化岩体与富水砂层，沿开挖方向还存在显著的地层变化，地质条件较为复杂。综上，针对这类富水复杂地质条件下矿山法隧道的数值模拟研究还不够深入。

理论分析方面，主要有极限平衡法^[15‒18]和极限分析法^[19‒23]。相较于极限平衡法，极限分析法不依赖于对破坏面力学条件的简化假设，计算结果更为严谨，因而被学者广泛采用。黄阜等^[24]构建了双模式盾构开挖过程中前方岩土体破坏机制，利用极限分析上限定理推导出极限状态下土仓压力的计算公式。Liu等^[25]基于极限分析上限定理，考虑拱顶与侧壁联合破坏的三维塌落机制，推导了相应的极限支护力表达式。Luo等^[26]推导了非线性破坏准则下浅埋隧道围岩压力计算公式。然而，已有研究大多聚焦于盾构隧道，相比之下，矿山法隧道由于在掌子面处缺乏盾构机提供的支护力，围岩应力状态和破坏机制均与盾构隧道有明显的不同。因此，现有研究成果难以直接应用于富水条件下的矿山法隧道。除此之外，目前的理论研究多针对不考虑地下水作用的单一地层，对于在富水条件下复合地层中隧道安全系数的解析解研究仍然较为缺乏。

针对上述问题，本文以广州地铁十一号线沙河站发生坍塌事故的1^#横通道为工程背景，建立考虑地层变化的三维数值模型，探究造成坍塌的关键致灾因素及其对掌子面稳定性的影响规律，揭示富水条件下复合地层中的矿山法隧道破坏机制。在此基础上，结合极限分析上限定理与强度折减法，提出考虑地下水作用和隔水层影响的矿山法隧道掌子面安全系数的计算方法。研究成果可为富水条件下矿山法隧道的稳定性分析与安全施工提供参考。

1 工程概况与事故分析

1.1 事故概况

勘察结果显示，沙河站事故横通道所处区域地质条件复杂，事故发生段拱顶围岩为强风化与微风化岩交互分布，其中，强风化岩围岩稳定性较差，岩体基本质量等级为Ⅴ类。1^#横通道H13地质剖面如图1所示。

1.2 坍塌隧道的三维数值建模

由图1可以看出，事故横通道上覆地层类型较为复杂，局部夹杂有少许粉质黏土层及中、全风化含砾粉砂岩层等。从空间分布特征来看，上述地层多为间断性分布，厚度较薄，在隧道开挖区内并未形成对整体稳定性具有决定性影响的控制层。相比之下，填土层、中粗砂层、强风化含砾粉砂岩层及微风化含砾粉砂岩层分布广泛，对力学与渗流特性起主导作用，故本文主要考虑这4类地层建立数值计算模型。

模型尺寸为100 m（x向）×100 m（y向）×65 m（z向），包含618 223个单元和111 051个节点，如图2（a）所示。地下水条件按线性分布设定，模型底部孔隙水压力为6.3×10⁵ Pa，梯度为-1×10⁴ Pa/m，即模型拱顶初始水头高度设为19 m。模型底部施加孔压和位移固定边界，模型四侧施加法向位移约束。实际施工时，上台阶最后1 m爆破开挖后，因喷浆设备故障，此部分支护结构未完成，因此数值模型中该段亦不设支护。值得注意的是，地质钻孔信息表明，坍塌区附近拱顶围岩由微风化含砾粉砂岩变为强风化含砾粉砂岩层，如图2（b）所示。

根据事故前后勘察工作中的土工试验结果，土层及支护结构的物理力学参数取值见表1。

1.3 隧道开挖过程及坍塌特征分析

实际施工中，上台阶预计推进50.6 m。在隧道掘进前49.6 m的过程中，初期支护按规范同步实施，未见突水、突泥等异常，围岩整体稳定。最后1 m爆破开挖后，由于喷浆机故障，初期支护未能及时完成，该无支护状态持续约10 h后，掌子面拱顶突现涌水且破坏在3 min内扩展至地表，导致地面塌陷。

图3为数值计算中上台阶开挖至49.6和50.6 m时拱顶竖向位移云图。由图3可以看出：当开挖进尺至49.6 m时，数值计算收敛，拱顶最大沉降为16.5 mm，隧道底部隆起11.2 mm；当开挖延伸至50.6 m且最后1 m处不设支护时，沉降值呈持续增大的不收敛趋势，最大沉降集中在未支护段拱顶，表明围岩显著失稳，在拱顶未支护区形成局部塌腔并沿最薄弱路径迅速贯通至地表。

图4（a）为地表的沉降云图。由图4（a）可以看出，最大沉降发生在掌子面上方，沉降区呈不规则圆形，直径约28 m，与事故现场（图4（b））的实测塌陷范围高度一致，表明所建立的数值模型能较好地反映实际情况。

1.4 关键致灾因素分析

富水地层中，渗流作用导致的隧道坍塌事故屡见不鲜。因此首先考察地下水在事故中的作用。保持其他条件不变，不考虑地下水作用时，开挖进尺至50.6 m的位移云图如图5所示。与图3（b）中考虑地下水作用的计算结果对比可见，此时数值计算可顺利收敛，沉降区分布相对均匀，最大沉降仅为12.9 mm，明显低于考虑地下水作用时的沉降值，隧道处于较为稳定状态。由此可见，地下水的存在显著加剧了隧道拱顶沉降和围岩失稳，是导致隧道坍塌的关键因素。

此外，图2（b）的土层分布显示，沿隧道掘进方向推进，拱顶上方微风化含砾粉砂岩层厚度逐步减小，最终过渡为强风化含砾粉砂岩层。结合拱顶竖向位移分布特征可以发现，随着拱顶微风化层厚度的减小，拱顶沉降值逐步增大，说明微风化层结构完整、力学性能良好，具备较强的承载能力，在开挖过程中有效分担了围岩应力集中，起到了保护层作用。而当拱顶揭露强风化含砾粉砂岩层后，因该岩层风化程度高、结构松散、黏聚力较低，在初期支护不及时的情况下，拱顶岩体无法有效抵抗开挖引起的应力扰动，极易诱发围岩失稳。当拱顶刚揭露强风化层围岩时，由于支护及时，隧道并未失稳；而在最后1 m的隧道掘进后未及时支护，隧道最终坍塌，因此，不及时支护也是该事故发生的主要诱因之一。综上，地下水作用、强风化含砾粉砂岩强度过低、支护不及时，以及拱顶保护层微风化层厚度不足，是本次坍塌事故的关键。

2 致灾因素对围岩变形的影响规律

在工程设计中通常将隧道的预留变形量作为其极限位移控制值^[27]。根据《铁路隧道设计规范》（TB 10003—2016）^[28]中复合式衬砌条件下不同级别围岩的变形控制建议值，结合1^#横通道拱顶所揭露的Ⅴ类强风化含砾粉砂岩层特性，以及工程处于浅埋、中跨隧道的实际情况，确定拱顶最大容许沉降为120 mm。因此，当计算位移超出容许沉降时，可认为隧道有失稳风险。在此基础上，量化分析水位高度、围岩强度参数及保护层厚度3个致灾因素对隧道稳定性的影响。

2.1 拱顶上方水位高度

为探究拱顶水位高度对隧道稳定的影响，保持其他参数不变，设置拱顶上方水头高度h_w分别为3、6、9、12和15 m的5组工况，结果如图6所示。

由图6可以看出，随水头高度增加：拱顶最大沉降值由22增加至146 mm；地表沉降槽宽度随之增加，其峰值由10.52增加至16.33 mm；掌子面位移由39.0增加至68.8 mm，水平收敛量由13.9增加至32.1 mm。当水头高度在13 m左右时，拱顶最大沉降值接近120 mm的规范容许极限，说明隧道结构已接近失稳临界状态。因此，当拱顶上方水头超过13 m时，应高度重视并及时采取降排措施。

2.2 强风化砂砾岩层强度参数

为模拟拱顶强风化层风化程度对隧道稳定性的影响，设置黏聚力c分别为1 200、1 000、800、600、400和200 kPa的6组工况，计算结果如图7所示。

由图7可知，随着c由1 200降至200 kPa，拱顶最大沉降值由42增加至158 mm，地表沉降量峰值由15.21增加至20.13 mm，掌子面位移由56.2增加至77.4 mm，水平收敛量由20.3增加至40.1 mm。当强风化层黏聚力与微风化层相同（1 200 kPa）时，拱顶最大沉降值仅为42 mm；而当黏聚力降至400 kPa时，拱顶最大沉降为125 mm，已超过规范容许值，表明围岩变形已达失稳边界。因此，在施工过程中，当超前钻孔显示拱顶上方揭露岩层风化程度过高时，应及时采取优化支护参数等控制措施，确保隧道安全施工。

2.3 保护层厚度

为研究保护层对隧道稳定性的影响，根据坍塌处剖面的地质情况，保持中粗砂层与强风化层厚度不变，将数值模型简化为均质水平层状地层进行计算。设置微风化保护层厚度d_s分别为0.5、1.0、1.5、2.0、2.5和3.0 m的6组工况，此时隧道整体位于微风化含砾粉砂岩层中，综合岩体特性与工程实际情况，根据规范^[28]确定拱顶极限沉降值为80 mm。

图8为拱顶最大沉降值随保护层厚度的变化曲线。由图8可知，随着保护层厚度由0.5增加至3.0 m，拱顶最大沉降值由215下降至68 mm。当保护层厚度小于2.0 m时，拱顶沉降值超过规范容许值且对厚度变化较为敏感；当保护层厚度增加到一定值后，随着其厚度进一步增大，对拱顶沉降的影响逐渐减弱。

不同保护层厚度下的围岩塑性区分布如图9所示。

由图9可见，开挖后拱顶上方形成“V”形塑性区且向地表扩展，其范围随保护层厚度增大而明显缩小。当保护层厚度小于2.0 m时，该区与上覆强风化层的塑性区贯通，形成连续破坏通道；当厚度大于等于2.0 m时，拱顶塑性区显著收缩，与强风化层之间逐渐形成分离趋势，这表明较厚的保护层有助于阻隔破坏发展路径，提升拱顶稳定性。

3 掌子面安全系数上限解

相较于数值模拟方法，理论方法无需建立复杂的三维数值模型即可快速便捷地判断隧道的稳定性。考虑到矿山法隧道掌子面不存在支护力的特点，可通过极限分析上限定理计算掌子面极限支护力，再根据强度折减法折减岩土层强度参数，得到恰好使得极限支护力为0时的折减系数，即为掌子面安全系数，以进一步量化分析矿山法隧道在富水地层中的稳定性。

3.1 破坏机制及安全系数求解过程

由于极限分析上限定理中的能量耗散计算依赖于预设的破坏机制，因此构建能够真实反映隧道破坏特征的合理破坏机制，是获得准确上限解的关键所在。基于相关研究成果^[29‒30]，结合本文数值模拟结果中隧道拱顶围岩塑性区呈现的“V”形贯通破坏特征，构建隧道拱顶上方岩土体在极限状态下的破裂机制模型，如图10所示。

沙河站1^#横通道上台阶开挖直径与开挖高度之比约为2，据此可将隧道上台阶开挖形状简化为半圆形。图10中：H为隧道埋深；B为隧道开挖直径；α为破裂面与水平面的夹角，\alpha =\frac{\mathrm{\pi }}{\mathrm{4}}+\frac{\phi }{\mathrm{2}}；d₀～d₄分别为破裂面与隧道轮廓相交处至隧道底部的高度、破裂面与隧道相交处至微风化含砾粉砂岩层顶的高度、强风化含砾粉砂岩层厚度、中粗砂层厚度和填土层厚度，d_{\mathrm{0}}=\frac{B}{\mathrm{2}}\cdot \mathrm{s}\mathrm{i}\mathrm{n}\alpha，且有H+\frac{B}{\mathrm{2}}={\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{4}}}{d}_i；\dot{v}为机动允许的速度场，方向竖直向下；τ_n和σ_n分别为土体的剪应力和正应力；d_s为保护层厚度，d_{\mathrm{s}}=H-d_{\mathrm{2}}-d_{\mathrm{3}}-d_{\mathrm{4}}；D₀～D₄分别为拱顶破裂面宽度、微风化层顶破裂面宽度、强风化层顶破裂面宽度、中粗砂层破裂面宽度、填土层破裂面宽度，D_i=\frac{\mathrm{2}\cdot {\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{4}}}{d}_i}{\mathrm{t}\mathrm{a}\mathrm{n}\alpha }；l(x)为破裂面曲线方程；L为破裂面长度，L=\frac{(H+\frac{B}{\mathrm{2}})}{\mathrm{s}\mathrm{i}\mathrm{n}\alpha }-\frac{B}{\mathrm{2}}。

根据极限分析上限定理，当结构处于极限状态时，外力所做功的功率应等于内部能量的耗散率，且速度间断线上任一点的内能耗散率\dot{D}可由法向正应力和切向应力计算。因此，有如下关系：

{\int }_V{\mathit{\sigma }}_{ij}{\dot{\mathit{\epsilon }}}_{ij}\mathrm{d}V={\int }_S{T}_i{\mathit{v}}_i\mathrm{d}S+{\int }_V{X}_i{\mathit{v}}_i\mathrm{d}V（1）

\dot{D}={\sigma }_{\mathrm{n}}{\dot{\epsilon }}_{\mathrm{n}}+{\tau }_{\mathrm{n}}{\dot{\gamma }}_{\mathrm{n}}（2）

式（1）、（2）中，T_i 和X_i 分别为作用于岩土体表面的边界力与内部的体积力， v_i 为机动许可的速度场矢量，S为速度场边界， σ_ij 为应力张量，{\dot{\mathit{\epsilon }}}_{ij}为应变率张量，V为研究区域体积，{\dot{\epsilon }}_{\mathrm{n}}和{\dot{\gamma }}_{\mathrm{n}}为法向应变率和法向剪应变率。

采用Mohr‒Coulomb强度准则描述岩土体剪切破坏行为，其对应的屈服函数f可表示为：

f={\tau }_{\mathrm{n}}-c-{\sigma }_{\mathrm{n}}\mathrm{t}\mathrm{a}\mathrm{n}\phi（3）

在塑性理论中，流动规则用以确定塑性应变增量的方向或塑性应变增量张量的各个分量之间的比例关系，塑性应变增量的方向由应力空间的等塑性势面g决定，两者关系如下：

{\dot{\epsilon }}_{ij}^{\mathrm{p}}=\dot{\lambda }\frac{\partial g}{\partial {\mathit{\sigma }}_{ij}}（4）

式中，{\dot{\epsilon }}_{ij}^{\mathrm{p}}为塑性应变率，\dot{\lambda }为一个贯穿于整个塑性加载历史的非负标量函数，梯度矢量\frac{\partial g}{\partial {\mathit{\sigma }}_{ij}}规定了塑性应变增量矢量的方向。

在极限分析中，通常假定材料为刚塑性体，弹性变形可忽略，因此应变率可视为塑性应变率，即{\dot{\epsilon }}_{ij}={\dot{\epsilon }}_{ij}^{\mathrm{p}}。

如果塑性势能面g与屈服面f是重合的，即g=f，说明流动法则与屈服准则是相关联的，被称为相关联的流动法则。据此建立与Mohr‒Coulomb强度准则相关联的流动法则，可以推导得出正应变率{\dot{\epsilon }}^{\mathrm{p}}与剪应变率{\dot{\gamma }}^{\mathrm{p}}之间的关系如式（5）所示：

\frac{{\dot{\epsilon }}^{\mathrm{p}}}{{\dot{\gamma }}^{\mathrm{p}}}=-\mathrm{t}\mathrm{a}\mathrm{n}\phi（5）

根据应变局部化理论，可以将土体破裂面视为具有较小厚度w的变形层，基于Mohr‒Coulomb 强度准则及其相关联流动法则，结合图10所示的几何关系，法向应变率大小可表示为：

{\dot{\epsilon }}_{\mathrm{n}}=\frac{v}{w\sqrt[]{\mathrm{1}+{\left|l^{\text{'}}(x)\right|}^{\mathrm{2}}}}（6）

式中：l^{\text{'}}(x)为破裂面曲线方程的斜率，l^{\text{'}}(x)=\frac{\mathrm{1}}{\mathrm{t}\mathrm{a}\mathrm{n}\phi }；v表示速度间断面两侧的相对速度。

将式（5）和（6）代入式（2）中，岩土体变形破坏产生的内能耗散率\dot{D}最终可表示为：

\dot{D}=c\frac{v}{w}\mathrm{c}\mathrm{o}\mathrm{s}\phi（7）

拱顶围岩向下发生失稳破坏时其内能耗散仅仅集中于速度间断线L上，根据破坏机制，围岩失稳时的内能耗散功率可表示为：

W_{\dot{D}}={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{W_{\dot{D}}}_i={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}\mathrm{2}{c}_i{\mathit{v}}_i{L}_i\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{a}\mathrm{v}\mathrm{g}}=\mathrm{2}{c}_{\mathrm{a}\mathrm{v}\mathrm{g}}vL\mathrm{c}\mathrm{o}\mathrm{s}{\phi }_{\mathrm{a}\mathrm{v}\mathrm{g}}（8）

式中：c_avg为综合等效黏聚力，c_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{c}_i{d}_i}{{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{d}_i}；{\phi }_{\mathrm{a}\mathrm{v}\mathrm{g}}为综合等效内摩擦角，{\phi }_{\mathrm{a}\mathrm{v}\mathrm{g}}=\frac{{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\phi }_i{d}_i}{{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{d}_i}；c_i 和φ_i 分别为对应第i（i=1,2,3, …）层土的黏聚力和内摩擦角；W_{\dot{D}i}为第i层土失稳时的内能耗散功率；L_i 为第i层土的破裂面长度。

外力做功主要包括隧道上方土体自重做功、水压力做功及隧道内部支护力做功这3个部分。其中，土体自重做功功率W_γ 为：

\begin{array}{l}{W}_{\gamma }={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{W_{\gamma }}_i={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}\left[{\gamma }_i\frac{\left(D_{i-\mathrm{1}}+D_i\right)}{\mathrm{2}}{d}_i\right]v-\\ {\gamma }_{\mathrm{1}}\left[\alpha {\left(\frac{B}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}{D}_{\mathrm{0}}{d}_{\mathrm{0}}\right]v\end{array}（9）

式中，{W_{\gamma }}_i为第i层土做功功率，γ_i 为各土层的重度。

微风化含砾粉砂岩地下水赋存条件相对较差，一般具微～弱透水性，富水性弱，因此，将其视为隔水保护层，则上方附加水压力做功功率W_w为：

W_{\mathrm{w}}={\gamma }_{\mathrm{w}}(h_{\mathrm{w}}-d_{\mathrm{s}})D_{\mathrm{1}}v（10）

式中，γ_w为水的重度。

支护力做功功率W_q为：

W_q=-q{D}_{\mathrm{0}}v（11）

式中，q为隧道内支护力。

由虚功率原理可知，从能量转换和平衡的角度出发，令外力所做的功率与岩体内部的能量耗散率相等，将式（8）、（9）、（10）与（11）代入式（12）中：

W_{\dot{D}}=W_{\gamma }+W_{\mathrm{w}}+W_q（12）

可以进一步求得极限状态下的支护力q_cr，如式（13）所示：

\begin{array}{l}{q}_{\mathrm{c}\mathrm{r}}=\left\{{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}\left[{\gamma }_i\frac{\left(D_{i-\mathrm{1}}+D_i\right)}{\mathrm{2}}{d}_i\right]\right.-\\ {\gamma }_{\mathrm{1}}\left[\alpha {\left(\frac{B}{\mathrm{2}}\right)}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}{D}_{\mathrm{0}}{d}_{\mathrm{0}}\right]+\\ \left.{\gamma }_{\mathrm{w}}(h_{\mathrm{w}}-d_{\mathrm{s}})D_{\mathrm{1}}-\mathrm{2}{c}_{\mathrm{a}\mathrm{v}\mathrm{g}}L\mathrm{c}\mathrm{o}\mathrm{s}\phi \right\}/D_{\mathrm{0}}\end{array}（13）

根据强度折减法，将原始的抗剪强度指标综合等效黏聚力c_avg和综合等效内摩擦角φ_avg按式（14）和（15）折减后代入式（13），并使折减后隧道极限支护力等于矿山法隧道施工过程中的实际支护力，即q_cr=q=0，此时的折减系数F_s即为安全系数的最优上限解。

c_{\mathrm{f}}=\frac{c_{\mathrm{a}\mathrm{v}\mathrm{g}}}{F_{\mathrm{s}}}（14）

{\phi }_{\mathrm{f}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(\frac{\mathrm{t}\mathrm{a}\mathrm{n}\phi }{F_{\mathrm{s}}}\right)（15）

式（14）、（15）中，c_f和φ_f分别为折减后的黏聚力和内摩擦角。

3.2 与已有计算方法的对比分析

由于目前关于矿山法隧道在富水条件下的安全系数研究较少，无法直接与已有成果进行对比。为验证本文方法的可靠性，选取考虑地下水作用的盾构隧道^[31]和不考虑地下水作用的矿山法隧道^[32‒33]的掌子面安全系数计算方法进行对比，最后再结合本文关于沙河站横通道的数值模拟结果，综合说明本文方法的准确性，结果见表2。由表2的对比结果可知，本文基于极限分析上限定理所求得的安全系数与文献及数值模拟计算结果均较为接近，误差最大为5.7%，验证了本文方法的合理性。造成误差的原因主要在于，黄阜等^[31]构建的破坏机制中破坏范围并未贯穿至地表，土体自重做功功率偏小，得到的安全系数更大。陈峥等^[32]与王海涛^[33]的方法考虑了锚杆和管棚的支护作用，因此安全系数均更大。本文计算方法针对复杂地层条件下贯穿地表塌陷的破坏形式，更接近现场破坏情况，因此更适用于对类似富水条件下矿山法隧道的稳定性评估。

3.3 算例分析

以沙河站1^#横通道的参数为背景，取隧道直径B=10～14 m，拱顶上方水位高度h_w=0～20 m，保护层黏聚力c=200～1 000 kPa，保护层厚度d_s=0～2.5 m，依托本文提出的安全系数计算方法，开展了多组算例分析，系统探讨了上述参数对隧道掌子面安全系数的响应规律，如图11所示。由图11可知，安全系数随着隧道开挖直径及水头高度的增加而逐渐减小，随着拱顶上方保护层黏聚力和厚度的增加而增大。

以图11（a）中接近施工情况的开挖直径B=12 m的算例为例，分析安全系数随拱顶上方水头高度变化的规律。当水头高度从0增加至20 m时，安全系数由1.21降低至0.85，呈持续下降趋势，其中：当水头高度达到10 m时，安全系数降至1.03，接近临界稳定状态；根据曲线的变化规律，当水头高度进一步增加至12 m时，安全系数将降至1以下，该结果与第2节数值模拟中拱顶上方水头高度超过 12 m时隧道趋于失稳的判断基本一致。

在图11（b）中B=12 m条件下，当黏聚力由200增加至1 000 kPa时，安全系数由0.81增加至1.48，呈显著增长趋势。其中，当黏聚力为400 kPa时，安全系数为0.97，低于稳定判别阈值，与数值模拟中拱顶岩土体黏聚力低于400 kPa时隧道趋于失稳的结论基本一致。

在图11（c）中B=12 m条件下，随着保护层厚度由0增加至2.5 m，安全系数由0.71增加至1.10；当保护层厚度为2.0 m时，安全系数达到1.02，表明拱顶围岩处于相对稳定状态，与第2节数值模拟中当保护层厚度小于2.0 m时隧道处于临近失稳状态的判断相一致。

4 结 论

以广州地铁十一号线沙河站1^#横通道的坍塌为背景，采用数值模拟与理论分析相结合的方法，系统研究了隧道坍塌破坏机制及稳定性影响因素，建立了考虑地下水作用的矿山法隧道掌子面安全系数上限解计算方法。主要研究结论如下：

1）在地下水作用、强风化含砾粉砂岩层强度过低及不支护条件下，拱顶上方保护层厚度不足是沙河站1^#横通道坍塌事故的关键致灾因素。数值模拟结果表明，当隧道开挖至事故段时，拱顶围岩完全揭露强风化含砾粉砂岩层且初期支护不及时，在地下水的渗流作用下引发拱顶围岩失稳，塑性区贯通至地表，模拟结果与实际塌陷特征高度一致。

2）随着水头高度升高、强风化层黏聚力降低，围岩变形明显加剧。现场条件下，当水头高度超过13 m或岩体黏聚力降至400 kPa以下时，拱顶最大沉降值超过规范要求的120 mm容许极限，隧道失稳风险较大。保护层厚度的增加对围岩稳定性具有显著提升作用，当厚度小于2.0 m时，塑性区自拱顶向地表贯通，围岩失稳；当厚度大于等于2.0 m时，拱顶塑性区范围有效减小，沉降也显著降低。因此，施工过程中应综合采取降排水、加强支护等措施，有效提升隧道的稳定性。

3）提出的隧道安全系数计算方法所得结果与数值计算和已有文献结果对比均较为吻合，最大误差仅为5.7%，表明了本文方法的可靠性。算例分析结果表明，隧道安全系数随着水头高度和开挖直径的增加而降低，随着拱顶保护层黏聚力和厚度的增加而增大。

需要说明的是，在计算掌子面安全系数上限解时，由于隔水层的存在，本文对拱顶上方水压力的影响分析是基于静水压力的分布假设，未考虑渗透系数的空间差异、降雨或突水工况下的渗流效应，所得结果主要反映了稳定地下水条件下水头变化对隧道稳定性的影响规律。在复杂水文地质环境中，若存在突涌水、降雨等情况，水压力分布特征及其对隧道失稳的作用机制可能更为复杂，需结合实际工况进一步开展更精细的渗流‒力学耦合分析。

参考文献

[1] Li Ao, Zhang Dingli, Huang Jun,et al.Progressive collapse mechanism and safety control of deep buried tunnel vault in soft and broken surrounding rock[J].Advanced Engineering Sciences,2022,54(6):85‒96.

[2] 李奥,张顶立,黄俊,等.软弱破碎围岩深埋隧道拱顶渐进性塌方机理及控制[J].工程科学与技术,2022,54(6):85‒96.

[3] Zhou Guannan, Yang Tengtian, Sun Ze,et al.Investigation of quantitative evaluation method and engineering application of shallow buried tunnel face stability[J].Applied Sciences,2022,12(13):6656. doi:10.3390/app12136656

[4] Liu Xingzong, Liu Fei, Song Kezhi.Mechanism analysis of tunnel collapse in a soft-hard interbedded surrounding rock mass:A case study of the Yangshan Tunnel in China[J].Engineering Failure Analysis,2022,138:106304. doi:10.1016/j.engfailanal.2022.106304

[5] Wang Jianxiu, Cao Ansheng, Wu Zhao,et al.Numerical stability analysis of large section tunnels using the double-side heading method:A case study of Xiamen Haicang evacuate-channel[J].Applied Sciences,2023,13(1):406. doi:10.3390/app13010406

[6] Cao Yang, Lin Xiangrong, Li Zilu.Refined numerical simulation of shallow shield construction in soft soil stratum[J].Advanced Engineering Sciences,2022,54(3):149‒158. doi:10.15961/j.jsuese.202100187

[7] 曹洋,林向荣,李子路.软土地层浅埋盾构施工的精细化数值模拟[J].工程科学与技术,2022,54(3):149‒158. doi:10.15961/j.jsuese.202100187

[8] Huang Feng, Shi Xiaoxiong, Wu Chuangzhou,et al.Stability analysis of tunnel under coal seam goaf:Numerical and physical modeling[J].Underground Space,2023,11:246‒261. doi:10.1016/j.undsp.2022.12.006

[9] Wang Chengwen, Liu Xiaoli, Yao Wenli,et al.Study on rock stability and structural response of shield-driven twin tunnels crossing fault fracture zone based on 3d numerical simulation[J].Bulletin of Engineering Geology and the Environment,2024,83(7):287. doi:10.1007/s10064-024-03777-9

[10] Liang Qiao, Xu Junjie, Wei Yuanguo.Three-dimensional stability analysis of tunnel face based on unified strength theory[J].Scientific Reports,2023,13:12326. doi:10.1038/s41598-023-39554-z

[11] Ding Xiaobin, Li Kang, Xie Yuxuan,et al.Face stability analysis of large shield-driven tunnel in rock-soil interface composite formations[J].Underground Space,2022,7(6):1021‒1035. doi:10.1016/j.undsp.2022.01.007

[12] Si Xuefeng, Zhang Zilong, Li Xibing,et al.Influences of maximum principal stress direction and cross-section sha-pe on tunnel stability[J].Journal of Rock Mechanics and Geotechnical Engineering,2025,17(4):2159‒2180. doi:10.1016/j.jrmge.2024.10.003

[13] Liu Yingnan, Lei Huayang, Ma Changyuan,et al.Three-dimensional stability analysis of shied tunnel face adjacent to existing tunnels[J].Chinese Journal of Geotechnical Engineering,2023,45(7):1374‒1383. doi:10.11779/CJGE20220439

[14] 刘英男,雷华阳,马长远,等.近接既有隧道盾构开挖面三维稳定性分析[J].岩土工程学报,2023,45(7):1374‒1383. doi:10.11779/CJGE20220439

[15] Li Ziqi, Lai Jinxing, Wang Meng,et al.Collapse mechanism and treatments of a deep tunnel in the weathered granite fault zone[J].Tunnelling and Underground Space Technology,2024,152:105891. doi:10.1016/j.tust.2024.105891

[16] Chen Longlong, Wang Yaqiong, Wang Zhifeng,et al.Characteristics and treatment measures of tunnel collapse in fault fracture zone during rainfall:A case study[J].Engineering Failure Analysis,2023,145:107002. doi:10.1016/j.engfailanal.2022.107002

[17] Xian Meiting, Xu Zengguang, Chen Bin.Numerical investigation of variable-mass seepage mechanism of broken rock mass in faults[J].Computers and Geotechnics,2023,160:105527. doi:10.1016/j.compgeo.2023.105527

[18] Cheng Cheng, Chen Yang, Zhao Chenyang,et al.Theoretical analysis of the shield tunnel face stability in dry sandy strata[J].European Journal of Environmental and Civil Engineering,2023,27(2):712‒732. doi:10.1080/19648189.2022.2062616

[19] Sun Yang, Che Yufei, Xu Zhenyang,et al.Soil destabilization analysis of shield tunnel excavation surface in composite strata considering arch effect[J].Chinese Journal of Underground Space and Engineering,2024,20(3):978‒986.

[20] 孙阳,车宇飞,徐振扬,等.考虑拱效应的复合地层盾构隧道开挖失稳分析[J].地下空间与工程学报,2024,20(3):978‒986.

[21] Song Qilong, Su Dong, Zhang Runqi,et al.Limit equilibrium model for active failure of large-size rectangular tunnel face considering effect of height-width ratio[J].Tunnelling and Underground Space Technology,2024,147:105712. doi:10.1016/j.tust.2024.105712

[22] Su Dong, Song Qilong, Lin Xingtao,et al.Limit equilibrium models for active failures of large-diameter shield tunnel faces in soft clay reinforced with soil-cement walls[J].Co-mputers and Geotechnics,2023,153:105104. doi:10.1016/j.compgeo.2022.105104

[23] Li Qingdong, Ning Zexu, Liu Jia,et al.Stability and economic evaluation of multi-step horizontal salt Caverns with different step distances in bedded salt formations[J].Journal of Energy Storage,2023,57:106192. doi:10.1016/j.est.2022.106192

[24] Hou C T, Yang X L.3D stability analysis of tunnel face with influence of unsaturated transient flow[J].Tunnelling and Underground Space Technology,2022,123:104414. doi:10.1016/j.tust.2022.104414

[25] Li Chuitian, Liu Yanhui, Lu Qingbi,et al.A calculation method for ultimate support force of excavation face of shield tunnels in water-rich composite strata based on upper limit method[J].Tunnel Construction,2024,44(S1):187‒196.

[26] 李垂忝,刘彦辉,卢清碧,等.基于极限上限法的富水复合地层盾构隧道开挖面极限支护力计算方法[J].隧道建设(中英文),2024,44(S1):187‒196.

[27] Guo Jiaqi, Wu Wenlong, Liu Xiliang,et al.Theoretical analysis on safety thickness of the water-resistant rock mass of karst tunnel face taking into account seepage effect[J].Geotechnical and Geological Engineering,2022,40(2):697‒709. doi:10.1007/s10706-021-01916-7

[28] Wang Wenpo, Deng Ronggui, Liu Hongfeng.Stability analysis of 3D tunnel face of shallow rectangular shield tunnel[J].KSCE Journal of Civil Engineering,2023,27(3):1368‒1382. doi:10.1007/s12205-023-1511-0

[29] Huang Fu, Chen Jingjing, Wang Yongtao,et al.Study on the stability of double-mode shield excavation face in composite strata based on limit analysis theorem[J].Journal of Civil and Environmental Engineering,2024,46(6):89‒99.

[30] 黄阜,陈晶晶,王勇涛,等.基于极限分析理论的复合地层中双模盾构开挖面稳定性研究[J].土木与环境工程学报(中英文),2024,46(6):89‒99.

[31] Liu Zhizhen, Cao Ping, Wang Fei,et al.Three-dimensional upper bound limit analysis of tunnel stability with an extended collapse mechanism[J].KSCE Journal of Civil Engineering,2022,26(12):5318‒5327. doi:10.1007/s12205-022-2065-2

[32] Luo Wei, Xiao Gequan, Tao Zhi,et al.Nonlinear stability analysis of shallow-buried bias tunnel based on failure mode improvement[J].Applied Sciences,2025,15(6):3153. doi:10.3390/app15063153

[33] 中华人民共和国交通运输部. 公路隧道施工技术规范:JTG/T 3660—2020 [S].北京:人民交通出版社,2020.

[34] 国家铁路局. 铁路隧道设计规范(2024、2026年局部修订):TB 10003—2016 [S].北京:中国铁道出版社,2017.

[35] Xie Jiaxiao.Discussion on critical depth and load of deep and shallow tunnels[J].Railway Standard Design,1986,30(11): 28‒39.

[36] 谢家烋.深、浅埋隧道临界深度及其荷载探讨[J].铁道标准设计通讯,1986,30(11):28‒39.

[37] Qu Xing, Li Ning.Analysis of calculation method of vertical pressure in loose rock mass and research on dividing line standard for deep- and shallow-buried tunnels[J].Chinese Journal of Rock Mechanics and Engineering,2011,30(S1):2749‒2757.

[38] 曲星,李宁.松散岩体竖向压力计算方法剖析及隧洞深浅埋划分方法研究[J].岩石力学与工程学报,2011,30(S1):2749‒2757.

[39] Huang Fu, Pan Qiujing, Ling Tonghua.Upper bound analysis of factor of safety for shield tunnel face subjected to underground water seepage[J].Chinese Journal of Geotechnical Engineering,2017,39(8):1461‒1469. doi:10.11779/CJGE201708013

[40] 黄阜,潘秋景,凌同华.地下水渗流作用下的盾构隧道开挖面安全系数上限分析[J].岩土工程学报,2017,39(8):1461‒1469. doi:10.11779/CJGE201708013

[41] Chen Zheng, He Ping, Yan Dumin,et al.Upper-bound limit analysis of tunnel face stability under advanced support[J].Rock and Soil Mechanics,2019,40(6):2154‒2162. doi:10.16285/j.rsm.2018.1903

[42] 陈峥,何平,颜杜民,等.超前支护下隧道掌子面稳定性极限上限分析[J].岩土力学,2019,40(6):2154‒2162. doi:10.16285/j.rsm.2018.1903

[43] Wang Haitao.Research on mechanism of pipe roof reinforcement and tunnel face stability[D].Dalian:Dalian University of Technology,2009.

[44] 王海涛.隧道管棚预支护体系的力学机理与开挖面稳定性研究[D].大连:大连理工大学,2009.