考虑地层空间变异性的岩溶区基础可靠度分析

吴高桥 聂建国

摘要：基于隨机场理论和Karhunen-Loeve展开构造随机场以模拟具有空间变异性的地层，并结合有限元极限分析法和蒙特卡洛模拟对含溶洞地层上的基础开展可靠度分析.本文重点是通过大量数值模拟进行全面参数分析，探究地层中溶洞特征（如溶洞埋深、溶洞水平偏移距离）以及土体空间变异性参数（如土体相关长度、土体变异系数）对基础极限承载力以及失效概率的影响；进而结合可视化溶洞-基础体系破坏模式深入揭示岩溶区基础失效概率变化机理；最后将有限元极限分析结果与已有研究进行对比，验证了数值模型的准确性及研究结论的可靠性.结果表明体系破坏线发散度越大，则基础在受荷破坏时可能发展出的破坏路径越多，从而导致失效概率提高；此外，下伏溶洞距基础越远或土体空间变异性越强则基础失效概率越高，且失效概率将随所选取安全系数的增加而大幅降低，在前期设计中若采用本文方法进行可靠度分析，为确保安全性，安全系数应尽量大于2并采用下限理论进行分析.

关键词：基础工程；溶洞；有限元极限分析；随机场理论；随机分析

中图分类号：U416.161文献标志码：A

Reliability Analysis of Footings Lying on Karst Area Considering Spatially Variability of Stratum

WU Gaoqiao，NIE Jianguo

（Department of Civil Engineering，Tsinghua University，Beijing 100084，China）

Abstract：This study simulates the soils with spatial variability by random field theory and Karhunen-Loeve expansion. Reliability analysis of strip footing lying on karst area is carried out based on finite element limit analysis and Monte Carlo simulation. The emphasis of this study is performing the parametric analysis through finite element modeling，in order to explore the effect of karst cave parameters （including depth and horizontal offset distance of karst caves）and spatial variability parameters（including correlation length and coefficient of variation）on the ultimate bearing capacity and failure probability of footings. Then，combined with the visualization of failure patterns of cave-footing system，the variation mechanism of footing failure probability is revealed. Finally，the finite element analysis results are compared with the existing research results to verify the accuracy of the numerical model and the reliability of the research conclusions. The results indicate that the higher the disperse degree of failure curve is，the more failure paths may develop under external loads，which leads to a higher failure probability of the footing. Fur- thermore，the footing failure probability increases with the increase of distance between the karst cave and the footing. The higher spatially variability of soils also leads to a higher failure probability of the footing. And the failure probability decreases significantly with the increasing factor of safety. If the method proposed in this paper is used for reliability analysis in the advanced design，the safety factor should be larger than 2 as far as possible and the lower limit theory should be used for analysis，in order to guarantee the safety.

Key words：foundation engineering；Karst cave；finite element limit analysis；random field theory；stochastic analysis

我国喀斯特地貌占地约344万平方千米，广泛分布于我国西部地区.随着近年城市化进程加速，土地紧缺问题日益严重，因此越来越多建筑将无可避免地选址在岩溶发育地区.然而近年来，因地层中空洞失稳引起的上部构筑物沉降和倾斜等灾害时有发生，据2020年中央级媒体、部委网站、公开出版物、中央重点新闻网站及地方重点报网等报道的数据，2020年地下空洞相关灾害与事故共237起，数量与2019年持平[1]，造成巨大的人员伤亡和经济损失.因此，对含溶洞地层上基础的可靠性进行分析，将有利于上部构筑物的安全性，并为上部施工提供指导，具备重要理论意义和实际工程价值.

