一种利用向量投影响应面法确定隧道顶板失效概率的方法

1.本发明涉及隧道工程技术领域，尤其涉及一种利用向量投影响应面法确定隧道顶板失效概率的方法。


背景技术：

2.伴随着地下空间的快速发展，隧道失稳破坏的事故越来越多，而隧道顶板的失稳垮塌是施工中常见的灾害形式。由于隧道所处的不良地质环境、自然环境等的影响，岩体被赋予天然的变异性和不确定性，在建设过程中，导致隧道顶板楔形体易发生变形破坏、失稳坍塌等事故，将对人身安全、建设发展等带来重大危害。为了隧道建设的稳定、安全和高质量的发展，进一步评估隧道顶板楔形体稳定性势在必行。在目前隧道工程的施工中多采用确定性分析方法来评估隧道顶板楔形体的稳定性，忽略了岩土体材料参数固有的不确定性，常常会出现在结构安全系数的合理范围内，仍出现隧洞坍塌的现象。故在隧道工程的设计、施工中应尽可能多地考虑岩体力学参数对隧道围岩稳定性的影响。近年来，发展起来的可靠度分析方法为充分考虑岩体力学参数不确定性来评估隧道稳定性提供了可能，因此，可以采用可靠性分析方法探究相关影响因素对结构稳定性的影响，并计算结构失效概率，为隧道建设提供安全评估的方法。对于隧道顶板楔形体，其受地应力、围岩应力、不良岩体结构及洞室断面形状和尺寸等多种因素影响，变形破坏过程较为复杂，准确评估稳定较为困难。因此，深入研究隧道顶板楔形体的稳定性对于研究隧道工程的稳定发展是非常有意义的。
3.目前，隧道工程可靠度分析方法主要包括：直接积分法、蒙特卡罗模拟法、可靠度指标法、响应面法等。其中可靠度指标法包括一次二阶矩、二次二阶矩、高阶矩法等，对于低阶且随机变量数量少的情况，矩法假定结构的极限状态方程服从某一分布状态，根据功能函数z，通过随机变量x的均值和标准差比值确定可靠度指标β，再将可靠度指标代入失效概率公式中确定失效概率pf，其概念清晰、计算简便、计算精度基本能满足一般的工程要求。但是对于高阶矩法，在评估隧道结构稳定性的过程中，需要计算各随机变量的偏导，这无疑大大提高了对数理统计掌握不好的工程师运用该方法的门槛。


技术实现要素：

