相关随机变量拉丁超立方抽样及隧道掌子面稳定性概率分析_章从旭.pdf

第 42 卷第 3 期重 庆 交 通 大 学 学 报(自 然 科 学 版)Vol 42No32023 年 3 月JOUNAL OF CHONGQING JIAOTONG UNIVESITY(NATUAL SCIENCE)Mar 2023DOI:103969/jissn1674-069620230307相关随机变量拉丁超立方抽样及隧道掌子面稳定性概率分析章从旭1，蓝元盛2，李斌2(1 中交路桥建设有限公司，北京 100027;2 武汉理工大学 交通与物流工程学院，湖北 武汉 430063)摘要:总结了4 种不同类型的相关随机变量拉丁超立方抽样原理，通过 Python 语言实现各种情况的相关随机变量拉丁超立方抽样;基于拉丁超立方抽样的样本，采用 FLAC3D软件，对隧道掌子面稳定性进行确定性分析，统计了失效概率。结果表明:采用拉丁超立方抽样可以在满足计算精度、提高计算效率的同时，大大减少不确定性分析所需的样本数量。关键词:隧道工程;岩土工程;相关随机变量;标准正态分布中图分类号:U452;TU432文献标志码:A文章编号:1674-0696(2023)03-052-07Latin Hypercube Sampling for elevant amdom Variables andProbabilistic Analysis of Tunnel Face StabilityZHANG Congxu1，LAN Yuansheng2，LI Bin2(1 oad Bridge International Co，Ltd，Beijing 100027，China;2 School of Transportation and Logistics Engineering，Wuhan University of Technology，Wuhan 430063，Hubei，China)Abstract:The Latin hypercube sampling principles of four different types of relevant random variables were summarized andLatin hypercube sampling of relevant variables in various situations was realized through Python language FLAC3Dwasutilized to perform deterministic analysis of tunnel face stability and calculate the failure probability based on the obtainedsamples generated by Latin hypercube sampling The results show that Latin hypercube sampling can greatly reduce thenumber of samples required for uncertainty analysis while meeting the required computational accuracy and improvingcomputational efficiencyKey words:tunnel engineering;geotechnical engineering;relevant random variables;standard normal distribution0引言随着可靠度理论的不断发展，可靠度设计在岩土工程的应用越来越广。可靠度设计方法包括一次二阶矩法、一阶可靠度法、响应面法、数据表法和蒙特卡洛法等。其中，蒙特卡洛法是最为直接和准确的方法1-2。蒙特卡洛法采集的样本可通过有限元计算得到失效概率，分析结果可作为其它方法计算结果的唯一验证3-9。然而，蒙特卡洛法的计算成本较高，对于失效概率较小的问题，有时几乎不可能实现。M D MCKAY 等10 在 1979 年指出拉丁超立方抽样技术可降低计算消耗，该技术具有分层均匀、抽样少的优点。具体计算步骤:确定抽样样本数 N;将变量的累积分布分成相同的 N 个小区间，即将变量累积分布(0，1)区间均分成 N 段，使得每个区间有相同的概率;从 N 个区间段中分别随机抽取一个值，将抽取的 N 个值分别通过变量的分布反函数映射为变量的样本，共计得到 N 个样本。随机变量的拉丁超立方抽样原理如图 1。图 1(a)表示将样收稿日期:2021-09-19;修订日期:2023-02-12基金项目:中交路建科技研发项目(ZJLJ-2019-14);国家自然科学基金项目(51608407)第一作者:章从旭(1974)，男，四川自贡人，高级工程师，主要从事大跨径桥梁及山区高速公路施工方面的研究。E-mail:452041358 qqcom通信作者:蓝元盛(1993)，男，江西赣州人，硕士研究生，主要从事隧道工程掌子面稳定性方面的研究。E-mail:yuans_l163com本的概率分布范围划分为 5 个概率相等的区间，每个区间采集一个样本点;图 1(b)表示 2 个独立样本经分层抽样后的组合情况，2 个变量在任何一个分割的概率区间只出现一次，以保证所采集的样本最具代表性。拉丁超立方抽样保证了尾部层中的样本也会被抽取，而这些样本对于失效概率计算来说是非常重要的;拉丁超立方抽取样本数量越多，分层就越多，极端的尾部值就更容易被抽取。图 1拉丁超立方抽样原理Fig 1The Latin hypercube sampling princple拉丁超立方抽样假定所有随机变量是相互独立的，而在岩土工程中很多参数往往具有相关性，因此，无法直接进行拉丁超立方抽样，需要进行相应的变换。具体变换方式可参考相关随机变量蒙特卡洛抽样。常见 4 种不同类型相关随机变量蒙特卡洛抽样的实现步骤如下:1)对于均为正态分布的多维相关随机变量，可基于协方差矩阵的 Choleshy 因子分解将其转换为独立标准正态分布的随机变量进行蒙特卡洛抽样11。2)对于均为对数正态分布的多维相关随机变量，可将其转换为正态分布的多维相关随机变量进行蒙特卡洛抽样，转换后协方差矩阵跟随改变12。3)对于非正态分布的多维相关随机变量，可采用 Nataf 变换分 2 步将其转换为独立标准正态分布的随机变量进行蒙特卡洛抽样13。4)对于已知秩相关矩阵的多维相关随机变量，可将其转换为具有相同秩相关矩阵及相同维度的标准正态分布的相关随机变量，再根据相关随机变量的边缘分布生成所需要的样本，使得每个维度上相关随机变量样本从大到小的排序与对应维度上标准正态分布相关随机变量蒙特卡洛抽样样本从大到小的排序一致，从而实现秩相关随机变量的蒙特卡洛抽样14。笔者参考相关随机变量蒙特卡洛抽样的实现，结合拉丁超立方抽样原理，采用 Python 语言编程实现了上述 4 种类型的相关随机变量拉丁超立方抽样。1相关随机变量拉丁超立方抽样原理及实现11相关随机变量 Pearson 相关111相关随机变量的正态分布正态分布进行线性变化后仍然是正态分布，且相关性不变，因此，任意满足正态分布的相关随机变量均可以通过标准正态分布的线性组合得到，表达式如式(1):x=A+(1)式中:x 为一组正态分布的相关随机变量;为各相关随机变量的均值矩阵;为一组相互独立的标准正态分布变量;A 为转换矩阵。G ZEOVNIK 等12 证明了 V=AAT，其中 V 为相关随机变量 x 的协方差系数矩阵。又知 V 是一个对称的正定矩阵，可以通过 Choleshy 分解成 2 个矩阵的乘积，因此，只要得到 V 的 Choleshy 分解矩阵就可以求得矩阵 A。在 Python 语言中，使用 numpy 库实现 Choleshy分解，得到 A=nplinalgcholesky(V)。如果是给出的相关系数矩阵，则 V=T，为对角阵，对角线上的元素为各相关随机变量的标准差。所以，对x 的拉丁超立方抽样可以转化为对独立标准正态分布 的拉丁超立方抽样。实现拉丁超立方抽样用到Python 的 scipystats 中的 qmc 和 norm。112相关随机变量的对数正态分布对数正态分布拉丁超立方抽样与正态分布拉丁超立方抽样具有相似性，若 xi服从对数正态分布，则 lnxi服从正态分布，即lnxi=nj=1Aijj+ui(2)式中:xi为第 i 个相关随机变量;Aij为转换系数;j为第 j 个独立标准正态分布变量;i为 lnxi的均值。令 yi=lnxi(3)35第 3 期章从旭，等:相关随机变量拉丁超立方抽样及隧道掌子面稳定性概率分析则，式(2)可以写成:y=A+(4)可见，式(4)的形式同式(1)。因此，一组满足对数正态分布的相关随机变量的拉丁超立方抽样，可以转换为一组满足正态分布相关随机变量的拉丁超立方抽样。根据 G ZEOVNIK 等12 研究，原对数正态分布 x 与转换后正态分布 y 之间的协方差系数的关系如式(5):lnVijxixj+1()=Vij?(5)式中:Vij为相关随机变量 xi与 xj之间的协方差;Vij?为相关随机变量 yi与 yj之间的协方差;xi、xj为相关随机变量 xi与 xj的均值。根据式(6)可以求得转换后 yi的均值:E(yi)=ln E(xi)12ln 1+D(xi)E(xi)2(6)式中:E(yi)为第 i 个 y 变量的均值;E(xi)为第 i 个x 变量的均值;D(xi)为第 i 个 x 变量的方差。得到了正态分布 y 的 A 与，后续进行拉丁超立方抽样的计算步骤与正态分布的拉丁超立方抽样一致。得到相关随机变量 yi的样本后，通过 xi=eyi转换，即得到相关随机变量 x 的样本。113相关随机变量的非正态分布LIU Peiling 等15-16 提出的 Nataf 转换法，将非正态分布转换到标准正态分布空间下，只需要知道变量的边缘分布以及变量间的相关系数，即可以把原始分布空间转换到独立标准正态空间。Nataf 转换分 2 步:将原始空间转换成相关的标准正态空间;将相关的标准正态空间转换成独立的标准正态的空间，这可通过正态分布相关随机变量的拉丁超立方抽样来解决。原始空间向相关的标准正态空间转换的操作为:输入一组相关随机变量 X=(x1，x2，xn)，若相关随机变量 xi(i=1，2，n)的概率密度函数fi(xi)和累积密度函数 Fi(xi)已知，通过等概率转换原则15-16 可得:xifi(xi)dxi=Fi(xi)=(yi)yi=1 Fi(xi)(7)式中:yi为一组标准正态相关随机变量 Y=(y1，y2，yn)中的变量;()、1()分别为标准正态累积密度分布函数和逆累积密度分布函数。设 X 中任意 2 个相关随机变量 xi、xj的相关系数为 xixj，则xixj=E(xixi)(xjxj)xixj=F1i(yi)xixiF1j(yj)xjxj121 2yiyj exp(y2i 2 yiyjyiyj+y2j)2(1 2yiyj)dyidyj，(i，j=1，2，n，i j)(8)式中:xi、xj分别为 xi、xj的均值;xi、xj分别为 xi、xj的标准差;yiyj为 yi、yj的相关系数。在二元标准正态分布的联合概率密度函数数值计算中，式(8)的双重积分上下限选择4，4 就足够精确了13。在 xixj、xi、xj已知的情况下，几乎不可能直接计算 yiyj，因此，需采用数值计算法求解yiyj。当 xixj0 时，具体步骤如下:Step 1Step 1在 0，1 区间假设一个 yiyj值。Step 2Step 2将假设的 yiyj值代入式(8)，计算得到xixj。Step 3Step 3若计算得到的 xixj与实际 xixj的误差不满足精度要求，则使用二分法调整 yiyj值，返回Step 2Step 2，直到计算得到的 xixj与实际 xixj的误差满足精度要求，此时 xixj对应的 yiyj即为所需值。当 xixj0 时，在 1，0 区间进行类似操作。114相关随机变量拉丁超立方抽样举例以岩土工程中土体常见的黏聚力 c、内摩擦角 为例进行拉丁超立方抽样，抽取 10 000 个样本。均值 c=10 kPa，=30，标准差 c=2 kPa，=2，Pearson 相关系数为05。分别进行正态分布相关随机变量抽样(c、均满足正态分布)、对数正态分布相关随机变量抽样(c、均满足对数正态分布)、非正态分布相关随机变量抽样(c 满足对数正态分布，满足正态分布)，抽样结果分布如图 2，抽样参数 c、概率密度如图 3。45重 庆 交 通 大 学 学 报(自 然 科 学 版)第 42 卷(a)正态分布抽样(b)对数正态分布抽样(c)非正态分布抽样图 2正态分布、对数正态分布、非正态分布抽样结果Fig 2Sampling result of normal distribution，lognormal normal distribution and nonnormal distribution图 3正态分布、对数正态分布、非正态分布抽样概率密度Fig 3Sampling of probability density of normal distribution，lognormal normal distribution and nonnormal distribution55第 3 期章从旭，等:相关随机变量拉丁超立方抽样及隧道掌子面稳定性概率分析由图 3 可见，当 Pearson 相关系数不等于 0 时，c 与 的抽样概率密度分布曲线与理论概率密度曲线存在些许差异。分析原因，一方面可能是拉丁超立方抽样样本数量相对较少造成的;另一方面，可能是每个区间抽样的随机性造成的，尤其是在均值附近。综上，随机变量 Pearson 相关的 3 种不同类型的拉丁超立方抽样保持了原变量的分布，在保持抽样精度的同时，大大降低了样本数量。12随机变量 Spearman 秩相关Spearman 秩相关系数只对变量的秩次大小进行线性相关分析，所以没有要求变量的分布，但秩相关系数与样本的排序有关，而与样本具体的值关联不大，因此，秩相关随机变量的拉丁超立方抽样可通过秩相关矩阵得到。相关随机变量 Z 如果满足 Z=PT，而 VS=PPT，VS为变量 Z 秩相关矩阵，则 Z 与 PT具有相同的秩相关系数，其中 为一组相互独立的随机变量14。既然秩相关系数与相关随机变量具体的分布没有关系，那么可以设 为标准正态分布，即 Z 为正态分布。假设一组包含 m 个秩相关的相关随机变量X=(x1，x2，xm)需抽取 N 个样本，秩相关矩阵为VS，通过 Python 语言实现拉丁超立方抽样过程如下:Step 1Step 1通过拉丁超立方抽样对随机变量 进行抽样，=(1，2，m)，其中:1、2、m为相互独立的标准正态分布随机变量，每个变量共抽取N 个样本。Step 2Step 2对 VS进行 Choleshy 分解得到 PT，通过Z=PT转换得到正态分布 Z 的 N 个样本，分别对Z=(z1，z2，zm)中每个变量的 N 个样本进行排序，得到排序序号。Step 3Step 3通过拉丁超立方抽样分别对 xi(i=1，2，m)抽取 N 个样本，并对这 N 个样本进行排序，排序序号与第 i 个变量 zi的 N 个样本排序号相同，从而得到秩相关随机变量 X=(x1，x2，xm)的 N 个拉丁超立方抽样样本。以岩土工程中土体常见的参数 c、为例进行拉丁超立方抽样，抽取 10 000 个样本。两个参数的均值 u 分别为 uc=10 kPa、u=30，标准差 c=2 kPa、=2，其中 c 满足对数正态分布、满足正态分布，Spearman 秩相关系数为05。抽样结果如图 4。图 4Spearman 秩相关抽样结果Fig 4Sampling result of spearman rank correlation2隧道掌子面失效概率分析21实例背景目前，隧道工程仍以确定性分析为主，但其围岩的材料参数往往具有不确定性，因此，有必要对某些情况进行不确定性分析或可靠度设计17。不确定性分析方法可以是蒙特卡洛抽样和数值计算方法的结合18。为了展示拉丁超立方抽样在工程应用中的优越性，笔者用一个直径 D=6 m，埋深 H=6 m 的圆柱形隧道掌子面进行失效概率分析。隧道周围土体采用摩尔-库伦模型，衬砌采用弹性模型。将土体的黏聚力 c 和内摩擦角 考虑为服从正态分布的随机变量，两个参数的均值 u 分别为uc=10 kPa、u=30，标准差 c=2 kPa、=2，Pearson 相关系数为05。土体的重度=18 kN/m3，弹性模量 E=240 MPa，泊松比=03。隧道采用全断面开挖，不考虑超前支护措施。22三维数值模型图 5 为概率分析所用的三维数值模型，考虑隧道结构和荷载的对称性，仅建立了 1/2 模型进行分析。为避免边界效应的影响，根据圣维南原理，隧道开挖边界距模型边界应为 35 倍洞径 D，因此分析时，隧道开挖边界至模型的水平边界以及模型底部边界的距离均取 3D，掌子面距模型纵向边界的距离也为 3D。在模型的纵向上，掌子面前方一倍洞径的土体网格单元为 05 m，其它部分的网格纵向尺寸为2 m;模型单元为 57 400 个，节点为 62 118 个;边界条件上，模型的左右和前后边界为法向约束，模型的底部为全约束。65重 庆 交 通 大 学 学 报(自 然 科 学 版)第 42 卷图 5数值模型Fig 5Numerical model23模拟过程要计算图 5 数值模型的隧道掌子面的失效概率，需抽取一定数量的样本进行确定性分析。为保证一定的精度，样本的数量 N 应满足N100/Pf(10)式中:Pf为失效概率。笔者采用批量采样的方式，进行确定性分析，统计失效样本的数量，每批次采集1 000个样本;重复这一步骤，直到满足 NPf100，终止采样，计算最终的失效概率。对采集的每个样本(ci，i)，传统方法是通过FLAC3D自带的强度折减法来计算安全系数 f 的。f1，表示隧道掌子面稳定;f1，表示隧道掌子面不稳定。为提高计算效率，笔者不计算安全系数具体值，而是通过运行命令“solve fos bracket 10 10”得到计算的最大不平衡力比值 r。当 r1105时，表明迭代计算是收敛的，即在不进行折减的情况下，隧道掌子面是稳定的;当 r1105时，掌子面不稳定。图 6 为数值模型在某一黏聚力 c 与内摩擦角 组合下的计算结果。可见，模型基本没有塑性区，r=996 106，因此，掌子面是稳定的。图 6三维数值模型计算结果Fig 6Calculation results of three-dimensional numerical model24结果分析图 7 为用蒙特卡洛法计算得到的隧道掌子面失效概率分析结果。图 7样本数量与失效概率Fig 7Number of samples and failure probability由图 7 可见，当样本数量 N=97 000 个时，掌子面失效样本的数量为 100 个，此时的失效概率 Pf=100/97 0000001。在蒙特卡洛模拟中，95%置信区间允许的误差为19:=2001PfNPf=20%(11)即，满足 95%置信区间的失效概率范围为:08Pf，12Pf=0000 82，0001 24。根据拉丁超立方抽样采样结果，当样本数量达到 16 000 次时，已能满足所需的计算精度。所需样本数量仅为蒙特卡洛抽样样本数的 165%。由此可见，拉丁超立方抽样更具代表性，可以显著降低概率分析所用的样本数量。3结论拉丁超立方抽样从多元分布中独立抽取最具代表性的样本，在保证计算精度的同时，可大大减少抽样样本数，但拉丁超立方抽样不适用于相关随机变量的直接抽样，需进行一定的转换。笔者通过Python 编程将 4 种不同类型的相关随机变量转换为独立标准正态分布随机变量进行拉丁超立方抽样。具体为:正态分布与对数正态分布用协方差矩阵通过 Choleshy 分解实现相关随机变量的转换;非正态分布通过 Nataf 变换实现相关随机变量的转换;秩相关通过秩相关矩阵分解实现相关随机变量的转换。同时，基于独立标准正态分布相关随机变量拉丁超立方抽样，开展了隧道掌子面稳定性概率分析，并与蒙特卡洛模拟结果进行了对比。研究得到以下主要结论:1)相关随机变量拉丁超立方抽样的抽样结果与75第 3 期章从旭，等:相关随机变量拉丁超立方抽样及隧道掌子面稳定性概率分析其理论分布结果基本吻合，验证了拉丁超立方抽样方法的准确性。2)隧道掌子面稳定性概率分析表明，拉丁超立方抽样所需样本数为蒙特卡洛抽样样本数的 165%即可满足计算精度要求，极大地提高了计算效率。参考文献(eferences):1仉文岗，王琦，刘汉龙，等岩体空间变异性对隧道拱顶失效概率的影响J 岩土力学，2021，42(5):1462-1472ZHANG Wengang，WANG Qi，LIU Hanlong，et al Influence ofrock mass spatial variability on probability of tunnel roof wedgefailure Jock and Soil Mechanics，2021，42(5):1462-1472 2仉文岗，王琦，陈福勇，等考虑岩体空间变异性的边坡可靠度分析及抗滑桩随机响应研究J 岩土力学，2021，42(11):3157-3268ZHANG Wengang，WANG Qi，CHEN Fuyong，et al eliabilityanalysis of slope and random response of anti-sliding pile consideringspatial variability of rock mass properties J ock and SoilMechanics，2021，42(11):3157-3268 3 LI Bin，FU Yong，HONG Yi，et al Deterministic and probabilisticanaly-sis of tunnel face stability using support vector machine J Geomechanics and Engineering，2021，25(1):17-30 4ZENG Peng，ZHANG Tianlong，LI Tianbin，et al Binary classificationmethod for efficient and accurate system reliability analyses of layeredsoil slopes J Georisk:Assessment and Management of isk forEngineered Systems and Geohazards，2020，16(3):435-451 5陈将宏，李建林，许晓亮，等 相关变量生成算法及边坡可靠度Monte Carlo 模拟 J 岩土力学，2017，38(11):3341-3346CHEN Jianghong，LI Jianlin，XU Xiaoliang，et al Algorithm forgeneration correlative variables and Monte Carlo simulation of slopereliability J ock and Soil Mechanics，2017，38(11):3341-3346 6吴帅兵，张坤，李典庆 相关非正态变量变换时相关性变化对可靠度的影响J 武汉大学学报(工学版)，2011，44(2):151-155WU Shuaibing，ZHANG Kun，LI Dianqing Effect of correlationchange from transformation of correlated abnormal variables onstructural reliabilityJEngineering Journal of Wuhan University，2011，44(2):151-155 7邢万波，周钟，巩满福 Monte Carlo 法和几何法的讨论和改进J 岩土力学，2007，28(增刊 1):211-216XINGWanbo，ZHOUZhong，GONGManfuDiscussionandimprovement on Monte Carlo method and geometric method J ockand Soil Mechanics，2007，28(Sup 1):211-216 8HUANG Jinsong，FENTON G，GIFFITHS D V，et al On the efficientestimation of small failure probability in slopes J Landslides，2017，14:491-498 9LI Tianzheng，DIAS D Tunnel face reliability analysis using activelearning Kriging modelCase of a two-layer soils J Journal ofCentral South University，2019，26(7):1735-1746 10MCKAY M D，BECKMAN J，CONOVE W J A comparison ofthree methods for selecting values of input variables in the analysisof output from a computer code J Technometrics，1979，21(2):239-24511 CAI G Q Generation of correlated random variables and stochasticprocesses J Probabilistic Engineering Mechanics，2018，52:40-46 12ZEOVNIK G，TKOV A，KODELI I A Correlated random samplingfor multivariate normal and log-normal distributions J NuclearInstruments Methods in Physics esearch Section A:Accelerators，Spectrometers，Detectors and Associated Equipment，2012，690:75-78 13文德智，卓仁鸿，丁大杰，等 蒙特卡罗模拟中相关变量随机数序列的产生方法J 物理学报，2012，61(22):220204WEN Dezhi，ZHUO enhong，DING Dajie，et al Generationrandomnumber sequence of related variables in Monte Carlo simulation J Acta Physica Sinica，2012，61(22):220204 14IMAN L，CONOVE W J A distribution-free approach to inducingrank correlation among input variables J Communications inStatistics-Simulation and Computation，1982，11(3):311-334 15Der KIUEGHIAN A，LIU Peiling Structural reliability under incom-plete probability information J Journal of Engineering Mechanics，1986，112(1):85-104 16LIU Peiling，Der KIUEGHIAN A Multivariate distribution modelswith prescribed marginals and covariances J Probabilisticengineering mechanics，1986，1(2):105-112 17仉文岗，王焱，宗梓煦，等 隧道施工与运营过程中的可靠度分析方法研究 J 重庆交通大学学报(自然科学版)，2020，39(3):1-13ZHANG Wengang，WANG Yan，ZONG Zixu，et al eliabilityanalysis methods in tunnel construction and operation JJournalof Chongqing Jiaotong University(Natural Science)，2020，39(3):1-13 18郑颖人 有限元极限分析法在隧洞工程中的应用 J 重庆交通大学学报(自然科学版)，2011，30(增刊 2):1127-1137ZHENG Yingren The application of FEM limit analysis in tunnelengineering J Journal of Chongqing Jiaotong University(NaturalScience)，2011，30(Sup 2):1127-1137 19LOW B K，PHOON K K eliability-based design and its comple-mentary role to Eurocode 7 design approach J Computers andGeotechnics，2015，65:30-44(责任编辑:田文玉)85重 庆 交 通 大 学 学 报(自 然 科 学 版)第 42 卷