国内外有众多学者通过理论分析、模型试验和数值模拟等手段针对空洞-基础体系稳定性开展研究.在理论分析方面，Wang和Hsieh[2]基于传统上限分析法研究了圆形空洞正上方条形基础的稳定性问题，提出基础稳定性理论计算框架和基础极限承载力计算公式，并总结3种典型破坏模式.随后Wang等[3]继续研究了空洞与条形基础水平距离对体系破坏模式的影响，将典型破坏模式扩展至10种，并归纳了不同破坏模式之间的演化规律.刘辉等[4]通过极限分析上限法合理构造速度场，进而通过剪切耗散和基础功率的关系推出溶洞上方条形基础的极限承载力计算方法，并详细分析了不同参数对体系破坏模式的影响.赵明华等[5]考虑自然溶洞形成过程中溶洞顶板的空间形态特征，采用结构力学分析理论建立了溶洞顶板最小抗弯厚度的计算公式，并揭示了该厚度随空洞跨度、地层物理力学性质以及基础所受荷载改变而变化的规律.在模型试验方面，Wood等[6]和Al-Tabbaa等[7]基于模型试验分析了下伏空洞对上部浅基础承载力的削弱作用.Kiyosumi等[8]研究了溶洞-基础体系的稳定性问题，通过对基础失去承载能力后溶洞周边岩体的破坏情况，总结了多溶洞工况下，不同排布位置的溶洞对上部基础极限承载力的影响.刘铁雄等[9]和王革立[10]将溶洞上方地层视为梁板，通过一系列试验分析了溶洞跨径、溶洞上覆岩土體厚度等对基础稳定性的影响.赵明华等[11]通过模型试验，研究了上部基础荷载作用下的岩溶顶板的承载机理，对诸如顶板厚跨比、地层力学性能等参数开展了参数分析，并总结了溶洞安全埋深的计算公式.

上述研究使研究工程人员对空洞-基础稳定性问题有了更深入的认知，但使用理论分析时，通常要事先对问题进行一定假设以降低求解难度，因此仅能分析相对简单的问题.而模型试验受限于时间成本和经济成本，无法开展大规模研究以进行参数分析.随着计算机算力的发展，数值模拟成为了一种被广泛使用的岩土工程分析方法.Baus和Wang[12]在1983年采用有限元法（finite element method，FEM）研究了粉质黏土中空洞上方条形基础的稳定性问题，并总结出能够反映空洞尺寸、空洞位置及条形基础埋深对基础极限承载力影响的设计图.Badie和Wang[13]采用三维有限元程序对基础-空洞体系稳定性进行分析，结果表明基础下方存在某临界区域，只有当下伏空洞处于该区域范围内时空洞才会对基础造成显著影响，且该临界区域的范围取决于土体物理力学性质以及空洞大小和位置.Jao等[14]使用有限元法研究了条形基础荷载作用下的有衬砌隧道稳定性问题，对不同地层类型、隧道大小、隧道位置以及衬砌厚度的影响进行了详细探讨.孙映霞等[15]采用有限差分法（finite difference method，FDM）研究下伏空洞埋深以及距基础距离对体系失稳机理的影响. 随后在其基础上，Kiyosumi等[16]使用有限元软件PLAXIS率先对多空洞上方条形基础的极限承载力进行研究，分析了不同位置空洞对基础造成影响强弱的关系，结果表明靠近基础的空洞对基础承载力造成的影响远大于较远处空洞.基于Kiyosumi的研究，Lavasan等[17]分析了位于两个平行圆形空洞上方的条形基础的极限承载力，深入研究了空洞的破坏机理并归纳了数种典型破坏模式.随后Zhou等[18]和Lee等[19]分别采用局部不连续优化法（discontinuity layout optimization，DLO）和有限元法分析了下伏双空洞对条形基础稳定性的影响.

上述研究基本全部假设空洞所处地层为均匀材质，然而现实中地层往往具有一定的空间变异性，这将导致基础承载机理发生改变，进而增加基础承载力的不确定性，因此有必要在基础稳定性分析中充分考虑土体空间变异性的影响.鉴于此，本研究结合随机场理论和Karhunen-Loeve展开构造实现随机场土体强度的随机分布，进而通过有限元极限分析和蒙特卡洛模拟（Monte Carlo simulation）对含溶洞地层上基础进行可靠度分析，着重分析下伏溶洞位置及土体空间变异程度对基础稳定性及失效概率的影响，以期为工程设计提供参考，并进一步确保岩溶区基础的安全性.

1问题描述