4.为了克服由于隧道顶板楔形体不同地质参数和几何参数存在的不确定性，导致计算出的安全系数远高于临界值时，仍然会出现部分隧洞出现坍塌破坏的问题。本发明提供了一种利用向量投影取样点的响应面法计算隧道顶板失效概率的方法，目的是为了解决现有技术中存在的缺点。
5.为了实现上述目的，本发明提供如下技术方案：一种利用向量投影响应面法确定隧道顶板失效概率的方法，包括如下步骤：
6.通过极限平衡法确定隧道顶板楔形体的破坏模式及其安全系数fs；
7.根据安全系数fs确定相关地质参数和几何参数x，将地质参数和几何参数x按照服
从正态分布xi(μ,σ)的形式赋予不确定性；
8.设定初始迭代点x，确定采样点范围控制参数f；
9.通过中心复合设计法生成展开点x并生成试验点，确定出试验点中的未知量；
10.通过基于向量投影取样点的响应面法计算新的展开点x，并用代替f；
11.以新的展开点x重复获取试验点x中的未知量与计算新的展开点x的步骤，直至前后两次||x||之差小于允许误差εq，确定展开点x；
12.将确定的展开点x输入功能函数中，确定可靠度指标β；
13.将可靠度指标β输入失效概率公式得到隧道顶板楔形体失效概率pf。
14.优选的，所述功能函数可以定义为：z＝g(x)＝fs-1，其中，x＝[x1,
…
xi…
,xn]
t
为相关地质参数和几何参数，其中i＝1,2,3,
…
,n。
[0015]
优选的，所述安全系数fs的计算公式为：
[0016][0017]
其中，s为切向力，n为法向力，w为楔形体的重力，α为隧道楔形体半顶角。
[0018]
优选的，所述确定试验点x中的未知量包括如下步骤：
[0019]
通过数值分析或实验，计算在展开点处的功能函数估计值
[0020]
计算功能函数估计值的系数λ与系数矩阵a，进而计算响应面函数的系数a0、a＝(a1,a2,
···
,an)
t
；
[0021]
计算结构可靠度指标β和设计点x
*
。
[0022]
优选的，所述计算功能函数估计值的系数λ与系数矩阵a的计算公式为：
[0023][0024]
式中，λ＝(a0,a1,a2,
···
,an)
t
，a为系数矩阵；
[0025]
若是按照展开点(x1,x2,
…
,xn)
t
，(x1,
…
,x
i-fσ
xi
,
…
,xn)
t
，以及(x1,
…
,xi+fσxi,
…
,xn)
t
的顺序计算的，其中i＝1,2,
…
,n，矩阵a为：
[0026][0027]
式中通过下式求得,并通过该式确定可靠度指标β及设计点x
*
：
[0028][0029]
式中，a0，a＝(a1,a2,
···
,an)
t
为系数；
[0030]
通过下式计算可靠度指标β：
[0031][0032]
其中，为岩体力学参数的均值，为岩体力学参数xi的标准差。
[0033]
优选的，所述计算新的展开点x的计算公式为：
[0034][0035]
式中，代表取样点与设计点的距离；ri为向量qi的第i个分量；表示在xi处的标准差；ei＝(δ
1i
,δ
2i
,
···
,δ
ni
)
t
为沿xi坐标轴的单位基向量，其中δ
ij
为kroneckerδ，即将n维零向量0中的第i个元素替换为1所得到的向量；qi表示单位投影向量。
[0036]
优选的，所述单位投影向量qi的计算公式为：
[0037]
qi＝qi/||qi||(i＝1,2,
···
,n)
[0038]
其中qi为：
[0039]
qi＝ω(ti+εqαr)+(1-w)ei，i＝1,2,
···
,n
[0040]
式中，εq为一个小数；0≤ω≤1是加权因子，ω＝1相应于全向量投影，ω＝0相应于基本响应面方法中的中心复合设计取样；αr为响应面函数在x
*
处的单位向量，公式如下；ti为第i个随机变量xi的取样点投影用到的单位向量：ti＝ti/||ti||(i＝1,2,
···
,n)，其中ti公式如下：
[0041][0042][0043]
式中，α
ri
为向量ar的第i个分量。
[0044]
优选的，所述述岩体力学参数包括节理摩擦角、楔形体的中心高度、楔形体半顶角、侧压力系数、垂直应力、容重、圆形隧道半径、节理的切向和法向刚度比。
[0045]
优选的，所述隧道顶板楔形体失效概率pf的公式为：
[0046]
pf＝φ(-β)＝1-φ(β)
[0047]
其中φ(β)查标准正态分布表可知。
[0048]
优选的，利用所述向量响应面法计算展开点之前，需进行如下步骤：
[0049]
统计所述影响隧道顶板楔形体失效概率相关参数的均值和标准差；
[0050]
根据所述相关参数的均值和标准差产生服从正态概率分布类型。
[0051]
本发明与现有技术相比具有以下有益效果：
[0052]
1、本发明的向量投影取样点的响应面法通过响应面的梯度投影选取取样点，以线性响应面函数为基础，利用一次二阶矩方法搜索设计点并计算可靠度指标，可以将非线性较高的真实的功能函数或极限状态面，用一个容易处理的函数或曲面代替，将其应用于计算失效概率的稳定性中，具有较高的准确性，并显著提高分析效率。
[0053]
2、与传统的计算隧道顶板楔形体稳定性的方法相比，本发明中的向量投影响应面
法将计算机技术引入到岩土工程中，在分析精度、广度、效率上都得到了明显提升。
附图说明
[0054]
图1为本发明实施例提供的利用向量投影取样点的响应面法确定隧道顶板楔形体失效概率方法流程图。
[0055]
图2为本发明实施例提供的隧道顶板楔形体破坏模式示意图。
[0056]
图3为本发明实施例提供的利用向量投影取样点的响应面法隧道顶板楔形体失效概率与蒙特卡罗随机抽样法mcs的隧道顶板失效概率方法的对比验证图。
具体实施方式
[0057]
下面结合附图，本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0058]
直接积分法是根据结构可靠度的原始定义发展而来的计算方法，需要采用多重数值积分求解系统结构的失效概率，导致在可靠度分析过程中计算繁琐，因此，在隧道工程中的应用相对较少。蒙特卡罗模拟法在抽样的样本足够大时计算的精度非常高，但是蒙特卡罗模拟法的主要缺点是为了达到满足要求的计算精度，往往需要大量的随机样本，从而导致可靠度分析的计算量非常大。可靠度指标法对于低阶且随机变量数量少的情况，其概念清晰、计算简便、计算精度基本能满足一般的工程要求。但是对于高阶矩法，在评估隧道结构稳定性的过程中，需要计算各随机变量的偏导，这无疑大大提高了对数理统计掌握不好的工程师运用该方法的门槛。
[0059]
为此，本发明提出了一种利用向量投影取样点的响应面方法，通过响应面的梯度投影选取取样点，以线性响应面函数为基础，利用一次二阶矩方法搜索设计点并计算可靠度指标，可以将非线性较高的真实的功能函数或极限状态面，用一个容易处理的函数或曲面代替，将其应用于计算失效概率的稳定性分析中，具有较高的准确性，并显著提高分析效率和准确性也更高。。
[0060]
参见图1～3，下面详细对一种利用向量投影响应面法确定隧道顶板楔形体失效概率的方法进行说明，该方法包括：：
[0061]
(1)功能函数可以定义为：z＝g(x)＝fs-1
[0062]
其中，x＝[x1,
…
,xn]
t
(i＝1,2,3,
…
,n)，为相关地质参数和几何参数；fs为隧道顶板楔形体失稳的安全系数。
[0063]
(2)获取相关地质参数和几何参数x随机量的统计特征，即服从正态分布xi(μ,σ)。
[0064]
(3)确定展开点取值，步骤为：
①
假定初始迭代点x，一般可取x＝μ
x
，确定系数f并展开点x生成试验点；
②
确定试验点x中的未知量：功能函数估计值系数λ、系数矩阵a、结构可靠度指标β和设计点x
*
；
③
计算新的展开点x并用代替f；
④
以新的展开点x重复步骤
②
、
③
，直至前后两次||x||之差小于允许误差εq，即可得到确定的展开点x。
[0065]
(4)将确定的展开点x代入功能函数：z＝g(x)＝fs-1。
[0066]
(5)将确定的结构可靠度指标β带入失效概率公式，即可得到隧道顶板楔形体失效概率pf：
[0067]
pf＝φ(-β)＝1-φ(β)
[0068]
其中φ(β)查标准正态分布表可知。
[0069]
具体实验数据和计算公式如下：
[0070]
(1)隧道顶板楔形体的地质参数和几何参数相关参数由现场或工地实验室获取。
[0071]
(2)隧道顶板楔形体的地质参数和几何参数相关参数包括：节理摩擦角楔形体的中心高度h、楔形体半顶角α、侧压力系数k0、垂直应力p、容重γ、节理的切向和法向刚度比ks/kn。
[0072]
(3)隧道顶板楔形体破坏的模式：以极限平衡法为理论基础的brady经典隧道顶板楔形体经典破坏模式；tonon在brady的基础上，提出将安全系数定义为被动力(有效抗剪强度)和主动力(抗剪强度)之比。
[0073]
(4)安全系数fs的计算公式为：
[0074][0075]
其中，s为切向力，n为法向力，w为楔形体的重力，α为隧道楔形体半顶角。
[0076]
安全系数fs中的n和w经由以下公式推导而来：
[0077][0078][0079][0080][0081][0082][0083][0084]
其中，h0为楔形体内部的水平力，h为楔形体的中心高度，k0为侧压力系数，p为垂直应力，γ为容重，ks/kn为节理的切向和法向刚度比。
[0085]
(5)利用向量投影取样点的响应面法确定功能函数为：z＝g(x)＝fs-1，其相应的功能函数为：
[0086][0087]
(6)基于向量投影取样点的响应面法计算设计点的计算如下：
[0088]
①
假定初始迭代点x，一般取其平均值为x＝μ
x
。
[0089]
②
选取f值，考虑到实际功能函数z的非线性，在选取f时可以先用f＝1.0和1.5进行试算，根据所得到的相应的可靠指标β
1.0
和β
1.5
，如定义δ＝|β
1.0-β
1.5
|/(n-1)，来计及z的非线性效应。如果δ较小，例如δ≤0.03，可以采用f＝1.0；如果0.03≤δ≤0.06，可以采用f＝1.2。一般f在1.0～2.0之间选取。
[0090]
③
利用试验设计展开点x生成试验点。
[0091]
④
通过结构的数值分析或试验产生在试验点上功能函数的估计值
[0092]
⑤
通过下式求解系数λ：
[0093][0094]
式中，λ＝(a0,a1,a2,
···
,an)
t
，a为系数矩阵。
[0095]
若是按照展开点(x1,x2,
…
,xn)
t
，(x1,
…
,x
i-fσ
xi
,
…
,xn)
t
(i＝1,2,
…
,n)，以及(x1,
…
,xi+fσ
xi
,
…
,xn)
t
(i＝1,2,
…
,n)的顺序计算的，矩阵a为：
[0096][0097]
式中通过下式求得,并通过该式确定可靠度指标及设计点x
*
：
[0098][0099]
式中，a0，a＝(a1,a2,
···
,an)
t
为系数；
[0100]
通过下式计算可靠度指标β：
[0101][0102]
其中，为岩体力学参数的均值，为岩体力学参数xi的标准差。
[0103]
⑥
计算新的展开点x以及以来代替f：
[0104]
展开点x计算公式为：
[0105][0106]
式中，代表取样点与设计点的距离；ri为向量qi的第i个分量；表示在xi处的标准差；ei＝(δ
1i
,δ
2i
,
···
,δ
ni
)
t
为沿xi坐标轴的单位基向量，其中δ
ij
为kroneckerδ，即将n维零向量0中的第i个元素替换为1所得到的向量；qi表示单位投影向量：qi＝qi/||qi||(i＝1,2,
···
,n)；其中qi为：
[0107]
qi＝ω(ti+εqαr)+(1-w)ei，i＝1,2,
···
,n
[0108]
式中，εq为一个小数；0≤ω≤1是加权因子，ω＝1相应于全向量投影，ω＝0相应于基本响应面方法中的中心复合设计取样；αr为响应面函数在x
*
处的单位向量，公式如下；ti为第i个随机变量xi的取样点投影用到的单位向量：ti＝ti/||ti||(i＝1,2,
···
,n)，其中ti公式如下。
[0109][0110][0111]
式中，α
ri
为向量ar的第i个分量。
[0112]
以新的展开点x重复步骤
③
～
⑥
，直至前后两次||x||之差小于允许误差εq。
[0113]
下面通过实施例结合附图进一步具体说明本发明的步骤：
[0114]
(1)根据图1的brady经典隧道顶板楔形体经典破坏模式，tonon提出将安全系数fs定义为被动力和主动力之比，其中，被动力为有效抗剪强度，主动力为抗剪强度，确定影响隧道顶板楔形体稳定性的地质参数和几何参数：节理摩擦角(φ)、楔形体的中心高度(h)、楔形体半顶角(α)、侧压力系数(k0)、垂直应力(p)、容重(γ)、节理的切向和法向刚度比(ks/kn)等岩体力学为研究的随机变量，以上相关参数均可以通过现场或工地实验室获得，整理出各个岩体力学参数的均值(μ)和标准差(σ)。具体的，岩土力学参数的取值范围如表1所示：
[0115]
表1不同岩土力学参数的取值范围
[0116][0117]
(2)为考虑各岩体力学参数对隧道顶板楔形体稳定性的影响，取各参数的均值设为初始迭代点x，将其输入功能函数z，确定系数f，同时计算新的展开点x与计算公式中的未知量：结构可靠度指标β，以新的展开点x重复上述步骤，直至误差小于εq＝1
×
10-6
，确定最终展开点x。
[0118]
初始迭代点x：
[0119]
x＝{35,5.1,25,0.5,0.5,0.027,0.1}
[0120]
在选取f时，先用f＝1.0和1.5进行试算，所得到的相应的可靠指标β
1.0
＝0.4628和β
1.5
＝0.3571，则δ＝|β
1.0-β
1.5
|/(n-1)＝0.035，0.03≤δ≤0.06，来计及z的非线性效应。如果δ较小，例如δ≤0.03，可以采用f＝1.0；如果0.03≤δ≤0.06，因此采用f＝1.2，一般f在1.0～2.0之间选取。
[0121]
结构可靠度指标β为：
[0122]
β＝0.4219
[0123]
最终展开点x：
[0124]
x＝{35.002,5.1226,25.491,0.4601.0.4956,0.0271,0.0996}
[0125]
将结构可靠度指标β输入失效概率公式，即可得到隧道顶板楔形体失效概率：
[0126]
pf＝1-φ(β)＝0.3365
[0127]
响应面法计算完成所用时间t：
[0128]
t＝0.007315s
[0129]
(3)采用matlab软件对服从正态分布的各相关参数生成随机场，然后，代入安全系数公式fs，统计fs《1的样本数量n，并利用蒙特卡罗的方法，计算失效概率pf。不同样本n下的隧道顶板楔形体失效概率见表2，在图3中可以清晰地看出样本量n变化下失效概率pf的变化。当n＝10时，n＝4，pf＝4/10＝0.4；当n＝100时，n＝33，pf＝33/100＝0.33；当n＝1
×
104时，n＝3394，pf＝3394/1
×
104＝0.3394；当n＝1
×
106时，n＝339700，pf＝339700/1
×
106＝0.3397。从图中可以看出当样本数量n较小时，pf会有一定的波动，当样本量n逐渐增大时，pf会逐渐收敛于固定的值，本发明的实施例中，p
fmcs
＝33.97％。
[0130]
表2不同样本量对应的隧道顶板失效概率
[0131][0132]
将利用向量投影取样点的响应面法和蒙特卡罗随机抽样的失效概率进行对比，两者之间的失效概率误差在仅为0.32％，说明了本发明专利提供的方法具有很高的精度，并且相较于蒙特卡罗计算更加简便、快捷，进一步验证了该方法计算隧道顶板楔形体的高效性和准确性。
[0133]
以上所述实施例仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的
简单变化或等效替换，均属于本发明的保护范围。