隧道可靠度BFOA逆向寻优方法及其工程应用

2022-05-07郑帅姜谙男董庆波罗国成

广西大学学报（自然科学版） 2022年1期

郑帅，姜谙男*，董庆波，罗国成

(1.大连海事大学交通运输工程学院，辽宁大连 1160262.中铁建大桥工程局集团第一工程有限公司，辽宁大连 116000)

0 引言

受地质演化等自然条件影响，隧道围岩环境大多表现为高度的不确定性，可靠度评价方法能够有效地实现对这种不确定性的定量化评价，即通过失稳概率等指标描述隧道稳定状态。孙明社等[1]将二次响应面函数与非概率可靠度分析方法相结合，把可靠度指标的解答转化为基于功能函数梯度信息的最优化问题求解。陈钒等[2]推导得出了能够考虑围岩膨胀性和腐蚀性的隧道初砌结构可靠度指标计算公式，并建立了不同模式下的时变可靠度分析模型。尹蓉蓉[3]建立了考虑衬砌结构抗力模糊性及荷载效应随机性的正常使用极限状态下模糊随机可靠度分析方法，并对某海底隧道海域Ⅲ级(Ⅳ类)围岩下衬砌结构最大裂缝宽度的可靠度进行了计算。这些方法均是在已知工程支护参数的前提下对结构失稳概率评价方法进行研究，这种研究方法为正向可靠度方法(forward reliability method)[4]，然而在实际工程中，往往需要根据工程特征在设计阶段定义目标工程的预期可靠度指标，并根据该指标进行相应的支护方案设计，近年来一些相关工程可靠度设计标准[5-7]的施行逐渐规范了这一过程，即逆向可靠度方法(inverse reliability method)[8]。李翔等[9]引入了逆向可靠度方法并归纳其具体计算实施流程，实现了在已知目标失稳概率条件下，反求工程参数的目的。方砚兵等[10]从风险水平的角度初步量化分级安全储备，提出在给定可靠性水平下直接计算隧道支护抗力的方法。这些研究实现了隧道可靠度的逆向计算过程，其解析式可靠度逆向求解方法收敛性较好；但对岩体本构方程存在一定的依赖性且各随机参数被认为满足于同一种分布特征，限制了该方法的应用效果，因此有必要针对这些问题探索更为合适的求解方法。

在已有研究结果的基础上，本文建立一种基于软弱围岩失稳变形极限状态方程与细菌觅食优化算法(BFOA)的隧道可靠度逆向求解方法。首先采用概率图法对现场所采集的围岩参数样本进行分布特征统计；然后基于围岩变形极限状态方程建立隧道建立映射关系，通过Monte Carlo抽样的方式实现可靠度求解并以此作为优化算法搜索过程中的响应面函数；最后通过BFOA寻优的方式根据工程的预期可靠度进行支护阻力的逆向求解。在此基础上，讨论了迭代次数、搜索区间开度对最终计算结果精度及计算效率的影响，同时分析各围岩参数相对于计算结果的敏感性以备工程应用所参考。

1 隧道可靠度逆向计算方法

1.1 围岩变形的基本力学模型

对围岩力学行为的合理描述是进行隧道可靠度评价的基础，在此采用文献[11-13]提出的针对软弱介质岩石隧道进行弹塑性力学分析工作所建立的围岩与支护相互作用关系模型，构建可靠度计算所需的功能函数，详细推导过程参照文献[9]，在此仅给出最终推导结果及部分关键参数表达方程。围岩径向塑性位移：

(1)

pcr=2(p0-σc)/(k+1)，k=(1+sinφ)/(1-sinφ)，

(2)

(3)

式中：uip为围岩径向塑性位移；r0为隧道开挖半径；rp为洞周塑性区半径；p0为围岩所受地应力；pcr为临界应力；ps为支护阻力；E为围岩弹性模量；μ为泊松比；c为黏聚力;φ为内摩擦角。