1.1基本假定

图1为待研究问题的平面模型示意图，问题关键影响因素已标示于图中.具体为：地表存在宽度为D的条形基础，在其中心位置施加有竖直集中荷载q_u.基础下方地层存在边长为D的正方形溶洞，溶洞中心到地表及基础中垂线的垂直距离分别为H和S.空洞周围的土体遵循特雷斯卡屈服准则（Tresca yield criterion）且存在一定空间变异性，其重度γ=20 kN/m³，平均剪切强度s_u=40 kPa.

1.2数值模型建立

为了尽量规避边界效应，需要建立尺寸足够大的数值模型，以保证在所有工况中塑性破坏区不会延伸至模型的侧面边界以及下边界.经测试，当模型尺寸宽度为30D、高度为15D时既可满足上述要求，也能保证较高的运算效率.此外，基于已有研究[19-21]，本文在模型底部边界设置完全约束，左右两侧边界设置水平约束，在地表边界则不设置约束.由于本文拟采用蒙特卡洛法分析基础失效概率，所以对每组工况均需开展上千次模拟以获得准确结果，这将耗费大量时间.因此为了尽可能提高分析效率和精度，采用网格自适应迭代技术基于分析域内能量耗散情况对网格进行重分布优化.针对本文数值模型，将网格初始数量和最终数量分别设置为1 000 和10 000，自适应迭代次数设为3次，经过重分布优化后的效果图如图2所示.

1.3随机场实现

本文假定地层土体遵循Tresca屈服准则，其中对基础稳定性影响最大的参数为不排水剪切强度s_u，因此在本文的随机分析中主要考虑土体剪切强度的空间变异性.为实现土体空间变异性，将土体不排水剪切强度在分析域内模拟为对数正态分布，其具体表达式为：

剪切强度的均值和标准差可以通过一种能简洁描述随机场土壤固有变异性的无量纲变异系数（coefficient of variation，COV）表示.该参数能够帮助设计者对土体中固有变异性的可能范围进行评估[22]，具体表达式为：

上述表达式中的对数正态分布参数可表示

如下：

对式（3）和式（4）进行转换可得：

则对数正态分布的累积分布函数（cumulativedistribution function，CDF）可表示为：

式中：erfc为互补误差函数.需要注意的是，CDF是随机分析中概率密度函数（probability density function，PDF）的积分形式.

除了上述对数正态分布参数，还有两个重要参数在随机场理论中至关重要，即用于描述空间变异性土体不同方向波动程度的水平相关长度（CL_x）以及垂直相关长度（CL_y）.相关长度被定义为土体中一临界距离，若土体中任意两点间距离小于该距离，则两点的剪切强度趋于密切相关，且相关长度越大则随机场中剪切强度的过渡越平滑；相反，相关长度越小则土体中各点剪切强度的关联性越弱，随机场趋于参差.由其定义可知当相关长度趋于无穷大时，分析域可视为均匀土体；而当其趋于0时，随机场中各点之间无任何关联.为了更直观地反映相关长度与随机场平滑度的关系，将水平、竖直相关长度不同的土体剪切强度分布模拟情况列于图3.对比图3（a）和图3（b）可以看出，当竖直相关长度较大时，土体剪切强度在竖直方向过渡较均匀；而当水平相关长度较大时，则土体剪切强度在水平方向过渡较均匀.由图3（c）可知若各向相关长度均较小，土体的剪切强度分布较离散且无规律.Phoon和Kulhawy[22]经过大量调研总结出土体在水平方向的随机性通常较弱，CL_x范围在46 m到60 m之间；而土体在竖直方向的随机性通常较强，CL_y的范围在0.8 m到6.1m之间.基于上述结论，本文中的CL_x固定为50 m，仅考虑CL_y，对基础稳定性的影响.为了方便表述，对相关长度进行无量纲化处理如下：

Θ_x=CL_x/D（8）

Θ=CL_y/D（9）

