基于谱表示法和时间衰变模型的岩体稳定性评价

吴川 廖坤阳 孙大齐

摘要：为了从时间尺度和空间变异性上研究石质文物危岩体的结构稳定性，建立了岩体强度参数的空间分布和时间退化模型。选取麦积山石窟13号窟危岩体为研究对象，利用强度折减法对石窟岩体稳定性进行模拟计算。模拟结果表明：① 麦积山13号窟现状条件下最大位移为1.8 mm，中部位移方向略微倾向石窟外侧，石窟顶部呈现局部塑性变形，但整体稳定性较好，模拟结果与现状调查较为一致。② 加入时间衰变模型后发现，13号窟现状条件下顶部位移为1.2 mm;预测50 a后洞窟最大位移位于洞顶，为1.3 mm;100 a后洞窟最大位移从洞顶转移至坡脚，为16.7 mm。③ 强度折减分析表明，100 a后模型产生大面积塑性区域贯通现象，洞顶发生拉张破坏，洞底剪切破坏。时空变异模型有效提高了传统静力学计算精度，通过时间退化模型真实地再现了石窟危岩体的动态破坏过程，在精确评价危岩体长期稳定性方面具有广泛应用价值，可为危岩治理提供参考。

关 键 词：

麦积山石窟; 时空变异性; 非平稳随机过程谱表示法; 危岩体稳定性; 石质文物

中图法分类号： TU47

文献标志码： A

DOI：10.16232/j.cnki.1001-4179.2021.06.013