联立式(1)—(3)，则在工程已知的情况下，可通过实验、勘察等方法获取围岩基本力学参数与隧道埋深、初始支护阻力等参数取值结果，进而求解方程获得塑性区半径rp及径向位移uip。

进一步地，根据rp与uip计算结果对围岩稳定状态进行评价，采用双评价指标，即

umax≥uip，rp/r0≤L，

(4)

式中：umax为目标工程设计结果中洞周围岩位移限值；L为与rp/ro有关的容许阈值[14-15]。当uip与rp计算结果同时满足式(4)时认为隧道围岩处于稳定状态，否则为失稳状态。

1.2 可靠度方法

考虑隧道力学参数随机性特点，引入可靠度理论，将最终的隧道稳定性计算结果的可靠或失效概率看为一个两类别的分类问题[16]，即隧道稳定性的概率结果，这就是隧道可靠度问题。

(5)

为了保证足够的抽样次数以满足计算精度，采用106次Monte Carlo方法进行隧道可靠度计算，具体计算过程可参照文献[17]。其中，通过前文式(4)构建可靠度功能函数

(6)

采用双评价指标，即当Z1(x)、Z2(x)同时不小于0时为隧道安全状态，其可靠度概率为

(7)

1.3 BFOA

BFOA是一个最小值求解过程，主要实现流程为：首先在目标区间内随机生成细菌种群，每一个细菌个体均代表了一个待求支护阻力Ps值。依据算法策略进行各细菌的趋向、复制、迁徙操作，这个过程模仿了人类大肠杆菌觅食特征，详细过程可参照文献[19]。各项操作上限次数分别限制为Nc、Nre、Ned，以搜索到的目标值与预期设计值误差最小为计算目的，BFOA实现流程如图1所示。

1.4 可靠度逆向计算流程

综上1.1—1.3节，实现了可靠度逆向求解过程，计算流程如图2所示，具体实现步骤为：

①依托背景工程，通过现场取样、室内试验等方式获取围岩基本力学参数样本组，参数包括围岩弹性模量、泊松比、黏聚力、内摩擦角。

② 根据得到的参数样本组，采用概率图法分别验证每个参数所属分布种类，分布验证目标种类包括正态分布、负指数分布、泊松分布、均匀分布、对数分布，统计获得各参数的特征函数。

③以初始地应力p0、支护阻力ps为不变量，随机抽取围岩力学参数E、μ、c、φ组成输入样本，通过式(1)计算获得围岩塑性区半径rp及径向位移uip，进一步根据式(4)判断当前稳定状态。大量重复这个过程(本研究设置为106次)，统计稳定性评价结果获得可靠度指标。

④ 根据现场实际情况确定待求支护阻力区间，随机生成初始种群后进行1.3节中BFOA寻优过程，以上一步过程作为响应面函数，直至满足初始设定的终止条件(最小适应值或最大迭代步)。

⑤ 记录并输出与预期可靠度相匹配的支护阻力值，完成可靠度逆向计算过程。

以上，实现了基于BFOA寻优的隧道可靠度逆向寻优计算流程。

图1 BFOA实现流程Fig.1 BFOA implementation process

图2 逆向可靠度求解流程图Fig.2 Reverse reliability solution process

2 工程应用

2.1 工程背景

图3 掌子面围岩破碎状态Fig.3 Broken rock of tunnel face

某隧道为高速公路隧道，受区域构造影响，区内存在多条异向断裂构造带，主要表现为压性断裂，以逆断层和挤压破碎带为主；围岩介质多为花岗质片麻岩，岩体完整性较差，见有玄武岩分布。受地质条件影响，隧道施工过程中面临的围岩状态多如图3所示，节理裂隙体系发育程度较为剧烈，岩石整体性、隧道自稳性较差。