本研究选择通过能够提供指数协方差函数解析解的Karhunen-Loeve展开式生成随机场.由于该展开式的项数将随着随机场相关长度的减小而大幅增加[23]，因此为了保证随机分析结果的可靠度，将展开项数设置为1 000.此外，本研究采用蒙特卡洛模拟分析基础失效概率，若模拟次数不足，将严重影响所得结果的准确性，因此在后续随机分析中各工况模拟次數均设为500次.

2确定分析

通过随机分析研究基础失效概率问题的原理是先不考虑土体空间变异性，对各工况进行传统确定分析（deterministic analysis）；随后再通过多次模拟随机场并总结小于经安全系数换算过的确定性分析基础承载力的次数，进而得出基础失效概率.则失效概率的表达式为：

p_f=p（N_ram≤N_det/FS）（10）

式中：p_f为基础失效概率，在本文中具体定义为500次蒙特卡洛模拟所得极限承载力中小于确定性承载力N_det与安全系数比值的比例；FS为设计中选取的安全系数（factor of safety），一般来说安全系数越大，则设计越保守；N_ram为每次蒙特卡洛模拟所得基础极限承载力因子；N_det为确定分析中的基础极限承载力因子，其具体定义为：

式中：q_u为基础极限承载力；s_u为土体剪切强度.

基础极限承载力与下伏空洞所处位置密切相关，参照已有研究，在本节中对H/D=1、2、3、4以及S/D=0、0.5、1、2、3、4对应的工况进行确定分析，所得结果为上限解及下限解的平均值，基础极限承载力因子总结于图4以及表1.由图4可知基础极限承载力将随着S/D以及H/D的增加而提升，这是因为下伏溶洞离基础越远，则其对基础造成的负面影响越弱. 同时，H/D越大则对应的承载力曲线越平缓，这说明溶洞埋深越大则水平位置的改变对基础的影响越弱.此外，图4中有多个工况基础极限承载力因子集中在5.20附近，该现象说明此时溶洞由于距离过远不再对基础产生影响.

3随机分析

本节将基于确定性分析结果开展随机分析，以期全面探究空洞埋深、空洞水平偏移距离以及土体竖直相关长度对基础平均承载力因子以及失效概率的影响.其中基础平均承载力因子为每个工况500 次蒙特卡洛模拟所得极限承载力因子的平均值.由于随机分析中每个工况均需构造500个不同的随机场并计算相应的基础极限承载力，导致运算量巨大，因此对所有影响因素进行排列组合再开展全面分析是不现实的.鉴于此，本文参考正交试验设计原理，在探究某些参数的影响时固定其他参数，以达到减少运算量且不影响所得规律的目的.本文随机分析所考虑工况的具体参数见表2.

3.1空洞埋深的影响

为了探究空洞埋深对基础稳定性的影响，对S/D=0，Θ_y=1以及H/D=1，2，3，4对应工况开展蒙特卡洛模拟，所得基础平均极限承载力因子和失效概率分别总结于图5和图6.由图5可以看出随机分析所得的平均承载力因子变化趋势基本和确定分析一致，但由于土体存在空间变异性，导致蒙特卡洛模拟所得的大部分基础极限承载力因子小于确定分析，所以平均承载力因子会小于确定分析所得结果.

图6为各工况通过上、下限分析所得的基础失效概率.可以看到对几乎所有工况而言，基础失效概率都随空洞埋深的增加而增加，该结论似乎与实际情况不符，因为通常认为下伏溶洞距离基础越远则其对基础承载能力造成的负面影响越小.实际上确如图5所示，基础的极限承载力因子随着H/D的增加而单调增加.但本文中的失效概率是指采用确定分析所得基础极限承载力进行设计时，由土体空间变异性引起失稳的概率，与基础极限承载力的绝对大小无必然联系.值得一提的是，在实际工程设计中一般会通过增加安全系数的方式以保证安全性，因此有必要在失效概率分析中考虑安全系数的影响.由图6可以看出随着安全系数的增加，基础失效概率显著下降，当安全系数为2时基础的失效概率将趋于0，同时由曲线上升幅度可知安全系数越小则空洞埋深对基础失效概率的影响越显著.此外，通过对比上、下限分析对应的失效概率，可以看出使用上限分析进行蒙特卡洛模拟时所得基础失效概率较低.为了进一步探究基础失效概率变化机理，将某次蒙特卡洛模拟所得溶洞-基础体系破坏模式总结于图7，可以看出当空洞埋深较浅时，破坏线较集中，说明体系的破坏模式较单一.而随着溶洞埋深的增加，主破坏线逐渐向周围发散且发散程度越来越明显，而这些较为发散的破坏线可能导致土体中失稳路径（failure path）增加，进而引起破坏模式的复杂化并导致失效概率的增加.

