黄 金，赵 淼

(1.智性纤维复合加固南通有限公司，江苏 南通 226000；2.北京城建中南土木工程集团有限公司，北京 100000)

隧道工程的风险分析有各种各样的方法，但是通常我们采用概率法并以试验、统计资料以及专家评定作为基础，对于不同的危害或事故的发生概率进行分析估计[1].为了进行实例隧道工程的围岩风险分析，我们需要先对隧道围岩进行模拟分析计算，得到相关参数变化规律，再利用概率分析工具对隧道围岩失效概率进行计算.

赵旭峰、王春苗[2]等人通过典型工程和实测信息对岩石力学理论分析的结果进行类比、修正，对服务隧道进行不断反演迭代，逐次修正待预测参数，直到满足要求，反演得到隧道初步支护后的等效侧压力系数与等效弹性模量.并利用反演参数对主隧道的围岩变形、破坏特性与支护效果进行快速分析及超前预报，实现隧道围岩稳定性动态反馈.卫军、李昊[3]等人利用有限元软件ABAQUS中基于Rankine准则的Molar Coulomb模型对四渡河悬索桥隧道锚进行了三维弹塑性有限元分析，模拟了初始地应力场、隧道锚与公路隧道的新奥法施工过程以及隧道锚隧洞周边接触软弱层，主要以锚碇围岩的稳定性和锚固系统的整体安全度作为主要研究内容.CAI M, KAISER P K[4]等人采用GSI系统来对围岩接壤处的残余应力进行分析决策，得出一些相关参数之间的数值关系，进而将其量化拟合分析.但是这个系统对于参数的输入就不像上面的那么有任意性.WU X, JIANG Z[5]等人采用了数值概率分析中的贝叶斯网络理论对隧道工程风险进行了分析决策，但是计算方法繁琐，不直观，在实际工程中难以广泛应用.

根据国内外研究现状，本文将采用隧道专用有限元分析软件PLAXIS 3D TUNNEL和贝叶斯网络工具Netica来建立隧道围岩失稳破坏的模型，系统的完成隧道围岩结构失稳风险分析，为实例隧道开挖工程项目提供数据支撑及技术支持.

1 贝叶斯网络简介与软件选择

1.1 贝叶斯网络简介

贝叶斯网络是描述数据变量之间依赖关系的一种图形模式，是一种用来进行推理的模型.在具体应用时，贝叶斯网络结构的构造是要根据具体问题的知识原理，将直接相互依赖的变量之间用有向线段连结，形成网络结构；贝叶斯网络的条件概率分布表可以通过专家指定或是通过学习得到[6].

贝叶斯决策就是在情况不完全已知下，对部分未知的状态用主观概率估计，然后用贝叶斯公式对发生概率进行修正，最后再利用期望值和修正概率做出最优决策.

1.2 贝叶斯网络工具—Netica

目前比较广泛的贝叶斯网络工具有MATLAB BNT工具箱、Bayesbuilder、Weka3.6、Winbugs、WinMineToolkit等等[6].但是这些软件对于一个初次接触贝叶斯网络的学者来说入门是比较难的.本文中将使用一个新型且相对比较简单的贝叶斯软件Netica来进行相应的建模分析.

Netica的工作原理相当简单，基于概率分析的贝叶斯网络，通过建立结构层次分析图(即各个影响因素)，并在贝叶斯网络结构中加入各个因素各自的概率分布规律，最终通过一个功能函数进行联结，然后将函数赋予整个网络结构，这样就可以进行最后的取样计算分析[7].

该软件建立贝叶斯网络结构模型的过程相对于其他贝叶斯网络工具来说，是比较简单的.它在风险分析决策方面的应用是极其广泛的，通过一个简单的贝叶斯网络图以及影响因素的概率分布，即可计算出在各因素影响下该项目的失效概率.通过改变单一条件，该工程项目的失效概率也随之变化，更有利于进行分析决策.

2 隧道风险参数计算分析

2.1 隧道衬砌轴应力计算

本文研究的实例隧道是连云港主体港区东疏港高速公路隧道.该隧址处基岩为中-微风化云台组浅粒岩、变粒岩，为硬质岩石，致密坚硬，抗压强度大，抗风化性能好，岩体完整性较好；局部地段分布有条带状绿泥石片岩，为软质岩，中-微风化；该处山体呈单面山构造；隧道工程区内未见断裂构造带，但局部分布有节理密集带，结构面发育，主要有两组，多属闭合-微张.在隧址区地面地质调绘中，中-微风化变粒岩呈块状、中厚层状，节理一般为2～5组，呈微张-闭合状，经统计表明，同一组节理间距一般为20～65 cm，岩体属较完整类.隧道围岩主体为Ⅲ～Ⅵ级，进出口段隧道围岩主体为Ⅴ～Ⅵ级.由地质情况分析得知拟模拟实例隧道的主体岩石考虑大多属于Ⅳ级围岩范畴，故实例隧道围岩等级取为Ⅳ级.