由于地质环境变化剧烈，传统的依靠施工经验确定支护方案的模式显然在此不宜适用,因无法准确判断围岩环境的真实状态，经验判断极易造成过度支护或支护不足状态进而导致施工成本浪费或失稳破坏两种极端状态。根据相关规范[5-7]，预留适当的安全余量，在保证工程可靠度的前提下进行支护方案设计是必要的，即通过可靠度的逆向求解过程求解所需支护阻力。

该隧道设计为左、右两幅，单幅隧道净宽10.69 m，净高7.10 m，两幅间距32 ～40 m，隧道间距超过3倍洞径，因此在研究过程中忽略两洞间相互影响，在此以右洞为例进行研究。

2.2 参数分布统计

2.2.1 样本采集

图4 围岩参数试验获取Fig.4 Test acquisition of surrounding rock parameterssurrounding rock parameters distribution

研究区域里程桩号为K95+080～K95+180，在区域内均匀采样获取了40个岩石样本，通过实验室三轴压缩试验测试获取岩样基本力学参数E、μ、c、φ：以1 MPa/min的速率将围压加载至5 MPa(根据围岩地应力情况确定)，保持围压不变，以1 mm/min的移动速率加载轴压，直至试样破坏获得岩样三轴瞬时强度，围岩参数试验获取如图4所示。根据实验结果，围岩力学参数统计结果见表1。

应当说明，三轴压缩至试样破裂试验过程所需时间较短，因此该样本采集过程是可以在施工过程中快速实现的。条件允许时亦可将项目管理方的岩体样本试验结果作为补充数据完善该采集过程。

表1 围岩力学参数统计结果Tab.1 Statistical of surrounding rock parameters

2.2.2 分布规律的概率图法验证

为了验证某个参数的统计样本x1,x2,…,xn是否满足某种分布特征函数F(x)，采用概率图法进行计算：首先对样本序列由小到大进行排列，得到递升序列y1,y2,…,yn，与之对应的分布函数则为F(yi)=i/n。

定义经验分位数取值为

(8)

依据所选定的分布特征函数F(x)，计算与每个经验分位数所匹配的特征分布值：

θi=F-1(qi) 。

(9)

以yi为横坐标，θi为纵坐标在坐标系中描绘出n个数据点，并拟合统计验证曲线。理论上，若统计样本x1,x2,…,xn满足分布特征函数F(x)，则该n个点将呈现为一条与坐标轴夹角为45°的直线，实际应用过程中，考虑测量取值误差等因素，仅需该n个点沿直线分布即可认为该组数据满足所选用的概型，而对截距和斜率不再做要求。

对表1中围岩4种基本力学参数采用概率图法进行分布特征验证，结果表明，围岩弹性模量服从正态分布、泊松比服从负指数分布、黏聚力服从对数正态分布、摩擦角服从正态分布。

2.2.3 分布特征统计及抽样函数

根据式(2)中确定的各参数分布特征函数，通过式(1)所述数据进行特征参数拟合，得到参数分布特征验证结果见表2，基于MATLAB平台实现Monte Carlo计算过程中相应的抽样函数见表2中最后一列。

表2 参数分布特征验证结果Tab.2 Parameter distribution feature verification result

2.3 结果分析

将研究区域等距离划分为4个区间，根据区间内平均埋深通过P0=γ·h计算隧道围压，其中γ为岩体重度，取值参照地质勘查结果；在单一区间内，支护阻力Ps1水平取为0.15、0.25、0.35 MPa，共形成12组试验组并采用前文1.2节正向可靠度方法计算各组可靠度指标βfor；令逆向可靠度计算的预期可靠度指标βtar=βfor，采用1.4节逆向可靠度方法计算待求支护阻力Ps2，计算精度设置为1×10-3，搜索终止条件为满足最小计算精度或最大迭代次数，可靠度逆向计算结果验证见表3。

表3 可靠度逆向计算结果验证Tab.3 Verification of inverse reliability calculation results