3.2空洞水平偏移距离的影响

为了探究空洞水平偏移距离的影响，将H/D=2，Θ_y=1以及S/D=0.5，1，2，3，4工况对应的平均承载力因子和基础失效概率总结于图8和图9.由图8可以看出随机分析所得的平均承载力因子变化趋势与确定分析基本一致，这与图5所展示的规律基本一致. 由此可以认为只要蒙特卡洛模拟次数达到一定数量，那么土体空间变异性几乎不会对随机分析中的平均承载力因子变化趋势造成影响.

从图9可看出除FS=1的下限分析以及FS=1.5的上限分析在S/D=0时失效概率略高于S/D=0.5，其余工况的基础失效概率随溶洞水平偏移距离的增加而增加.图10为各工况对应的体系破坏模式，可以看出当溶洞水平偏移距离较小时，破坏主要集中在洞顶，且破坏线较细.随着空洞偏移距离的持续增加，破坏模式逐渐变为洞顶-洞侧联合破坏，且溶洞顶部将受到来自其正上方土体的挤压.此外可以看到S/D=0和S/D=0.5时破坏线离散度较低，因此其所对应的失效概率接近.而随着溶洞水平偏移距离的增加，破坏线发散程度加剧进而导致基础受荷时土体中可能发展出的失稳路径增多，引起失效概率的增加，这与图9所展示的规律相吻合.结合3.1与3.2节的研究，可以认为对大部分工况而言，除了溶洞水平偏移距离极小的工况，下伏溶洞与基础的间距（包括埋深及水平偏移距离）越大，则对应的基础失效概率越高.此外结合图6和图9可知采用极限分析下限法进行随机分析时，所得基础失效概率相对于上限法明显偏高，而承载力偏低，亦即采用下限法进行施工设计时总体偏于保守.

3.3土体空间变异性参数的影响

为了探究土体空间变异性参数对基础稳定性的影响，在本节对H/D=2、S/D=1时不同Θ_y和COV對应的基础平均承载力因子和失效概率进行分析.由图11可知基础平均承载力因子随竖直相关长度的增加而单调递增，这是因为土体强度相关性越高则地层中剪切强度的过渡更平滑，从而使得由于土体强度分布不均导致的极低承载力工况数量减少.且承载力增强速度将逐渐平缓，当Θ_y=4时平均承载力因子已与确定性分析所得结果（图中水平虚线所示）相差无几，这说明此时土体空间变异性的影响较弱，已趋近于均质土.此外，COV越大则平均承载力越低，这是因为土体变异系数越大则土体空间变异性越强，这将导致蒙特卡洛模拟中低承载力工况数量的增多，进而引起平均承载力因子的下降.

图12为Θ_y和COV对基础失效概率的影响.可以看到失效概率随Θ_y的增加而单调递减，且当Θ_y很小时，失效概率相对很高，此时Θ_y只需略作提高则基础失效概率将大幅降低，而随着Θ_y的进一步增加，失效概率的降低速度减缓.结合上文所得结论，可知空间变异性较强时，空间变异性参数对基础承载力以及失效概率的影响较显著，而当空间变异性较弱时，空间变异性参数的改变仅对基础产生微弱影响.