通过采用PLAXIS 3D TUNNEL(三维隧道软件)对实例隧道进行仿真计算分析，以衬砌轴应力σ的最大值σmax为研究对象，得到衬砌轴应力σmax关于随机变量c(粘聚力)、φ(摩擦角)、E(弹性模量)、D(衬砌厚度)的变化数据表，见表1、2、3、4所示.根据文献[8]的研究内容，在用有限元软件计算单因素变化条件下的衬砌轴应力值时，其他随机变量均取最低值[8]，例如：在粘聚力c变化时，摩擦角φ、弹性模量E、衬砌厚度D分别取31°、3 000 MPa、0.3 m.

表1 粘聚力变化下的衬砌轴应力

表2 摩擦角变化下的衬砌轴应力

表3 弹性模量变化下的衬砌轴应力

表4 衬砌厚度变化下的衬砌轴应力

2.2 隧道衬砌轴应力函数拟合

得到了上述4个表的对应数据后，采用数据分析软件Origin拟合成相应的函数，即单因素变化下的隧道衬砌轴应力的变化表达式.拟合结果如下图1～4所示.

通过数值分析拟合，得到了4条单因素变化下的衬砌轴应力函数.

σmax=-5.57+0.14×e2.35c

(1)

σmax=-5.31+0.002 2×e0.13φ

(2)

σmax=-1.85-6.65×e-0.000 232E

(3)

σmax=-1.82-14.37×e-4.83D

(4)

图1 单因素下衬砌轴应力随粘聚力的变化规律图Fig.1 Change diagram of lining axis stress withcohesive force under single factor

图2 单因素下衬砌轴应力随摩擦角的变化规律图Fig.2 Change diagram of lining axis stress withfriction angle under single factor

图3 单因素下衬砌轴应力随弹性模量的变化规律图Fig.3 Change diagram of lining axis stress withelastic modulus undersingle factor

图4 单因素下衬砌轴应力随衬砌厚度的变化规律图Fig.4 Change diagram of lining axis stress withlining thickness under single factor

为了对衬砌结构进行风险分析，在得到这些函数后，需要进一步拟合得到多因素下的应力变化函数.根据单因素拟合的特征规律，可以作如下假定，其多因素下的衬砌轴应力函数为以下形式：

σmax=A1+A2×e2.35c+A3×e0.13φ
+A4×e-0.000 232E+A5×e-4.83D

(5)

要求出待定系数A1,A2,A3,A4,A5，需要对其进行多因素回归分析.虽然从式(5)中我们可以发现衬砌轴应力σmax关于基本变量c、φ、E、D并不是线性关系，但我们可以令：c′=e2.35c、φ′=e0.13φ、E′=e-0.000 232E、D′=e-4.83D.通过上述的定义，衬砌轴应力σmax与各参数成为线性关系，可以进行线性多元回归分析.本文中多元线性回归以概率论与数理统计中的方法为基础，采用Matlab程序进行计算，得出最后的系数，则回归后的衬砌轴应力的表达式如式(6)所示.

σmax=1.126 3+0.136 4×e2.35c-
0.02e0.13φ-6.831×e-0.000 232E-
14.367 8×e-4.83D

(6)

2.3 隧道围岩功能函数的建立

本文最终要得到的结果是该实体隧道围岩的失效概率，需要根据其承载能力极限状态建立相应的功能函数.根据现行的隧道设计规范中的衬砌截面抗压和抗拉检验算式，来建立相应的衬砌截面抗拉强度和抗压强度的极限状态方程[9].

通过计算分析可知，本实例隧道分析的截面由抗压强度控制承载能力，即偏心矩eo≤0.2d(d为衬砌厚度)，相应的抗压强度极限方程如下所示：

za=N极限-N=KPRαbdRa-N=0

(7)

式中：N极限为衬砌混凝土所能承受的极限轴力(即抗力)；N为计算所得的截面轴力(即荷载效应)；KPR为抗力计算模式不确定性；b为纵向宽度，计算时取1 m；d为截面厚度；Ra为混凝土极限抗压强度；α为偏心影响系数.

实例隧道衬砌材料为钢筋混凝土材料，不考虑材料强度本身的随机变化，根据隧道设计资料，其混凝土的轴心抗压强度为12.5 MPa，故可按式(7)进行计算，由于σ为负值，因此该功能函数可表示为