由表3计算结果可见，根据预期可靠度指标βtar，通过逆向可靠度求解所得的待求支护阻力Ps2与正向可靠度计算过程中支护阻力Ps1一致性较好。其中最大绝对误差出现在第3试验组，为9.2×10-4，此时相对误差率为0.712‰，说明所建立的可靠度逆向求解方法能够实现待求设计参数的准确求解。

3 讨论

计算过程中发现，逆向可靠度计算结果精度与计算效率受搜索区间开度、迭代次数与施工环境等多条件影响，在此对各因素影响规律进行总结分析，为相应的工程应用提供参考。

本节所述计算效率问题是指在常规工作电脑上的运行结果，参考配置为：i5-8250U型处理器，NVIDIA-940MX型显卡，8GB-DDR4型运行内存。

3.1 搜索区间开度影响

对于表3所示验证过程，由于在正向可靠度计算过程中已知待求阻力值，逆向计算时将搜索区间设置为Ps±0.5，即区间开度为1。在实际工程中应用逆向可靠度方法根据预期可靠度求解支护阻力时，待求目标参数无法细化认定准确搜索区间，需要适当扩大搜索范围，因此有必要探讨搜索区间开度对计算过程的影响。

以表3中第1、4、7三组为例，分别将搜索区间以0.5为单位向两侧扩展，在区间开度为1～5的条件下进行逆向可靠度计算，搜索终止条件为结果精度达到1×10-3，记录不同区间开度条件下、不同试验组的计算时间与所需BFOA迭代次数，不同搜索区间效率对比如图5所示。由图5可见，扩大搜索区间会对计算效率形成一定影响，区间开度分别为1和5时求解时间最大差值为252 s，显然这在实际应用过程中是可以接受的；而随着搜索区间的扩大，最终逆向求解结果精度并未受到明显影响。因此，在工程建设中应用逆向可靠度方法求解支护阻力等参数时可适当放大搜索区间保证其能够覆盖最优解所在位置。

图5 不同搜索区间效率对比Fig.5 Calculation efficiency of different search intervals

3.2 迭代次数影响

这里的迭代次数分为两方面：一为BFOA寻优过程所需进行的正向目标函数寻优次数k1；二为单次正向计算时Monte Carlo抽样次数k2。

图6 寻优迭代次数对计算精度影响曲线Fig.6 Calculation accuracy curve of optimization times

对于1.3节所述BFOA寻优算法，其完成计算所需调用1.2节正向可靠度计算模型的迭代次数k1为g次趋化操作、n次复制操作、m次迁徙操之乘积：k1=g×n×m。前文采用k1=50×8×4，k2=1×106，在此基础上以表3中误差最大的第3试验组为例进行不同迭代次数对逆向可靠度计算效率与精度影响规律探索，搜索终止条件设置为满足预设迭代次数。寻优迭代次数对计算精度影响曲线如图6所示，抽样迭代次数对计算精度影响曲线如图7所示。

由图6和图7可见，计算精度会随着迭代次数的增加快速增大，当迭代次数达到某一定值后趋于稳定，即迭代次数的设置存在阈值，当超过阈值后计算精度才能得到有效保证。综合计算结果，选取迭代次数k1=1 600，抽样次数k2=1×105，此时能够在保证足够计算精度的前提下获得最高的计算效率。

3.3 围岩参数敏感性分析

实际应用中发现，不同围岩环境下逆向可靠度求解所得支护阻力波动剧烈，说明围岩参数对支护阻力的设计存在影响。在此以逆向可靠度的形式对围岩力学参数敏感性进行分析，基于前述逆向可靠度计算过程，待分析参数不再参与抽样，而是在表2所示取值区间内等距选取10个试算点，将待分析参数依次取为这10个定值进行支护阻力值的逆向求解，以此类推，完成全部参数的试算后得到围岩参数敏感性如图8所示。