图13为不同Θ_y和COV对应的破坏模式图.对比图13（a）（c）（e）（g）和图13（b）（d）（f）（h）可以看到当Θ_y较小时，破坏线的离散度很高，因此所对应的基础失效概率很高；而当Θ_y从2增加到4时，体系破坏模式及破坏线几乎没有变化，说明此时基础失效概率趋于稳定，与图12所示规律相吻合.而对比COV=25%（左侧纵列）和COV=60%（右侧纵列）对应的破坏模式，可以看出COV=60%对应的地层空间变异性强度明显高于COV=25%时.以图13（b）为例，在基础和溶洞间的主要破坏线之外还发育有多条破坏路径，这将导致基础失效概率剧增及承载力下降，因此COV越大则基础失效概率越高，平均承载力因子越低，与图11和图12所示规律一致.

4模型验证

为了验证有限元极限分析模型的可靠性，将本文确定分析的结果与已有有限元法"9]和不连续局部优化法结果118]进行对比，如图14所示.通过对比可以发现本文上限分析结果略高于已有研究结果，下限分析结果则略低于已有研究结果，但总体变化趋势基本吻合，且最大误差小于5%.此外，有限元极限分析的真值是介于上、下限分析结果之间的，因此可以认为本文数值模型所得结果是准确可靠的.

5结论

本文通过随机场理论和Karhunen-Loeve展开的方式构建随机场，对岩溶发育地层上基础开展了一系列随机分析.采用有限元极限分析法以及蒙特卡洛模拟全面探究了溶洞埋深、溶洞水平偏移距离、土体相关长度以及土体变异系数对基础稳定性以及失效概率的影响，并结合可视化体系破坏模式深入分析基础失效概率变化机理.本研究所得结论如下：

1）在随机分析中，若蒙特卡洛模拟样本数足够，则随机分析所得基础承载力变化趋势与确定分析所得趋势几乎一致；

2）基础的承载力将随溶洞埋深、溶洞偏移距离、土体相关长度的增大以及土体变异系数的减小而提升，且上述参数较小时对基础的影响更显著；

3）溶洞与基础的间距越远，则基础的失效概率越高.且基础的失效概率随土体相关长度的增大而降低、随土体变异系数的增大而升高；

4）在施工设计时采用极限分析下限法对基础稳定性进行随机分析时将偏于保守，所得基础失效概率偏高且极限承载力偏低，且基础失效概率将随所取安全系数的增加而大幅降低；

5）破坏模式中离散度较高的破坏线将导致基础受荷后可能发生破坏的失稳路径增多，从而导致更高的失效概率.

参考文献

[1]中国工程院战略咨询中心，中国岩石力学与工程学会地下空间分会，中国城市规划学会.2021中国城市地下空间发展蓝皮书[R].上海，2021：12.

Strategic Consulting Center of Chinese Academy of Engineering，Underground Space Branch of Chinese Society of Rock Mechanics，and Engineering，Urban Planning Society of China. 2021 Blue book of urban underground space development in China[R]. Shanghai，2021：12.（In Chinese）

[2]WANG M C，HSIEH C W. Collapse load of strip footing above circular void[J]. Journal of Geotechnical Engineering，1987，113 （5）：511-515.

[3]WANG M C，KIM Y U，JUN J T. Cavity effect on collapse load of strip footings [C]//Proceedings of the International Conference on Soil Mechanics and Geotechnical Engineering. Lisse，The Netherlands Exton，PA：AA Balkema Publishers，2001：317-320.

[4]刘辉，杨峰，阳军生.空洞上方浅基础地基破坏模式与极限承载力分析[J].岩土力学，2010，31（11）：3373-3378.

LIU H，YANG F，YANG J S. Failure mechanism and ultimate bearing capacity of strip footing above void[J]. Rock and Soil Mechanics，2010，31（11）：3373-3378.（In Chinese）

[5]赵明华，卢晓明，张锐.下伏空洞地基极限承载力与破坏模式上限有限元法[J].岩土力学，2017，38（1）：229-236.

ZHAO M H，LU X M，ZHANG R. Upper bound finite element method for ultimate bearing capacity and failure mechanism of subgrade above void [J]. Rock and Soil Mechanics，2017，38（1）：229-236.（In Chinese）

[6]WOOD L A，LARNACH W J. The behavior of footings located above voids[C]// Proc 11th Int Conf Soil Mech Found Eng. San Francisco：AA Balkema，1985：273-276.

