工程科学学报，第41卷，第1期：78-87,2019年1月 Chinese Journal of Engineering,Vol.41,No.I:78-87,January 2019 D0L:10.13374/j.issn2095-9389.2019.01.008;htp:/journals.ustb.edu.cm 深埋硬岩隧道围岩参数概率反演方法 吴忠广12)四，吴顺川13) 1)北京科技大学土木与资源工程学院，北京1000832)交通运输部科学研究院标准与计量研究中心，北京100029 3)昆明理工大学国土资源学院，昆明650093 ☒通信作者，E-mail:kinliwu@163.com 摘要在贝叶斯理论框架下，提出了一种基于多源数据融合的深埋硬岩隧道围岩参数概率反演方法.首先，分析硬岩隧道 常用的启裂一剥落界限本构模型中围岩单轴抗压强度、启裂强度与抗压强度比及抗拉强度三个参数不确定性来源，确定其概 率统计特征；其次，利用粒子群算法优化多输出支持向量机，建立反映反演参数与隧道监测数据间非线性映射关系的智能响 应面：最后，结合贝叶斯分析方法构建概率反演模型，运用马尔科夫链蒙特卡洛模拟算法实现了围岩参数的动态更新.将该方 法应用到某深埋硬岩隧道中，利用反演的围岩参数计算隧道拱顶下沉点、周边收敛点变化值及开挖损伤区深度，与监测数据 吻合较好.结果表明，该方法可以实现围岩多参数快速概率反演，更新后的参数可用于硬岩隧道施工安全风险评估与结构可 靠性设计. 关键词硬岩隧道：概率反演：多源数据融合：贝叶斯理论：多输出支持向量机 分类号TU452 Probabilistic back analysis method for determining surrounding rock parameters of deep hard rock tunnel WU Zhong-guang.2)回，WU Shun-chuan',3) 1)School of Civil and Resource Engineering,University of Science and Technology Beijing,Beijing 100083.China 2)Research Center for Standards and Metrology,China Academy of Transportation Sciences,Beijing 100029,China 3)Faculty of Land Resource Engineering,Kunming University of Science and Technology,Kunming 650093.China Corresponding author,E-mail:kinliwu@163.com ABSTRACT A large number of tunnel projects are being constructed or will be constructed in the mountainous areas of western Chi- na.However,they are several safety challenges in the construction of deep hard rock tunnels because of the complex topographic and geological conditions,strong geological tectonic activities,large burial depth,and high in situ stress level.Uncertainty of tunnel wall parameters is one of main factors that contribute to tunnel construction risk.The traditional deterministic back analysis method cannot reflect the uncertainty characteristics of tunnel wall parameters;therefore,within the framework of Bayesian theory,a probabilistic back analysis method based on integrating multi-source monitoring information was proposed for determining the surrounding rock parameters of deep hard rock tunnel.First,the uncertainty sources of three parameters-uniaxial compressive strength(UCS),crack initiation stress to UCS ratio,and tensile strength for the widely used damage initiation and spalling limit approach-were analyzed,and their probabilistic statistical characteristics were determined.Second,a multi-output support vector machine(MSVM)was optimized by par- ticle swarm optimization (PSO)algorithm,and an intelligent response surface model was established to reflect the nonlinear mapping relationship between back-analyzed parameters and field monitoring data.Last,by combination with the Bayesian (B)analysis meth- od,the B-PSO-MSVM model was established,and surrounding rock parameters were dynamically updated by applying the Markov 收稿日期：2018-05-29 基金项目：国家重点研发计划专项资助项目(2017YF008053003):国家自然科学基金资助项目(51174020)
工程科学学报,第 41 卷,第 1 期:78鄄鄄87,2019 年 1 月 Chinese Journal of Engineering, Vol. 41, No. 1: 78鄄鄄87, January 2019 DOI: 10. 13374 / j. issn2095鄄鄄9389. 2019. 01. 008; http: / / journals. ustb. edu. cn 深埋硬岩隧道围岩参数概率反演方法 吴忠广1,2) 苣 , 吴顺川1,3) 1) 北京科技大学土木与资源工程学院, 北京 100083 2) 交通运输部科学研究院标准与计量研究中心, 北京 100029 3) 昆明理工大学国土资源学院, 昆明 650093 苣 通信作者, E鄄mail:kinliwu@ 163. com 摘 要 在贝叶斯理论框架下,提出了一种基于多源数据融合的深埋硬岩隧道围岩参数概率反演方法. 首先,分析硬岩隧道 常用的启裂—剥落界限本构模型中围岩单轴抗压强度、启裂强度与抗压强度比及抗拉强度三个参数不确定性来源,确定其概 率统计特征;其次,利用粒子群算法优化多输出支持向量机,建立反映反演参数与隧道监测数据间非线性映射关系的智能响 应面;最后,结合贝叶斯分析方法构建概率反演模型,运用马尔科夫链蒙特卡洛模拟算法实现了围岩参数的动态更新. 将该方 法应用到某深埋硬岩隧道中,利用反演的围岩参数计算隧道拱顶下沉点、周边收敛点变化值及开挖损伤区深度,与监测数据 吻合较好. 结果表明,该方法可以实现围岩多参数快速概率反演,更新后的参数可用于硬岩隧道施工安全风险评估与结构可 靠性设计. 关键词 硬岩隧道; 概率反演; 多源数据融合; 贝叶斯理论; 多输出支持向量机 分类号 TU452 收稿日期: 2018鄄鄄05鄄鄄29 基金项目: 国家重点研发计划专项资助项目(2017YFC08053003);国家自然科学基金资助项目(51174020) Probabilistic back analysis method for determining surrounding rock parameters of deep hard rock tunnel WU Zhong鄄guang 1,2) 苣 , WU Shun鄄chuan 1,3) 1) School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China 2) Research Center for Standards and Metrology, China Academy of Transportation Sciences, Beijing 100029, China 3) Faculty of Land Resource Engineering, Kunming University of Science and Technology, Kunming 650093, China 苣 Corresponding author, E鄄mail: kinliwu@ 163. com ABSTRACT A large number of tunnel projects are being constructed or will be constructed in the mountainous areas of western Chi鄄 na. However, they are several safety challenges in the construction of deep hard rock tunnels because of the complex topographic and geological conditions, strong geological tectonic activities, large burial depth, and high in situ stress level. Uncertainty of tunnel wall parameters is one of main factors that contribute to tunnel construction risk. The traditional deterministic back analysis method cannot reflect the uncertainty characteristics of tunnel wall parameters; therefore, within the framework of Bayesian theory, a probabilistic back analysis method based on integrating multi鄄source monitoring information was proposed for determining the surrounding rock parameters of deep hard rock tunnel. First, the uncertainty sources of three parameters———uniaxial compressive strength (UCS), crack initiation stress to UCS ratio, and tensile strength for the widely used damage initiation and spalling limit approach———were analyzed, and their probabilistic statistical characteristics were determined. Second, a multi鄄output support vector machine (MSVM) was optimized by par鄄 ticle swarm optimization (PSO) algorithm, and an intelligent response surface model was established to reflect the nonlinear mapping relationship between back鄄analyzed parameters and field monitoring data. Last, by combination with the Bayesian (B) analysis meth鄄 od, the B鄄PSO鄄MSVM model was established, and surrounding rock parameters were dynamically updated by applying the Markov
吴忠广等：深埋硬岩隧道围岩参数概率反演方法 .79 Chain Monte Carlo simulation algorithm.The method was applied to a deep hard rock tunnel,and parameters from probabilistic back analysis were utilized to calculate the point change of the tunnel vault settlement and peripheral displacement convergence as well as the depth of excavation damage zones,and the results agreed well with the actual monitoring data.It is shown that this method can be used to back analyze multi parameters of surrounding rock quickly and probabilistically,and parameters updated can be applied for risk as- sessment in construction safety and structural reliability design for the hard rock tunnel. KEY WORDS hard rock tunnel;probabilistic back analysis;multi-source information fusion;Bayesian theory;multi-output support vector machine 岩土工程的不确定性导致风险已成为业界共 系，根据监测数据的变化反演计算围岩参数，实现参 识，不确定性因素包括岩体性质、荷载、几何尺寸、初 数的动态更新，响应面方法常近似作为表征这种复 始条件、边界条件、计算模型、破坏机理等，Baecher 杂映射关系的替代模型】.传统的多项式响应面 与Christian)、Der Kiureghian与Ditlevsen、Lang- 方法随着阶数的增加所需的计算量增大[24]，对于复 ford与Diederichst3]等根据数据信息的掌握程度将 杂的非线性地下工程问题，多项式响应面拟合效果 不确定性分为随机不确定性和认知不确定性两类， 不佳]，甚至可能会产生伪极限状态面[2].基于神 这两类不确定性分类方法在岩土工程领域逐渐得到 经网络或者支持向量机的智能响应面具有利用较少 广泛认可.深埋硬岩隧道围岩力学参数具有以上两 的样本数据进行高阶逼近的优点，逐渐被应用于岩 种不确定性特征，围岩力学性质与岩体中的结构面、 石工程智能分析领域中[2”-2].神经网络方法采用 结构体及赋存环境密切相关[，常表现出非线性、 经验风险最小化准则，存在收敛速度慢、泛化能力较 各向异性、尺寸效应等特征，导致围岩力学参数存在 差、过学习、局部极小点等缺陷[]，而支持向量机基 随机不确定性[].然而，由于现场数据缺失、测量不 于结构风险最小化准则，能较好的解决小样本、高维 准确、计算模型误差等原因，造成围岩力学参数也存 数、非线性、局部极小点等实际问题03).经典的 在认知不确定性).因此，合理确定围岩力学参数， 支持向量机方法主要是为解决二分类问题提出的， 科学量化围岩参数不确定性对于研判深埋硬岩隧道 是单输出模型，对于求解多输出的岩土工程问题，通 工程结构安全性能及施工安全风险具有重要的 常采用建立各测点独立的支持向量机方法，不但会 意义 增加计算量，而且忽略各测点数据间相关性的做法 概率反演是不确定性分析的有效方法，它考虑 会带来反演误差，多输出支持向量机(multi-output 围岩参数先验信息，使反演更能反映实际的不确定 support vector machine,MSVM)有效解决了以上不 性情况.概率反演方法在岩土工程中应用较为广 足[2-4].但是，在应用多输出支持向量机方法时， 泛，解决了很多工程实际问题，如Gilbert等s-6)、 要首先解决模型的参数优化问题，参数选择是否合 Zhang等-]、Zhang等1、Wang等o,i等利用 适对于模型的学习和推广能力影响很大.采用人工 概率反演方法分析了边坡工程不确定性问题，M- 搜索最优参数的方法费时耗力，粒子群(particle randa等[2]、Miro等[]通过对岩土参数进行概率反 swarm optimization,PS0)优化算法[3s]具有高效的全 演，研究其不确定性对隧道结构及周围环境的影响. 局寻优能力，对解决不同类型的支持向量机参数优 贝叶斯方法能够融合岩土工程位移、应力等多源监 化问题应用效果较好36] 测数据，利用概率反演计算参数的后验分布，合理表 本文结合贝叶斯方法(Bayesian)与粒子群优化 征参数的不确定性.目前，贝叶斯概率反演已成为 的多输出支持向量机模型，提出一种新的深埋硬岩 岩土力学参数不确定性问题研究的热点.Haas与 隧道围岩参数动态概率反演贝叶斯-粒子群-多输 Einstein14)、Miranda等[is]、Zhang等[6]、Peng等】 出支持向量机（以下简称B-PSO-MSVM)方法.将该 均利用马尔科夫链蒙特卡洛方法开展了贝叶斯分 方法用于反演硬脆性岩体常用的启裂-剥落界限 析，Feng与Jimenez[is]、Wang等[u9-o】、Wang与 (damage initiation and spalling limit,DISL)本构模型 Aladejare2]等提出了等效岩土力学参数与分析参 三个基本参数，即围岩单轴抗压强度(uniaxial com- 数间相关系数的贝叶斯计算方法，Contreras等22]利 pressive strength,UCS)、抗拉强度(tensile strength, 用贝叶斯方法对Hoek-Brown强度准则参数开展了 T)以及启裂强度与抗压强度比(the ratio of crack in- 反演研究与不确定性分析.贝叶斯概率反演方法通 itiation and UCS,CL/UCS),并量化其不确定性，为 过建立围岩参数与工程监测数据间的非线性映射关 分析隧道围岩脆性力学行为及结构安全性能提供依
吴忠广等: 深埋硬岩隧道围岩参数概率反演方法 Chain Monte Carlo simulation algorithm. The method was applied to a deep hard rock tunnel, and parameters from probabilistic back analysis were utilized to calculate the point change of the tunnel vault settlement and peripheral displacement convergence as well as the depth of excavation damage zones, and the results agreed well with the actual monitoring data. It is shown that this method can be used to back analyze multi parameters of surrounding rock quickly and probabilistically, and parameters updated can be applied for risk as鄄 sessment in construction safety and structural reliability design for the hard rock tunnel. KEY WORDS hard rock tunnel; probabilistic back analysis; multi鄄source information fusion; Bayesian theory; multi鄄output support vector machine 岩土工程的不确定性导致风险已成为业界共 识,不确定性因素包括岩体性质、荷载、几何尺寸、初 始条件、边界条件、计算模型、破坏机理等,Baecher 与 Christian [1] 、Der Kiureghian 与 Ditlevsen [2] 、Lang鄄 ford 与 Diederichs [3] 等根据数据信息的掌握程度将 不确定性分为随机不确定性和认知不确定性两类, 这两类不确定性分类方法在岩土工程领域逐渐得到 广泛认可. 深埋硬岩隧道围岩力学参数具有以上两 种不确定性特征,围岩力学性质与岩体中的结构面、 结构体及赋存环境密切相关[4] ,常表现出非线性、 各向异性、尺寸效应等特征,导致围岩力学参数存在 随机不确定性[2] . 然而,由于现场数据缺失、测量不 准确、计算模型误差等原因,造成围岩力学参数也存 在认知不确定性[2] . 因此,合理确定围岩力学参数, 科学量化围岩参数不确定性对于研判深埋硬岩隧道 工程结构安全性能及施工安全风险具有重要的 意义. 概率反演是不确定性分析的有效方法,它考虑 围岩参数先验信息,使反演更能反映实际的不确定 性情况. 概率反演方法在岩土工程中应用较为广 泛,解决了很多工程实际问题,如 Gilbert 等[5鄄鄄6] 、 Zhang 等[7鄄鄄8] 、Zhang 等[9] 、Wang 等[10] 、Li 等[11] 利用 概率反演方法分析了边坡工程不确定性问题,Mi鄄 randa 等[12] 、Miro 等[13]通过对岩土参数进行概率反 演,研究其不确定性对隧道结构及周围环境的影响. 贝叶斯方法能够融合岩土工程位移、应力等多源监 测数据,利用概率反演计算参数的后验分布,合理表 征参数的不确定性. 目前,贝叶斯概率反演已成为 岩土力学参数不确定性问题研究的热点. Haas 与 Einstein [14] 、Miranda 等[15] 、Zhang 等[16] 、Peng 等[17] 均利用马尔科夫链蒙特卡洛方法开展了贝叶斯分 析, Feng 与 Jimenez [18] 、 Wang 等[19鄄鄄20] 、 Wang 与 Aladejare [21] 等提出了等效岩土力学参数与分析参 数间相关系数的贝叶斯计算方法,Contreras 等[22] 利 用贝叶斯方法对 Hoek鄄Brown 强度准则参数开展了 反演研究与不确定性分析. 贝叶斯概率反演方法通 过建立围岩参数与工程监测数据间的非线性映射关 系,根据监测数据的变化反演计算围岩参数,实现参 数的动态更新,响应面方法常近似作为表征这种复 杂映射关系的替代模型[23] . 传统的多项式响应面 方法随着阶数的增加所需的计算量增大[24] ,对于复 杂的非线性地下工程问题,多项式响应面拟合效果 不佳[25] ,甚至可能会产生伪极限状态面[26] . 基于神 经网络或者支持向量机的智能响应面具有利用较少 的样本数据进行高阶逼近的优点,逐渐被应用于岩 石工程智能分析领域中[27鄄鄄29] . 神经网络方法采用 经验风险最小化准则,存在收敛速度慢、泛化能力较 差、过学习、局部极小点等缺陷[25] ,而支持向量机基 于结构风险最小化准则,能较好的解决小样本、高维 数、非线性、局部极小点等实际问题[30鄄鄄31] . 经典的 支持向量机方法主要是为解决二分类问题提出的, 是单输出模型,对于求解多输出的岩土工程问题,通 常采用建立各测点独立的支持向量机方法,不但会 增加计算量,而且忽略各测点数据间相关性的做法 会带来反演误差,多输出支持向量机( multi鄄output support vector machine,MSVM) 有效解决了以上不 足[32鄄鄄34] . 但是,在应用多输出支持向量机方法时, 要首先解决模型的参数优化问题,参数选择是否合 适对于模型的学习和推广能力影响很大. 采用人工 搜索最优参数的方法费时耗力,粒子群( particle swarm optimization,PSO)优化算法[35] 具有高效的全 局寻优能力,对解决不同类型的支持向量机参数优 化问题应用效果较好[36] . 本文结合贝叶斯方法(Bayesian)与粒子群优化 的多输出支持向量机模型,提出一种新的深埋硬岩 隧道围岩参数动态概率反演贝叶斯鄄鄄 粒子群鄄鄄 多输 出支持向量机(以下简称 B鄄PSO鄄MSVM)方法. 将该 方法用于反演硬脆性岩体常用的启裂鄄鄄 剥落界限 (damage initiation and spalling limit, DISL)本构模型 三个基本参数,即围岩单轴抗压强度( uniaxial com鄄 pressive strength, UCS)、抗拉强度( tensile strength, T)以及启裂强度与抗压强度比(the ratio of crack in鄄 itiation and UCS, CI/ UCS),并量化其不确定性,为 分析隧道围岩脆性力学行为及结构安全性能提供依 ·79·
·80. 工程科学学报，第41卷，第1期 据.首先，分析启裂-剥落界限模型参数的不确定性 Perras与Diederichs[]等进一步提出了以岩石启裂 特征，进而提出围岩参数动态概率反演方法，最后， 强度CI为基础的抗拉强度T经验确定方法.因此， 利用工程案例证明该方法的适用性与可行性 由于实验误差或者模型不确定性等导致完整岩石的 抗拉强度T取值存在不确定性 1启裂-剥落界限模型的不确定性 2概率反演B-PSO-MSVM方法 1.1启裂-剥落界限模型 启裂-剥落界限模型是Diederichs]在H-B强 2.1贝叶斯概率反演方法 度准则3]的基础上提出的，常用于隧道围岩脆性破 假设隧道工程计算模型可用函数[)表示为： 坏深度预测[90]与层裂破坏风险分析].HB强 =g(9) (2) 度准则可由公式(1)表示： 式中，0为模型的计算参数：y为模型响应的计算值， 。=g+Ucsa晨+广 (1) 如位移、破裂深度等，则实际监测值y为： y=y+s (3) 式中，σ和σ分别为破坏时的最大和最小有效主应 式中，ε为模型残差，假设ε服从均值为0，方差为 力；UCS为完整岩石的单轴抗压强度；m为完整岩 σ2的正态分布，即e~N(0,o2) 石的材料常数m,的折减值，m的取值范围为 设日为随机变量，此处，0为启裂-剥落界限模 0.0000001~25:s为岩石的完整性系数，s的取值范 型的三个基本参数{UCS、C/UCS、T}的集合，可通 围为0~1：a为包络线曲率参数. 过实验或历史数据统计获得先验概率密度函数 DISL模型利用UCS、CI及T三个参数确定岩体 强度曲线的峰值参数与残余参数[7]，具体取值见下 f代).由贝叶斯公式可知，0的后验概率密度函数 f(ly)为： 表1. f(0ly)=Af(yle)f(0) (4) 表1DISL模型参数取值表 式中：入为归一化系数：f(ly)为似然函数 Table 1 DISL model input parameters 对于多源监测数据，如已知y=(P。,P,EDZ), 参数 峰值参数 残余参数 其中，P,和P,分别为拱顶下沉点与周边收敛点变化 0.25 0.75 值，EDZ为开挖损伤区深度(excavation damage CI/UCS) 0.001 zone,EDZ),假设监测信息间彼此独立，则 s(UCS/ITI) 612 f(0ly)=f(Po10)f(P,10)f(EDZI0)(5) 1.2模型参数不确定性来源 由于式(4)中后验概率分布是复杂、高维、非标 由表1可知，使用启裂-剥落界限模型需要事 准的形式，直接积分计算十分困难，马尔科夫链蒙特 先确定UCS、CI及T三个参数值，这些参数具有概 卡洛(markov chain monte carlo,MCMC)模拟方 率统计特征.原因如下：(1)UCS是由完整岩石的 法[6-是处理复杂统计学问题的有效工具，其核心 室内单轴抗压强度实验得到，由于岩石中物质成分 思想是通过构造一条马尔科夫链，然后按照转移核 的不同与微结构面的存在导致岩石具有各向异性， 规则引导马尔科夫链扰动过程使其逼近目标分布， 岩样制作及实验测量误差等使得岩石单轴抗压强度 抽取逼近后的样本来近似计算后验分布.MCMC的 具有很大的变异性：(2)启裂强度C1主要是由实验 算法很多，Metropolis-Hasting[-9]算法应用广泛，主 室声发射试验来确定，CI对于实验室测试误差影响 要计算步骤如下： 不是很敏感，并没有随岩样的UCS变大而变化2]， (1)选取初始值0。，满足f(0。)>0. Martin等]通过大量实验和统计分析得出启裂强 (2)对于i=1,2,…,n: 度与抗压强度比值C/UCS为0.3~0.5.作为启裂- 1)从转移概率分布f(0·19-1)中产生候选样本 剥落界限模型的输入参数，C/UCS作为一个独立 0°，其中，转移概率函数要满足对称性f(0·10-1)= 参数存在，比值确定的差异与不确定性对于计算结 f0-110°)； 果具有重要的影响：(3)抗拉强度T的确定通常利 2)计算概率密度比 用室内实验和经验公式两种方式，室内实验不确定 r=fa2=y106) 性来源于UCS类同；经验公式法主要以H-B强度准 (6) f(0:-1ly)f(yl8-1)f0:-1) 则[】与Gif价h理论[4]为基础，Murrellt4s]、Cai〔6]、 3)在(0,1)均匀分布间随机产生一个“，使得
工程科学学报,第 41 卷,第 1 期 据. 首先,分析启裂鄄鄄剥落界限模型参数的不确定性 特征,进而提出围岩参数动态概率反演方法,最后, 利用工程案例证明该方法的适用性与可行性. 1 启裂鄄鄄剥落界限模型的不确定性 1郾 1 启裂鄄鄄剥落界限模型 启裂鄄鄄剥落界限模型是 Diederichs [37] 在 H鄄B 强 度准则[38]的基础上提出的,常用于隧道围岩脆性破 坏深度预测[39鄄鄄40] 与层裂破坏风险分析[41] . H鄄B 强 度准则可由公式(1)表示: 滓忆1 = 滓忆3 + UCS ( mb 滓忆3 UCS + s ) a (1) 式中,滓忆1和 滓忆3分别为破坏时的最大和最小有效主应 力;UCS 为完整岩石的单轴抗压强度;mb为完整岩 石的材 料 常 数 mi 的 折 减 值, mb 的 取 值 范 围 为 0郾 0000001 ~ 25;s 为岩石的完整性系数,s 的取值范 围为 0 ~ 1;a 为包络线曲率参数. DISL 模型利用 UCS、CI 及 T 三个参数确定岩体 强度曲线的峰值参数与残余参数[37] ,具体取值见下 表 1. 表 1 DISL 模型参数取值表 Table 1 DISL model input parameters 参数 峰值参数 残余参数 a 0郾 25 0郾 75 s (CI/ UCS) 1 / a 0郾 001 m s(UCS / | T | ) 6 ~ 12 1郾 2 模型参数不确定性来源 由表 1 可知,使用启裂鄄鄄 剥落界限模型需要事 先确定 UCS、CI 及 T 三个参数值,这些参数具有概 率统计特征. 原因如下:(1) UCS 是由完整岩石的 室内单轴抗压强度实验得到,由于岩石中物质成分 的不同与微结构面的存在导致岩石具有各向异性, 岩样制作及实验测量误差等使得岩石单轴抗压强度 具有很大的变异性;(2) 启裂强度 CI 主要是由实验 室声发射试验来确定,CI 对于实验室测试误差影响 不是很敏感,并没有随岩样的 UCS 变大而变化[42] , Martin 等[43]通过大量实验和统计分析得出启裂强 度与抗压强度比值 CI/ UCS 为 0郾 3 ~ 0郾 5. 作为启裂鄄 鄄剥落界限模型的输入参数,CI/ UCS 作为一个独立 参数存在,比值确定的差异与不确定性对于计算结 果具有重要的影响;(3) 抗拉强度 T 的确定通常利 用室内实验和经验公式两种方式,室内实验不确定 性来源于 UCS 类同;经验公式法主要以 H鄄B 强度准 则[38]与 Griffith 理论[44] 为基础,Murrell [45] 、Cai [46] 、 Perras 与 Diederichs [47]等进一步提出了以岩石启裂 强度 CI 为基础的抗拉强度 T 经验确定方法. 因此, 由于实验误差或者模型不确定性等导致完整岩石的 抗拉强度 T 取值存在不确定性. 2 概率反演 B鄄PSO鄄MSVM 方法 2郾 1 贝叶斯概率反演方法 假设隧道工程计算模型可用函数[13]表示为: y^ = g(兹) (2) 式中,兹 为模型的计算参数;y^ 为模型响应的计算值, 如位移、破裂深度等,则实际监测值 y 为: y = y^ + 着 (3) 式中,着 为模型残差,假设 着 服从均值为 0,方差为 滓 2 着 的正态分布,即 着 ~ N(0,滓 2 着 ). 设 兹 为随机变量,此处,兹 为启裂鄄鄄 剥落界限模 型的三个基本参数{UCS、CI/ UCS、T} 的集合,可通 过实验或历史数据统计获得先验概率密度函数 f(兹). 由贝叶斯公式可知,兹 的后验概率密度函数 f(兹 | y)为: f(兹 | y) = 姿f(y | 兹)f(兹) (4) 式中:姿 为归一化系数;f(兹 | y)为似然函数. 对于多源监测数据,如已知 y = (P0 ,P1 ,EDZ), 其中,P0和 P1分别为拱顶下沉点与周边收敛点变化 值, EDZ 为 开 挖 损 伤 区 深 度 ( excavation damage zone,EDZ),假设监测信息间彼此独立,则 f(兹 | y) = f(P0 | 兹)f(P1 | 兹)f(EDZ | 兹) (5) 由于式(4)中后验概率分布是复杂、高维、非标 准的形式,直接积分计算十分困难,马尔科夫链蒙特 卡洛 ( markov chain monte carlo, MCMC ) 模 拟 方 法[16鄄鄄17]是处理复杂统计学问题的有效工具,其核心 思想是通过构造一条马尔科夫链,然后按照转移核 规则引导马尔科夫链扰动过程使其逼近目标分布, 抽取逼近后的样本来近似计算后验分布. MCMC 的 算法很多,Metropolis鄄Hasting [48鄄鄄49] 算法应用广泛,主 要计算步骤如下: (1)选取初始值 兹0 ,满足 f(兹0 ) > 0. (2)对于 i = 1,2,…,n: 1)从转移概率分布 f(兹 * | 兹i - 1 )中产生候选样本 兹 * ,其中,转移概率函数要满足对称性 f(兹 * | 兹i -1 ) = f(兹i - 1 | 兹 * ); 2)计算概率密度比 r = f(兹 * | y) f(兹i - 1 | y) = f(y | 兹 * )f(兹 * ) f(y | 兹i - 1 )f(兹i - 1 ) (6) 3)在(0,1)均匀分布间随机产生一个 u,使得 ·80·
吴忠广等：深埋硬岩隧道围岩参数概率反演方法 81· (0°，r≥u 径向基函数.本文选择性能较好的径向基函数： 0:= (7) (0:-1,r<u /1x-x:12) 4)确定是否收敛：如果不满足要求，重复步骤 K(x,;)=exp2 (14) 1)~3),直至产生稳定序列 核函数参数σ2和惩罚因子C是影响多输出支 对于转移概率分布函数的协方差矩阵取值，参 持向量机建模精度的两个重要因素，为高效确定最 照Zhang等)研究成果，当转移概率分布函数的协 优参数值，采用粒子群搜索的方法优化（σ2，C)的取 方差矩阵为先验概率分布协方差矩阵值的0.5时， 值.粒子群算法[]主要思想是首先初始化一群随 马尔科夫链计算效率较高且能够得到合理的接受 机粒子，然后通过迭代寻找最优解.在每次迭代中， 率.对于模型残差e取值，参照Peng等)]等推荐 粒子通过跟踪两个极值来更新自己.一个是粒子本 值，假设E~N(0,0.25u),其中，u为P。、P1与EDZ 身所找到的最优解，即个体极值P,另一个是整体 三个测试值各自对应的平均值. 种群目前找到的最优解，称为全局极值g·粒子 2.2基于粒子群-多输出支持向量机的智能响应 找到这两个极值后，用下式更新自己的速度和位置： 面方法 V=wV+crand()(PBea -p)+czrand()(gBea-p) 本文采用多输出支持向量机模型建立待反演参数 (15) 与隧道监测数据之间的非线性映射关系，假设有组 P=p+V (16) 训练样本数据{x,}(i=1,2,…,n),x:∈R"为待反演 式中：V是粒子的速度，p为粒子的当前位置，c,和c2 参数，y,∈R为监测数据计算值，可建立如下关系： 是学习因子，通常在0~2间取值，rand()为(0,1) y:R"→R (8) 之间的随机数，为加权因子.针对传统粒子群算 根据多输出支持向量机理论[)，相应的映射模 法易早熟以及算法后期易在全局最优解附近产生振 型与对应的优化问题分别见式(9)与式(10)： 荡的现象，采用权重线性递减的粒子群算法[s0]解决 (x)=W(x)+b (9) 此问题.权重随算法迭代次数的变化公式为： IwI2+C()(10) 0=0m-(0s-0ea) (17) 式中：u:=le:‖=√ee,e=y-(x)rW-b, 式中：w和wn分别表示0的最大值和最小值，t W=[w2,w2,…,w]:b=[b,b2,…,b9] 为当前的迭代步数，tm为最大迭代步数，通常取 0mx=0.9,0min=0.4. 为了得到式(10)对应的最优解，采用迭代加权 在优化过程中，目标函数选用通用的均方差函 最小二乘算法(RWLS)求解2-]，进而得到下式： 数，建立的粒子群-多输出支持向量机模型作为代 L,(W,b)=g(W,b)= 替有限元方法的智能响应面用于MCMC模拟中估 a,4+T (11) 算似然函数 式中：T是与W或b无关的常数项的和； 2.3模型计算过程 概率反演B-PSO-MSVM方法具体计算流程见 0, u:≤B 下图1. a:=2C(u-E) (12) u>8 具体流程如下： 步骤一：通过室内试验与现场测量统计，获取 根据L”(W,b)相对w与驻点条件转化得到 UCS、CI/UCS和T的先验分布信息，如均值、标准 下式： 差、变异系数、分布类型等： K+D j=1,2,…,Q) 步骤二：利用均匀设计建立参数样本集，结合 Phase22数值模拟计算隧道变形及损伤区深度值，构 (13) 建多输出支持向量机学习样本； 式中：D。=diag(a1,a2,…,a1）;=[a,a2,…, 步骤三：利用粒子群方法优化模型参数，通过式 a1]'y=[yny2,…,y];p,=[B,f2,…,B]. (9)建立多输出支持向量机智能响应面模型代替数 其中，K=是核函数矩阵，Φ=[p(x1), 值模型； p(x2),…,p(x)]. 步骤四：将Bayesian方法与多输出支持向量机 常见的核函数分别有线性函数、多项式函数与 模型相结合，基于实际监测数据，构建B-粒子群-
吴忠广等: 深埋硬岩隧道围岩参数概率反演方法 兹i = 兹 * , r逸u 兹{ i - 1 , r < u (7) 4)确定是否收敛;如果不满足要求,重复步骤 1) ~ 3),直至产生稳定序列. 对于转移概率分布函数的协方差矩阵取值,参 照 Zhang 等[7]研究成果,当转移概率分布函数的协 方差矩阵为先验概率分布协方差矩阵值的 0郾 5 时, 马尔科夫链计算效率较高且能够得到合理的接受 率. 对于模型残差 着 取值,参照 Peng 等[17] 等推荐 值,假设 着 ~ N(0, 0郾 25滋),其中,滋 为 P0 、P1与 EDZ 三个测试值各自对应的平均值. 2郾 2 基于粒子群鄄鄄 多输出支持向量机的智能响应 面方法 本文采用多输出支持向量机模型建立待反演参数 与隧道监测数据之间的非线性映射关系,假设有 n 组 训练样本数据{xi,y^ i}(i =1,2,…,n),xi沂R n 为待反演 参数,yi沂R 为监测数据计算值,可建立如下关系: y^:R n寅R (8) 根据多输出支持向量机理论[32] ,相应的映射模 型与对应的优化问题分别见式(9)与式(10): y^(x) = W T渍(x) + b (9) Lp (W,b) = 1 2 移 Q j = 1 椰w j椰2 + C 移 l i = 1 L(ui) (10) 式中:ui = 椰ei椰 = e T i ei,e T i = y T i - 渍( xi ) TW - b T , W = [w 1 ,w 2 ,…,w Q ];b = [b 1 ,b 2 ,…,b Q ] T . 为了得到式(10)对应的最优解,采用迭代加权 最小二乘算法(IRWLS)求解[32鄄鄄33] ,进而得到下式: Lp (W,b) = L义p (W,b) = 1 2 移 Q j = 1 椰w j椰2 + 1 2 移 l i = 1 aiu 2 i + 子 (11) 式中:子 是与 W 或 b 无关的常数项的和; ai = 0, ui臆着 2C(ui - 着) ui , ui > ì î í ïï ïï 着 (12) 根据 L义p (W,b)相对 w j 与 b j 驻点条件转化得到 下式: K + D - 1 琢 1 琢 T 1 é ë ê ê ù û ú T ú 茁 j b é ë ê ê ù û ú j ú = y j 琢 T y é ë ê ê ù û ú j ú (j = 1,2,…,Q) (13) 式中:D琢 = diag ( 琢1 , 琢2 ,…, 琢1 ); 琢 = [ 琢1 , 琢2 , …, 琢1 ] T ;yj = [yj1 ,yj2 ,…,yj1 ] T ;茁j = [茁j1 ,茁j2 ,…,茁j1 ]. 其中, K = 椎椎 T 是核函数矩阵,椎 = [渍( x1 ), 渍(x2 ),…, 渍(xl)] T . 常见的核函数分别有线性函数、多项式函数与 径向基函数. 本文选择性能较好的径向基函数: K(x,xi) = exp { - | x - xi | 2 滓 2 } (14) 核函数参数 滓 2和惩罚因子 C 是影响多输出支 持向量机建模精度的两个重要因素,为高效确定最 优参数值,采用粒子群搜索的方法优化(滓 2 ,C)的取 值. 粒子群算法[35] 主要思想是首先初始化一群随 机粒子,然后通过迭代寻找最优解. 在每次迭代中, 粒子通过跟踪两个极值来更新自己. 一个是粒子本 身所找到的最优解,即个体极值 pBest,另一个是整体 种群目前找到的最优解,称为全局极值 gBest . 粒子 找到这两个极值后,用下式更新自己的速度和位置: V = wV + c1 rand()(pBest - p) + c2 rand()(gBest - p) (15) p = p + V (16) 式中:V 是粒子的速度,p 为粒子的当前位置,c1和 c2 是学习因子,通常在 0 ~ 2 间取值,rand( )为 (0,1) 之间的随机数,w 为加权因子. 针对传统粒子群算 法易早熟以及算法后期易在全局最优解附近产生振 荡的现象,采用权重线性递减的粒子群算法[50] 解决 此问题. 权重随算法迭代次数的变化公式为: w = wmax - t(wmax - wmin ) tmax (17) 式中:wmax 和 wmin分别表示 w 的最大值和最小值,t 为当前的迭代步数, tmax 为最大迭代步数,通常取 wmax = 0郾 9,wmin = 0郾 4. 在优化过程中,目标函数选用通用的均方差函 数,建立的粒子群鄄鄄 多输出支持向量机模型作为代 替有限元方法的智能响应面用于 MCMC 模拟中估 算似然函数. 2郾 3 模型计算过程 概率反演 B鄄PSO鄄MSVM 方法具体计算流程见 下图 1. 具体流程如下: 步骤一:通过室内试验与现场测量统计,获取 UCS、CI/ UCS 和 T 的先验分布信息,如均值、标准 差、变异系数、分布类型等; 步骤二:利用均匀设计建立参数样本集,结合 Phase2 数值模拟计算隧道变形及损伤区深度值,构 建多输出支持向量机学习样本; 步骤三:利用粒子群方法优化模型参数,通过式 (9)建立多输出支持向量机智能响应面模型代替数 值模型; 步骤四:将 Bayesian 方法与多输出支持向量机 模型相结合,基于实际监测数据,构建 B鄄鄄 粒子群鄄鄄 ·81·
.82 工程科学学报，第41卷，第1期 3.2基于均匀设计方法的训练样本集确定 确定UCS、CI/UCS、T的统计特征 步骤1 (均值、标准差、变异系数、分布类型等) 将UCS、CI/UCS与T三个参数作为随机变量， 在参数三倍标准差范围内取值，即UCS取值范围为 步骤2 利用均匀设计构造测试参数样本集 Phase2数值模型 [uues -30ucs,uucs +3ucs]=[50,170],CI/UCS 值范围为[as-3 Gcvues,Hevucs+30cwCs]= PSO优化 [0.337,0.469],T取值范围为[u-307,r+ 步骤3 3o]=[-1.5,-9.9],利用均匀设计方法[s2]构 建立MSVM响应面 构建MSVM学习样本 建UCS、C/UCS与T的参数样本集，选取U,(372) 构造输入训练样本集，并计算拱顶下沉点P。和周边 步骤4 构建B-PSO-MSVM模型 实际监测值 收敛点P,变化值，与开挖损伤区深度EDZ作为训练 输出样本，具体见表3.以试验序号19参数值为例， ICS.CVUCS、T动态更新， 利用Phase2软件建立算例模型见图2，计算结果见 步骤5 估测参数不确定性 表3训练样本集 ------- Table 3 Set of training samples 结束 试验 样本输人 样本输出 图1B-PSO-MSVM计算流程 UCS/ 序号 CI Po/ P/ EDZ/ Fig.1 Calculation flowchart of B-PSO-MSVM method MPa UCS MPa mm mm m 50 0.3722 -6.5413.513.5 3.808 多输出支持向量机概率反演方法，利用MCMC算法 54 0.414 -3.4612 3.263 实现参数动态更新： 58 0.4558 -8.78 15 13.5 2.966 60 0.359 -5.42 9.5 12 2.963 步骤五：确定UCS、C/UCS和T的后验分布信 62 0.3986 -2.2 9 12 2.898 息，进行参数不确定性分析 6 66 0.4382 -7.66 8.95 12 2.066 个 70 0.3458 -4.3 10.5 13.5 2.948 3算例分析 8 74 0.3854 -9.62 8.5 12 2.052 9 78 0.425 -6.4 1.35 9.5 0.556 3.1工程背景 10 80 0.4646 -3.18 8.8 12 1.929 计算模型为圆形地下试验隧道[39,1]，隧道直径 11 82 0.37 -8.5 95 12 2.027 86 0.4118 -5.14 8.5 12 1.842 为6.5m,长度为46m,埋深约680m.沿整个隧道长 13 90 0.4514 -2.06 7.45 10.5 1.791 度上，岩体由石灰岩组成，围岩坚硬完整，表现出高 94 0.3546 -7.38 12 g 2.252 地应力特点.经现场地应力实测，最大主应力为33 15 98 0.3942 -4.02 8.8 12 1.825 16 100 0.436 -9.34 12 1.788 MPa,中间主应力为17MPa,最小主应力为15MPa, 17 102 0.3414 -5.14 9.5 12 1.827 中间主应力与隧道轴线近似平行，最大与最小主应 18 106 0.381 -2.9 10.5 12 1.787 力位于与洞轴线垂直的平面内，其中，最大主应力方 9 110 0.4206 -8.22 12 10.5 1.378 2 114 0.4602 -5 8.55 10.5 1.236 向按近似水平计算，岩体力学参数按照Langford与 118 0.3678 -1.78 22 11.5 1.423 Diederichs【a]、吴成与张平Is]案例取值，具体见表2. 22 120 0.4074 -7.1 8.75 10.5 1.774 23 122 0.447 -3.74 10.5 10.5 1.367 表2岩体力学参数 24 126 0.3502 -9.2 1.785 Table 2 Rock mass mechanical parameters 25 130 0.392 -5.98 6.3 10.5 1.375 26 134 0.4338 -2.62 2 10 0.718 参数 取值 27 138 0.337 -7.94 12 1.399 地质强度指标，GSI 90 28 140 0.3766 -4.86 8.1 12 1.772 岩石单轴抗压强度，UCS/MPa 1ogN(110,20) 29 142 0.4162 -1.5 7.65 10.5 1.111 启裂强度/抗压强度，C/UCS N(0.403,0.022) 30 146 0.458 -6.82 9.6 10 0.868 岩石抗拉强度，T/MPa N(-5.7.1.4) 31 150 0.3634 -3.6 10.5 1.235 残余m值，me 32 154 0.403 -9.06 15 10.5 1.282 岩石弹性模量，E/GPa 19 33 158 0.4426 -5.7 4.35 10 0.675 160 0.348 8.55 10.5 0.866 泊松比，” 34 -2.34 0.12 35 162 0.3898 -7.8 8.95 10.5 1.123 注：N(4,,G,)指服从均值为41，标准差为σ，的正态分布：logN 36 166 0.4294 -4.58 3 9.5 0.674 (山，02)指服从均值为山2，标准差为2的对数正态分布. 37 170 0.469 -9.9 3.5 9 0.556
工程科学学报,第 41 卷,第 1 期 图 1 B鄄PSO鄄MSVM 计算流程 Fig. 1 Calculation flowchart of B鄄PSO鄄MSVM method 多输出支持向量机概率反演方法,利用 MCMC 算法 实现参数动态更新; 步骤五:确定 UCS、CI/ UCS 和 T 的后验分布信 息,进行参数不确定性分析. 3 算例分析 3郾 1 工程背景 计算模型为圆形地下试验隧道[39,51] ,隧道直径 为 6郾 5 m,长度为 46 m,埋深约 680 m. 沿整个隧道长 度上,岩体由石灰岩组成,围岩坚硬完整,表现出高 地应力特点. 经现场地应力实测,最大主应力为 33 MPa,中间主应力为 17 MPa,最小主应力为 15 MPa, 中间主应力与隧道轴线近似平行,最大与最小主应 力位于与洞轴线垂直的平面内,其中,最大主应力方 向按近似水平计算,岩体力学参数按照 Langford 与 Diederichs [41] 、吴成与张平[51]案例取值,具体见表2. 表 2 岩体力学参数 Table 2 Rock mass mechanical parameters 参数 取值 地质强度指标,GSI 90 岩石单轴抗压强度,UCS / MPa logN(110,20) 启裂强度/ 抗压强度,CI/ UCS N(0郾 403,0郾 022) 岩石抗拉强度,T / MPa N( - 5郾 7,1郾 4) 残余 m 值,mres 8 岩石弹性模量,E/ GPa 19 泊松比,淄 0郾 12 注:N( 滋1 ,滓1 )指服从均值为 滋1 ,标准差为 滓1 的正态分布;logN (滋2 ,滓2 )指服从均值为 滋2 ,标准差为 滓2的对数正态分布. 3郾 2 基于均匀设计方法的训练样本集确定 将 UCS、CI/ UCS 与 T 三个参数作为随机变量, 在参数三倍标准差范围内取值,即 UCS 取值范围为 [滋UCS - 3滓UCS , 滋UCS + 3滓UCS ] = [50,170],CI/ UCS 取 值范围 为 [ 滋CI/ UCS - 3滓CI/ UCS , 滋CI/ UCS + 3滓CI/ UCS ] = [0郾 337, 0郾 469 ], T 取值范围为 [ 滋T - 3滓T , 滋T + 3滓T ] = [ - 1郾 5, - 9郾 9],利用均匀设计方法[52] 构 建 UCS、CI/ UCS 与 T 的参数样本集,选取 U37 (37 12 ) 构造输入训练样本集,并计算拱顶下沉点 P0和周边 收敛点 P1变化值,与开挖损伤区深度 EDZ 作为训练 输出样本,具体见表 3. 以试验序号 19 参数值为例, 利用 Phase2 软件建立算例模型见图 2,计算结果见 表 3 训练样本集 Table 3 Set of training samples 试验 序号 样本输入 样本输出 UCS / MPa CI/ UCS T / MPa P0 / mm P1 / mm EDZ / m 1 50 0郾 3722 - 6郾 54 13郾 5 13郾 5 3郾 808 2 54 0郾 414 - 3郾 46 12 12 3郾 263 3 58 0郾 4558 - 8郾 78 15 13郾 5 2郾 966 4 60 0郾 359 - 5郾 42 9郾 5 12 2郾 963 5 62 0郾 3986 - 2郾 2 9 12 2郾 898 6 66 0郾 4382 - 7郾 66 8郾 95 12 2郾 066 7 70 0郾 3458 - 4郾 3 10郾 5 13郾 5 2郾 948 8 74 0郾 3854 - 9郾 62 8郾 5 12 2郾 052 9 78 0郾 425 - 6郾 4 1郾 35 9郾 5 0郾 556 10 80 0郾 4646 - 3郾 18 8郾 8 12 1郾 929 11 82 0郾 37 - 8郾 5 9郾 5 12 2郾 027 12 86 0郾 4118 - 5郾 14 8郾 5 12 1郾 842 13 90 0郾 4514 - 2郾 06 7郾 45 10郾 5 1郾 791 14 94 0郾 3546 - 7郾 38 12 14 2郾 252 15 98 0郾 3942 - 4郾 02 8郾 8 12 1郾 825 16 100 0郾 436 - 9郾 34 12 12 1郾 788 17 102 0郾 3414 - 5郾 14 9郾 5 12 1郾 827 18 106 0郾 381 - 2郾 9 10郾 5 12 1郾 787 19 110 0郾 4206 - 8郾 22 12 10郾 5 1郾 378 20 114 0郾 4602 - 5 8郾 55 10郾 5 1郾 236 21 118 0郾 3678 - 1郾 78 22 11郾 5 1郾 423 22 120 0郾 4074 - 7郾 1 8郾 75 10郾 5 1郾 774 23 122 0郾 447 - 3郾 74 10郾 5 10郾 5 1郾 367 24 126 0郾 3502 - 9郾 2 18 12 1郾 785 25 130 0郾 392 - 5郾 98 6郾 3 10郾 5 1郾 375 26 134 0郾 4338 - 2郾 62 2 10 0郾 718 27 138 0郾 337 - 7郾 94 15 12 1郾 399 28 140 0郾 3766 - 4郾 86 8郾 1 12 1郾 772 29 142 0郾 4162 - 1郾 5 7郾 65 10郾 5 1郾 111 30 146 0郾 458 - 6郾 82 9郾 6 10 0郾 868 31 150 0郾 3634 - 3郾 6 12 10郾 5 1郾 235 32 154 0郾 403 - 9郾 06 15 10郾 5 1郾 282 33 158 0郾 4426 - 5郾 7 4郾 35 10 0郾 675 34 160 0郾 348 - 2郾 34 8郾 55 10郾 5 0郾 866 35 162 0郾 3898 - 7郾 8 8郾 95 10郾 5 1郾 123 36 166 0郾 4294 - 4郾 58 2 9郾 5 0郾 674 37 170 0郾 469 - 9郾 9 3郾 5 9 0郾 556 ·82·