危岩体的稳定性评价和分析一直是岩土工程、地质工程领域的重要研究方向[1]。目前危岩体稳定性分析方法主要分为两类：第一类以刚体极限平衡法为基础，如陈洪凯等[2]将三峡库区危岩体失稳破坏分为滑塌式、倾倒式和坠落式，并在考虑自重、地震和裂隙水压力的情况下给出了3种危岩体的静力学计算方法。第二类评价方法以数值方法（如限元法、有限差分和离散元）为基础，如郑颖人等[3]将有限元强度折减法应用于岩土边坡稳定性评价中，开创了求解岩质边坡滑动面与稳定安全系数的先河。前一类方法过于简便。第二类评价法凭借计算机强大的运算效率划分单元体进行稳定性评价，却对于岩体性质的变异性考虑不足，如岩体参数（凝聚力、内摩擦角、重度等）在空间分布的不均匀性和时间上是不断变化的。Symbol`@@石质文物由于遭受长年累月的风化侵蚀作用，整体稳定性问题是文物保护工作的重点，同时由于其文物属性无法采取大量破坏性的工程防护措施，因此其在时间和空间尺度上的稳定性研究具有现实意义。

目前，考虑时空变异性的石质文物危岩稳定性评价的成果较少，大部分研究成果仍处于现象描述阶段[4-7]。本文以麦积山石窟13号窟为研究对象，综合以上研究成果，利用谱表示法所生成的随机场来模拟岩体抗剪强度参数空间变化特征。作为一种尝试，建立基于强度衰减规律的时间变化模型。借助有限元强度折减法动态分析了石窟危岩体稳定性，并针对性地提出了适用于石质文物保护的防治方案。

1 空间参数分布模型

各种成因的环境地质作用会导致岩体在空间分布上出现差异性[8]，如发育于岩体中结构面的几何分布特征，对于岩体强度、变形和渗透性的影响不容忽视。由于岩体强度参数既具有变异性又具有相关性，无法用简单线性函数表示，所以选择非平稳随机过程谱表示法进行模拟。谱表示法于1990年被Shinozuka等用于模拟生成随机场。本文采用了梁建文[9]推导的非平稳随机过程蒙特卡洛模拟的一个谱表示法，该方法产生的样本函数能够较准确反映岩体参数在水平方向的延续性与垂直方向上的变异性。

岩体强度参数是一个切实存在的变量（不能小于零）。依据摩尔-库伦剪切强度准则可知，对于岩体稳定性影响最为显著的参数是凝聚力和内摩擦角，而其他参数如容重、弹性模量以及摩擦角等空间变异性不明显，对于岩体稳定性分析贡献很小。考虑到数值计算效率，本文假定其他参数恒为常量，选用二维正态分布的非平稳随机场来表示基本符合模拟要求。二维对数正态分布表达式如下：

fx，y=exp{ln1+V2×

M-1i=0N-1j=0

σij[Vijθcosω1ix+ω2jy+

Wijθsinω1ix+ω2jy]+

lnμ-ln1+V22}（1）

式中：x和y为单元体空间坐标;V，μ为别为凝聚力和内摩擦角的标准差和平均值;σij为（iM+j+1）项的标准差;Vij（θ），Wij（θ）分别为相对独立的标准正态分布随机变量;ω1i，ω2j则分别为频率坐标。

其中，σij的计算公式如下：

σij=4Sω1i，ω2jΔω1Δω2（2）

式中：S為单边功率谱密度函数;Δω1，Δω2分别为频率坐标轴的离散区间。

由于相关函数对于岩体结构稳定性指标影响不明显[9]，对石质文物稳定性模拟中选择计算量较小的二维指数函数作为相关函数：

ρx，y=exp-2xδx+2yδy（3）

式中：ρ代表自相关系数;x，y代表水平和垂直方向距离，m;δx和δy分别代表水平方向和垂直方向的波动范围，m。

现场测定样本越多，就越能模拟出与实际情况相近的空间分布情况，理论上此种模拟方法可以获得无限接近于真实情况的解，但在实际中是不可能的。通过谱表示法，根据有限个关键点出测得的强度参数（c，φ）统计值进行随机场模拟，从而反演出可能的岩土体强度参数空间分布情况是可行的。由于文章篇幅限制，具体过程可参见文献[8]。

2 时间衰变模型

描述岩石力学特性随时间变化的规律是危岩体稳定性分析中的难题，从时间尺度上评价岩体的稳定性研究尤为重要。许宏发等[10]通过实验定量研究了软岩的长期强度特征，认为岩体存在一组初始强度和最终强度，随时间推移初始强度逐渐降低并无限接近于最终强度，导致这一变化过程是岩体内壁微观损伤积累的结果。本文根据李连崇[11]关于岩体细微损伤随时间变化试验研究，尝试建立了以指数衰减的力学退化模型（见图1），用以描述危岩体的长期强度。

从图1的退化模型中可知，岩体的强度（σt）随时间逐渐衰退，初始强度（σ0）和长期强度（σl）固定不变;当t=0时，σt=σ0，当t→∞时，σt=σl。

由此将岩体微观表征单元强度的时间退化模型用式（4）表示：

σt=σl+σ0-σle-kt（4）

式中：系数k和长期强度σl与岩体赋存环境有关。k反映强度衰减的快慢程度，若岩体所处环境有风化剥蚀、地震等劣化岩体结构的因素存在，那么系数k势必会较大，同时σl的大小受环境条件与岩性条件共同控制。对于以上两种量值的选取需根据现场和室内相关标定试验来进行。

3 基于時空变异模型的评价方法

在将时间和空间变异模型加载之前，用强度折减法自动搜索滑动面，且应用数值模拟软件对其稳定性进行分析[12-13]。本文将新鲜岩体和强风化岩体的强度指标作为初始强度σ0和长期强度σl。对于空间分布参数随机场的模拟，选取麦积山石窟关键点的强度指标输入谱表示计算法中，根据关键点强度参数分布规律，谱表示法反演出拟合最好的岩体整体空间参数分布情况。

（1） 基础模型建立。首先在ANSYS中对模型进行剖分。通过程序转换将建立好的三维模型导入FLAC3D。再在FLAC3D中对模型区域、边界条件、材料属性和模型等进行详细设置。最后将完整的模型文件保存为后缀为xml的接口文件。

（2） 空间参数设定。根据所推演出的谱表示法模拟出不同方向（水平和垂直）波动范围的随机场，得到c和φ的n组随机组合。然后将模型文件中的强度参数用以上n组数据替代，从而得到基于空间变异的新模型。

（3） 时间参数设定。根据强度随时间衰变模型，将上一步所推导出的n组强度参数根据衰变公式折减后得n组新数据，基于这n组新数据的模型即为t时刻的危岩体三维模型，其具有时间和空间上的变异性。

（4） 稳定性计算。利用WinbatchTM（批处理软件）赋值，在FLAC3D中对新模型文件的稳定性进行动态分析，得出计算结果。

4 算例分析

4.1 模型尺寸与网格划分

采用三维激光扫描建立麦积山石窟的三维地形曲面。为了确保建模过程中模型不受网格离散化和建模范围的影响，本文选取的模型范围超过研究洞窟尺寸的2.5～3.0倍。依据三维激光扫描确定的模型曲面较为复杂和贴合实际情况，本文通过ANASYS采用四面体单元进行建模处理（见图2），模型建好后导入到FLAC3D中进行数值分析。

由于麦积山石窟13号窟为南北向布置，为了便于计算，本文的模型坐标系设定竖直方向为Z轴，水平方向为X和Y轴，其中X轴正方向朝南，Y轴正方向朝东，石窟边坡岩体走向与Y方向平行。最终确定模型尺寸为垂直方向（Z）高50 m，南北方向（X）长75 m，东西向（Y）长35 m。另外，考虑到石窟岩体表面性状复杂，可能对后续的模型计算效率产生影响，因此本文对石窟中造像进行了一定程度简化，并不影响石窟整体稳定性计算结果，简化后的模型见图2。

4.2 边界条件与模型参数

模型计算过程中设定顶面和底面为自由边界，模型南侧（前侧）同样为自由边界，侧面和后面为固定约束。破坏准则参考摩尔-库仑破坏准则。现场选取石窟基岩砂岩进行室内物理力学试验。由于研究对象为石质文物，不宜进行大规模破坏性勘探，本文利用轻型便携式钻机沿石窟顶面进行钻探，共布设钻探剖面6条，钻孔数量24个，钻孔孔深2.0～5.0 m，通过对岩芯进行力学试验发现弱风化岩体饱和单轴抗压强度一般为30.0～40.0 MPa，强风化岩体一般为12.6～18.7 MPa，相较而言强风化岩体强度较弱风化强度衰减约38.1%～47.5%。现场调查发现石窟岩体以弱风化和强风化岩体为主。模型建模所用岩体力学试验结果见表1。

石窟模型在空间上的力学参数通过谱表示法来表征，根据有限个关键点测出的抗剪强度参数统计值进行随机场模拟。由于岩体内变形模量、泊松比和密度空间变化不大，因此将其简化为常量。本文通过现场钻探取样过程中确定的关键点的物理力学特性指标，将其换算成有关二维正态分布的函数值，其计算参数见表2。

4.3 模拟成果分析

4.3.1 现状位移分布规律

将空间分布模型导入麦积山石窟模型力学参数中，用FLAC3D进行麦积山石窟的位移和应力分布规律分析，其模型的总位移矢量和等值线云图见图3。

由图3可知，在自重应力作用下，石窟总体的位移呈现从顶部到底部均匀降低的规律：在山顶部最大位移为1.8 mm，中部位移方向略微倾向石窟外侧，底部位移基本为0，说明该模型自重应力分布规律基本合理。

4.3.2 时间尺度预测

在石窟整体位移分析的基础上，为了更全面分析和研究洞窟的时间衰变规律，对洞窟中间位置切割一处典型剖面，依据时间衰变模型设定强度衰减系数k为0.02，预测了13号窟随时间的应力及位移变化规律。

图4（a）为t=0现状条件下的石窟中轴线剖面位移云图，可以很明显发现剖面的总位移方向向模型下后方倾斜，并且最大位移位于石窟顶部，其值为1.8 mm，平台底部位移为0。13号洞窟顶部位移为1.2 mm，底部位移为0.6 mm，位移方向均为下后方。

图4（b）代表t=50 a即50 a后的石窟剖面位移云图。对比发现模型总位移呈增大趋势，最大值0.3 mm，同时平台底部位移依旧为0，说明该段时间石窟稳定性依旧较好。

当t=100 a时，模型顶部位移几乎无变化，最大位移从石窟顶部转移到底部平台，石窟山脚的最大位移达到16.7 mm，13号洞窟顶板的沉降同样达到13.7 mm，说明石窟已经发生了不同程度的局部失稳破坏。现场调查情况发现，13号洞窟砂岩岩体受风化和裂隙渗水侵蚀较为严重，已经存在不同大小的风化凹槽，若不采取措施，可以预见100 a后石窟必然会发生一定程度的稳定性破坏问题。

模拟发现50 a内石窟山顶和山脚处应力集中尚不存在，但是随着时间推移应力逐渐在两个位置集中。在100 a后模型中山体向下应力逐渐增大，洞窟顶板和前壁洞口均发生了应力集中现象。

4.3.3 稳定性时空评价

加入空间参数分布模型和时间衰变模型后，强度拆减理论计算量大大超过传统计算方式。因此本文将折减系数限定在1.0～3.0进行试运算。研究发现，当t=0，R=1.42时，模型顶部的局部区域开始发生塑性破坏，说明现状条件下13号窟整体稳定性尚可，但是石窟顶部可能发生部分破坏。当t=100 a，R=1.10时，模型开始发生大面积塑性区域贯通现象。如图5所示，在洞口上部黄色区域为大面积的拉张破坏，在洞窟底部红色区域为剪切破坏，因此可以预测100 a石窟整体处于临界破坏，随时可能发生大规模的顶部岩体崩塌现象。

对于石质文物岩体的稳定安全性系数，目前并没有统一的认识。参照现阶段工程岩体边坡相关规范，考虑到文物的价值属性，石窟属于重要保护文物，防治工程等级按一级考虑，危岩体防治安全系数应为1.40～1.60。目前來看13号窟整体处于稳定状态。

5 结 论

（1） 岩体由于各种成因地质作用，其强度参数（c，φ）在空间分布上具有不确定性，而空间分布参数模型的引入克服了传统静力学计算中这一缺点，使得计算更加精确和符合实际情况。石质文物岩体由于岩体流变特性，其强度参数随时间会逐渐降低，导致强度降低的主要因素是软弱结构面受风化剥蚀等内外地质作用。时间衰变模型基于实际勘测成果从整体上把握了岩体强度衰变规律，预测了麦积山石窟13号窟长时间的稳定性，取得良好效果。

（2） 麦积山石窟现状条件下最大位移为1.8 mm，石窟顶部呈现局部塑性变形，但整体稳定性较好，模拟结果与现状调查较为一致;加入时间衰变模型后发现洞窟顶部位移为1.2 mm，预测50 a后洞窟最大位移位于洞顶，为1.3 mm，100 a后洞窟最大位移从洞顶转移至坡脚，为16.7 mm;强度折减分析表明100 a后模型产生大面积塑性区域贯通现象，洞顶发生拉张破坏，洞底发生剪切破坏。

（3） 目前，将强度参数的时间和空间变异性结合起来评价岩体稳定性的研究尚处于初始阶段。本文只是考虑了凝聚力和摩擦角的不确定性，对于重度、变形模量、边界条件以及理论模型等假定其恒定不变，因此对于重要程度较高的复杂岩体，今后的稳定性评价应朝更精细、更复杂的方向发展。

参考文献：

[1] 唐辉明.工程地质学基础[M]北京：化学工业出版社，2008：15-21.

[2] 陈洪凯，唐红梅，王蓉.三峡库区危岩稳定性计算方法及应用[J].岩石力学与工程学报，2004，23（4）：614-619.

[3] 郑颖人，赵尚毅.有限元强度折减法在土坡与岩坡中的应用[J].岩石力学与工程学报，2004，23（19）：3381-3388.

[4] 姜效玺，吴琼，邓小龙，等.复杂形态立柱式石质文物的病害调查及防治对策：以浙江省柯岩云骨岩柱为例[J].工程勘察，2016，44（10）：1-5.

[5] 甄刚，秦立科，马宏林.乾陵无字碑裂隙分布情况调查及其稳定性分析[J].文博，2015（1）：94-98.

[6] 唐智亮.固结灌浆在保护云岗石窟石质文物遗址的应用研究[D].长春：吉林大学，2013.

[7] 田小甫.太原晋阳大佛边坡岩体风化分级及地震稳定性评价研究[D].北京：中国地质大学（北京），2009.

[8] 潘别桐.岩体结构面网络模拟及应用[M].武汉：中国地质大学出版社，1987：1-2.

[9] 梁建文，肖笛.非平稳随机过程模拟的一个谱表示方法[J].应用概率统计，2005，21（4）：375-386.

[10] 许宏发.软岩强度和弹模的时间效应研究[J].岩石力学与工程学报，1997，16（3）：246-246.

[11] 李连崇，徐涛，唐春安，等.单轴压缩下岩石蠕变失稳破坏过程数值模拟[J].岩土力学，2007，28（9）：1978-1982，1986.

[12] 赵尚毅，郑颖人，时卫民，王敬林.用有限元强度折减法求边坡稳定安全系数[J].岩土工程学报，2002，24（3）：343-346.

[13] 向娟，陈占锋，范文臣.考虑岩体节理密度影响的高陡边坡稳定性研究[J].中国科技论文，2018，13（24）：2790-2795.

（编辑：郑 毅）

Stability analysis on dangerous rockmass based on spectral representation

method and time degradation model

WU Chuan1，LIAO Kunyang2，SUN Daqi3

（1.Chongqing Industry Polytechnic College，Chongqing 401120，China; 2.Fujian Forestry Vocational & Technical College，Longyan 353000，China; 3.China Academy of Railway Sciences，Beijing 100081，China;）

Abstract：

In order to study the structural stability of stone cultural relic dangerous rock masses from time scale and spatial variability，a model of spatial distribution and time degradation of rock mass strength parameters was established.The dangerous rock mass of Grotto 13 in the Maiji Mountain Grottoes was selected as the research object.The strength reduction method was used to simulate and calculate the rock mass stability.The simulation results show that the maximum displacement under the current conditions of the Maiji Mountain Grotto 13 is 1.8 mm，displacement direction of the middle part is slightly inclined to the outside of the cave，and the top of the cave shows local plastic deformation，however the overall stability is good，and the simulation results are consistent with the current survey.By using the time degradation model，it is found that the top displacement of Grotto 13 is 1.2 mm under current conditions，and the maximum displacement of the cave is predicted to be 1.3 mm after 50 years，still locating at the cave top，however the maximum displacement of the cave will be transferred from the cave top to the slope foot after 100 years，reaching 16.7mm.It shows that the model will have a large plastic area penetration phenomenon after 100 years，top of the cave shows tensile failure，and the bottom of the cave shows shear failure.The spatio-temporal variation model effectively improves the accuracy of traditional static mechanicalcal culations，and the time degradation model truly reproduces the dynamic failure process of dangerous rock masses in caves.It has extensive application values in accurately evaluating the long-term structural stability of dangerous rock masses，and can be used as references for dangerous rock treatment project.

Key words：

Maiji Mountain Grottoes;spatio-temporal variability;spectral representation method of nonstationary random process;stability of dangerous rock masses;stone cultural relics