[7]AL-TABBAA A，RUSSELL L，O，REILLY M. Model tests of footings above shallow cavities [J]. Ground Engineering，1989，22 （7）：39-42.

[8]KIYOSUMI M，KUSAKABE O，OHUCHI M. Model tests and analyses of bearing capacity of strip footing on stiff ground with voids[J]. Journal of Geotechnical and Geoenvironmental Engineering，2011，137（4）：363-375.

[9]刘铁雄，彭振斌，安伟刚，等.岩溶地区桩基特性物理模拟[J].中南工业大学学报（自然科学版），2002，33（4）：339-343.

LIU T X，PENG Z B，AN W G，et al. Simulated research of pile foundation character in the Karst area[J]. Journal of Central South University of Technology（Natural Science），2002，33（4）：339-343.（In Chinese）

[10]王革立.巖溶地基嵌岩桩桩基特性分析与试验研究[D].长沙：中南大学，2002.

WANG G L. The characteristic analysis and model test of socketed pile on karst area[D]. Changsha：Central South University，2002. （In Chinese）

[11]赵明华，吴高桥，肖尧，等.不同厚度溶洞顶板与基桩作用机理室内模型试验研究[J].水文地质工程地质，2017，44（6）：29-36.

ZHAO M H，WU G Q，XIAO Y，et al. Laboratory model tests on punching properties of karst cave roof under pile tip considering various thickness-span ratios[J]. Hydrogeology and Engineering Geology，2017，44（6）：29-36.（In Chinese）

[12]BAUS R L，WANG M C. Bearing capacity of strip footing above void[J]. Journal of Geotechnical Engineering，1983，109（1）：1-14.

[13]BADIE A，WANG M C. Stability of spread footing above void in clay[J]. Journal of Geotechnical Engineering，1984，110（11）：1591-1605.

[14]JAO M，WANG M C. Stability of strip footings above concrete- lined soft ground tunnels[J]. Tunnelling and Underground Space Technology，1998，13（4）：427-434.

[15]孫映霞，张智浩，张慧乐.岩溶区桩基稳定性影响参数敏感性分析[J].地下空间与工程学报，2013，9（2）：297-303.

SUN Y X，ZHANG Z H，ZHANG H L.Analysis on parameter sensitivity of pile foundations stability in Karst areas based on the theory of grey relation[J]. Chinese Journal of Underground Space and Engineering，2013，9（2）：297-303.（In Chinese）

[16]KIYOSUMI M，KUSAKABE O，OHUCHI M，et al.Yielding pressure of spread footing above multiple voids[J]. Journal of Geotechnical and Geoenvironmental Engineering，2007，133（12）：1522- 1531.

[17]LAVASAN A A，TALSAZ A，GHAZAVI M，et al. Behavior of shallow strip footing on twin voids[J]. Geotechnical and Geological Engineering，2016，34（6）：1791-1805.

[18]ZHOU H Z，ZHENG G，HE X P，et al.Bearing capacity of strip footings on c-φsoils with square voids [J]. Acta Geotechnica，2018，13（3）：747-755.

[19]LEE J K，JEONG S，KO J. Undrained stability of surface strip footings above voids[J]. Computers and Geotechnics，2014，62：128-135.

[20]WU G Q，ZHAO M H，ZHAO H，et al.Effect of eccentric load on the undrained bearing capacity of strip footings above voids[J]. International Journal of Geomechanics，2020，20（7）：04020078.

[21]WU G Q，ZHAO M H，ZHANG R，et al.Ultimate bearing capacity of eccentrically loaded strip footings above voids in rock masses [J]. Computers and Geotechnics，2020，128：103819.

[22]PHOON K K，KULHAWY F H. Characterization of geotechnical variability[J]. Canadian Geotechnical Journal，1999，36（4）：612-624.

[23]ALI A，LYAMIN A V，HUANG J S，et al.Undrained stability of a single circular tunnel in spatially variable soil subjected to surcharge loading[J]. Computers and Geotechnics，2017，84：16-27.