za=KPRαbdRa-N=KPRαbdRa+σbd=0

(8)

公式(8)中同时略去截面厚度以及纵向宽度，即可得到式(9)的推导.

zα=12.5×KPR×α+σ=12.5×KPRα+

(1.126 3+0.136 4×e2.35c-0.02e0.13φ-

6.831×e-0.000 232E-14.367 8e-4.83D)

(9)

3 基于netica的隧道围岩失稳风险分析

3.1 隧道围岩安全状态简介

在得到了该隧道围岩的功能函数zα后，我们将结构功能函数的取值进行严格地划分，具体分为如下三种不同的状态(即可靠状态、极限状态、失效状态)[9]，见式(10).隧道围岩的失效概率即功能函数小于0的概率.

(10)

3.2 贝叶斯网络模型建立与分析

采用Netica软件建立隧道风险分析决策模型，并用该模型进行概率计算分析，得出最后隧道围岩的失效概率，进行相应的决策分析.出于对计算结果准确性的考虑，本隧道贝叶斯网络模型计算取样数设定为10 000 次[10]，贝叶斯网络计算模型见图5所示.

这时我们可以看到各个随机变量的数值概率分布情况，只要改变其中一个数值，该隧道围岩的失效概率将会发生改变.对于Ⅳ级围岩，衬砌厚度D的设计值一般大约为0.4～0.6 m左右[11]，通过贝叶斯网络得到衬砌厚度为0.4 m、0.5 m和0.6 m的情况下隧道围岩的失效概率.在贝叶斯网络中衬砌厚度选择0.4 m，结果如图6所示.由图可知在衬砌厚度为0.4 m时，该隧道围岩的失效概率f=2.77% ，根据相关规定，该失效概率满足实例隧道施工要求(本实例隧道安全可靠度控制级别为优时，失效概率≤5%).

图5 贝叶斯网络模型计算完成示意图Fig.5 Completion schematic of bayesian network model calculation

图6 衬砌厚度为0.4 m时围岩失效概率Fig.6 Failure probability of surrounding rock with lining thickness of 0.4 m

通过贝叶斯网络，可得到其他随机变量变化下的隧道围岩失效概率及其变化规律，如表5～8所示：

表5 不同衬砌厚度时得到的隧道围岩失效概率

由表5可知，Ⅳ级围岩时，当衬砌厚度越厚时，该隧道围岩的失效概率将会越低.因此，在进行隧道设计时，增加衬砌厚度，可降低隧道围岩的失效概率.但一味的增加衬砌厚度是不经济的，需通过计算对比分析经济性及安全性来确定最优的衬砌厚度.通过表6～8，可知隧道围岩的失效概率随着岩体粘聚力的增加而减小，随着摩擦角的增大而不断增大且变化幅度越来越大，随着弹性模量的增加而不断减小[12-13].

得到了相应的隧道围岩的失效概率随着各随机变量的变化规律后，可以通过不同的条件选择不同的施工对策来减小隧道开挖施工的风险，为隧道工程的安全施工提供一定的参考性[14-15].

表6 不同粘聚力时得到的隧道围岩失效概率

表7 不同摩擦角时得到的隧道围岩失效概率

表8 不同弹性模量时得到的隧道围岩失效概率

4 结论

本文提出了贝叶斯网络的概念，确立了本文基于贝叶斯网络的隧道围岩风险分析决策的理念.通过采用有限元分析软件得到隧道风险分析参数，并选用贝叶斯网络工具Netica软件来对本文隧道风险分析结构进行贝叶斯网络的建模.通过理论计算和实例分析得到以下结论：

(1)采用有限元软件建模可以得到隧道围岩衬砌轴应力关于隧道各个随机变量的数据变化曲线，并可拟合得到多因素下的应力变化函数，再根据其承载能力极限状态建立相应的功能函数，为隧道风险分析模型的建立提供数据支持.

(2)当衬砌厚度越厚时，该隧道围岩的失效概率将会越低.因此，在进行隧道设计时，增加衬砌厚度，可降低隧道围岩的失效概率.但一味的增加衬砌厚度是不经济的，需通过计算对比分析经济性及安全性来确定最优的衬砌厚度.

(3)隧道围岩的失效概率随着岩体粘聚力的增加而减小，随着摩擦角的增大而不断增大且变化幅度越来越大，随着弹性模量的增加而不断减小.

(4)得到了相应的隧道围岩的失效概率随着各随机变量的变化规律后，可以通过不同的条件选择不同的施工对策来减小隧道开挖施工的风险，为隧道的安全施工提供一定的参考性.