抛填块石土是一种物质组成复杂、结构分散、块石分布不均、透水性强、粒间结合力弱的岩土体，其力学特性区别于一般的土体或岩体^[1]。抛填块石土边坡因其不均匀性导致其物理力学性质复杂，非线性强，受到块石–土体空间分布影响大，利用传统的均质土体边坡研究手段难以准确地对其进行研究。近年来，国内外学者对块石土各方面特性进行了诸多研究，取得了很多具有实际工程价值的成果。徐文杰等^[2]基于前人研究，提出具体的块石土概念、划分细则和工程意义。Medley和Yue等^[3-5]利用剖面成像技术研究了块石分布对块石土边坡稳定性和混凝土剪切强度的影响。周中等^[6]采用现场试验、室内试验、理论分析及数值模拟相结合的研究方法，对块石土边坡的流–固耦合特性及其预测进行了探讨。油新华等^[7]利用野外大面积水平推剪试验，得出块石土的变形特点和相关的抗剪强度参数。何春灿等^[8]基于数字图像技术，建立软硬岩石模板库，并利用PFC2D颗粒流软件，对不同硬度岩石存在时的土石混合体细观损伤模式与力学特性进行了分析。徐文杰等^[9]基于数字图像技术，提出土石混合体边坡有限元模型的建立方法，并对其稳定性进行了研究。张森等^[10]运用圆形颗粒填充的方法研究了不同含石量对边坡稳定性的影响。张抒等^[11]基于图像处理技术和最大似然估计法，研究了块石长短尺寸和长短轴比的分布规律。

现阶段块石土边坡稳定性研究主要有模型计算法和模型试验法。模型试验法因成本高、随机性差等因素，无法进行大规模重复性试验，故要获取块石土边坡特性，有限元模型计算法是一种更为合理且有效的计算方法。徐文杰^[9]和张森^[10]等均采用有限元分析法，得到了较好的分析结果，但由于分析样本较少，且未分析形状参数等对边坡稳定性及破坏形式的影响。本文基于作者团队提出的随机生成方法，对多种块石几何（形状、粒径和空间）参数对块石土边坡的稳定性和破坏形式的影响进行有限元模拟研究。

作者团队提出的块石形状数据库随机生成方法^[12]基于数字图像处理技术采集块石形状，并利用数据库方法进行随机生成，是一种更高效、更准确的块石土细观模型随机生成方法。

为满足随机生成时形状参数的多样性，首先，采集大量块石形状数字图像并进行预处理；然后，将其批量导入MATLAB进行边界坐标化、光滑化、规则化等处理，并获取块石长短轴比等形状参数；最后，输出块石形状参数和处理后的块石坐标，建立数据库。基于数据库，根据需要设定参数阈值，进行随机生成。图1为基于块石数据库块石土模型生成技术流程^[12]。

图1(Fig. 1)

图1 随机生成实现流程图[12] Fig. 1 Flowchart of random generation[12]

该方法满足了土石混合体随机建模所要求的空间位置随机性与块石形状多样性，可直接生成符合现实块石形状并具有良好随机性的土石混合体模型；生成效率较高，所需时间根据含石量不同，一般耗时为5～120 min；同时，输出结果为块石节点坐标数据.txt文件集，可对应多种分析程序建立导入接口，兼容性极强。

在块石数据库中，包含3种块石形状表征参数，分别为长短轴比S_ab、面积差异比S_cp、棱角指数 $AI $ 。各参数计算方法如式（1）～（3）所示，图中描述见图2，部分块石形状及其参数见图3。

图2(Fig. 2)

图2 不同岩石形状参数计算图 Fig. 2 Calculating sketch of various rock’s shape parameters

图3(Fig. 3)

图3 块石模型及其形状参数 Fig. 3 Block model and its shape parameters

长短轴比S_ab：

式中，X_a为块石的长轴，X_b为块石的短轴。长短轴比S_ab可用于描述块石的条状性。

面积差异比S_cp：

式中： $A_{\rm{c}} $ 为块石边界与面积当量椭圆边界围成的面积； $A_{\rm{p}} $ 为面积当量椭圆面积，其值等于块石面积。S_cp可用于描述块石形状与同面积椭圆的差异。

棱角指数 $AI $ ：

式中，l_i为第i–1个角点和第i+1个角点连线的长，d_i为角点i到相邻角点之间连线的距离。棱角系数越大，表示该角点越尖锐。

除以上参数外，块石土模型颗粒参数（含石量、最大粒径）和空间参数（块石倾角、空间分布）也是块石土边坡稳定性及破坏形式的重要影响因素。根据所采用的随机生成算法特点，生成多个不同含石量、不同块石最大粒径及不同块石倾角的块石土边坡模型。

边坡模型采用文献[13]中的算例，边坡尺寸如图4所示。

图4(Fig. 4)

图4 边坡模型尺寸 Fig. 4 Size of slope

根据文献[14]，本文采用的工程特征尺寸L_c为坡高10 m，土–石阈值d_thr=0.05L_c=0.5 m（d_thr为土–石粒径界限值，当颗粒粒径小于该值视为土体，当颗粒粒径大于该值视为块石）。块石最大粒径视情况而定：当进行长短轴比分析时，为减小空间分布的影响，最大粒径d_max为0.1L_c，即1.0 m；当进行块石最大粒径分析时，d_max随分析情况改变，其余情况d_max取0.2L_c，即2.0 m。块石颗粒尺寸分布如图5所示。

图5(Fig. 5)

图5 块石颗粒尺寸分布 Fig. 5 Particle size distribution

选用ABAQUS进行有限元建模与分析，该软件较强的非线性计算能力，对于块石土这种强非线性结构，可靠性较高。由于输出文件集为块石节点坐标.txt文件，需利用python语言建立子程序将坐标导入CAE进行建模。

岩土体本构模型均采用摩尔库伦本构模型，具体参数如表1所示。模型网格采用ABAQUS内置网格CPE4，即四边形平面应变网格。网格播种尺寸为0.15 m，局部应力集中区域增大网格密度，减小网格尺寸。总单元数根据模型的不同而变化，一般为20 000～30 000个单元。文献[15]认为，ABAQUS中摩擦接触模型可良好地应用于块石土中土–石界面模拟；根据文献[16]对土体外摩擦系数的研究，本文认为块石土土–石接触面模型可定义为黏结模型^[17]，界面黏聚力主要由土体自身黏聚力发挥，而摩擦系数设置为0.363。各种不同参数模型及其网格划分如图6所示。

表1(Tab. 1)

表1 岩土体参数 Tab. 1 Parameters of rock and soil

表1 岩土体参数 Tab. 1 Parameters of rock and soil 成分 密度/ (g·cm–3) 弹性模量/ MPa 泊松比 $\;\mu $ 内摩擦角/ (°) 黏聚力/ kPa 块石 2.42 1 090 0.2 36 1 200 土体 1.81 50 0.3 20 20

图6(Fig. 6)

图6 块石土边坡有限元模型 Fig. 6 Soil-rock mixture slope FEM models

采用强度折减法计算边坡稳定性，通过添加场变量^[18]，ABAQUS能很好地实现有限元强度折减法，且计算结果较为准确。

计算考虑重力荷载，并做地应力平衡，即在重力作用下完成沉降过程。模型两侧固定水平位移，底部固定水平、竖向位移。可靠性分析选取50%含石量边坡进行。为了对比分析，同形状均质土边坡也进行了计算。

根据文献[19]中对边坡有限元稳定性计算中安全系数取值的讨论，本文于坡面底部、中部和顶部各一点提取水平位移，计算平均值作为判定位移，并选取折减系数判定位移曲线拐点为极限点，对应折减系数为安全系数。选取极限点下一计算步对应点为破坏点（failure point），进行横向位移及塑性区云图提取。图7为均质土边坡和50%含石量块石土边坡计算结果及对应模型破坏点处最大主应力场和横向位移场。

图7(Fig. 7)

图7 均质土边坡和50%含石量块石土边坡折减法计算结果 Fig. 7 Strength reduction calculation results of homogeneous soil slope and 50% rock content soil-rock mixture slope

均质土边坡极限点安全系数为1.36，与理论计算结果相符。如图7所示：均质土边坡最大主应力场呈现明显的分层分布，由坡面向坡脚呈现递减趋势，且于坡脚处出现应力集中情况；50%含石量块石土边坡最大主应力场总体上仍呈分层递减趋势，但内部分布极不均匀，且在土–石界面出现应力集中，与文献[10]中结果相似。根据文献[14]，横向位移场能很好地反映边坡滑动面位置。均质土边坡滑动面较为光滑，与理论计算滑面基本相符；块石土边坡滑面呈现明显的折线型，滑面形状较均质土边坡有明显变化。根据文献[20]水平推剪试验结果可知，块石土边坡滑面模拟情况与实际情况相符合，说明采用数据库随机生成模型，并利用ABAQUS进行折减计算的方法可靠度较高。

含石量是影响块石土边坡稳定性及破坏形式的重要因素。本文采用多模型计算结果取平均值的方法研究块石土边坡含石量对稳定性的影响。随机生成20%、30%、40%、50%、60%、70%共6级含石量块石土边坡，每级含石量随机生成5个边坡模型，分别编号为1^#、2^#、3^#、4^#、5^#，共30个块石土边坡模型进行分析。

图8为块石土边坡安全系数平均值与块石含石量关系曲线。随着边坡含石量的增加，边坡安全系数平均值呈上升趋势，且上升幅度在20%～60%含石量区间随着含石量增大呈上升趋势。其中：当含石量从50%增大到60%时，安全系数上升幅度最大；当含石量从60%增大到70%时，安全系数上升幅度大幅降低，此时边坡安全系数趋于稳定，随含石量上升变化不大。

图8(Fig. 8)

图8 块石土边坡安全系数平均值与块石含石量关系 Fig. 8 Relationship between average values of soil-rock mixture slope safety factor and the content of rock

在含石量相同的情况下，5个模型的安全系数也有一定的差异，如图9所示。本文采用同含石量模型安全系数标准差描述同含石量块石土边坡稳定性差异（图10）。传统的边坡安全系数计算方法受滑面形状影响较大，滑面形状发生小的改变，安全系数可能有较大差别，所以块石分布对于滑动面形状、边坡稳定性的影响较大。由图9、10可知：当含石量较小时，如20%～30%，块石难以集中分布，对滑动面形状影响较小，安全系数与均质土边坡相近，且安全系数标准差较小，同含石量边坡稳定性差异较小；当含石量为40%～50%时，块石分布影响逐渐增大，安全系数标准差也逐渐增大；当含石量增加到60%时，块石已逐渐代替均质土体作为边坡骨架存在，此时块石对于滑面形状影响达到最大值，安全系数标准差也达到峰值。

图9(Fig. 9)

图9 不同含石量块石土边坡安全系数点状图 Fig. 9 Scattergram of soil-rock mixture slope safety factor variation according to the rock block content

图10(Fig. 10)

图10 不同含石量块石土边坡安全系数标准差 Fig. 10 Standard deviation of soil-rock mixture slope safety factor variation different rock block content

为分析块石最大粒径对块石土边坡稳定性的影响，建立4种不同最大粒径的50%含石量抛填块石土边坡模型各5组，最大粒径d_max分别为1.0、1.5、2.5、3.0 m，颗粒级配控制方法与第3.1节模型相同，命名规则为“1.0 m（最大粒径） 1^#（模型编号）”，以此类推。将上述模型计算结果结合粒径为2 m的50%含石量边坡模型分析结果进行对比。图11为不同d_max时块石土边坡安全系数平均值分布，图12为不同d_max时块石土边坡安全系数点状图，图13为不同d_max时边坡安全系数标准差分布图。

图11(Fig. 11)

图11 不同dmax时块石土边坡安全系数平均值分布 Fig. 11 Mean value distribution of safety factor of soil-rock mixture slope with different dmax

图12(Fig. 12)

图12 不同dmax时块石土边坡安全系数点状图 Fig. 12 Scattergram of soil-rock mixture slope safety factor variation according to the dmax

图13(Fig. 13)

图13 不同dmax时块石土边坡安全系数标准差 Fig. 13 Standard deviation of safety coefficient of soil-rock mixture slope with different dmax

由图11～13可知：边坡稳定性随块石d_max的增大呈先减小后增大的趋势。当块石d_max较小，对比较大d_max边坡，其块石数量更多，块石分布较为均匀，块石整体分布对于边坡稳定性影响较小，故d_max为1.0与1.5 m的边坡安全系数标准差较小。随着d_max增大，块石分布影响逐渐增大，安全系数标准差开始增大，但此时单个粒径块石不足以决定滑动面形状，故安全系数最小值减小，平均值也下降，且于d_max=2.0 m时达到最小。当d_max继续增大，此时块石分布影响继续增加，单个“巨石”已足以影响或改变滑面形状，故相同d_max组安全系数标准差继续增大，安全系数极小值却开始增大，且安全系数平均值也达到最大。

块石空间分布是影响块石土边坡稳定性和破坏形式的重要因素。根据第3节分析，随着块石粒径或含石量增大，边坡安全系数标准差呈增大的趋势，主要是由于块石空间分布对安全系数的影响随着块石最大粒径或含石量的增大而增大。但影响的具体情况，仍需详细分析，故基于第3节计算模型，针对块石空间分布对块石土边坡影响的具体形式进行分析。

图14为部分边坡极限点下一计算步的横向位移分布。

图14(Fig. 14)

图14 部分块石土边坡破坏时的水平位移分布 Fig. 14 Failure horizontal displacement distribution of partial soil-rock mixture slope

根据含石量和块石最大粒径的变化，分3种情况进行讨论：

1）大粒径块石影响。如图14（a）、（b）所示，此类边坡滑动面上皆有大粒径块石阻挡滑动，存在“绕石效应”^[14]，且没有块石剪穿情况，滑面向边坡内部变化，滑动面形状折线化，不利于滑动发生，边坡整体稳定性大幅提高。

2）小粒径块石影响。如图14（c）、（d）所示，此类边坡滑动面上一般分布有较密集的小粒径块石，此时小粒径块石形成整体，滑动面形状折线化，不利于边坡滑动；由于小粒径块石形成的整体抗剪强度小于同体积单个大粒径块石，故边坡稳定性较情况1）略低。

3）低（无）影响。如图14中（e）、（f）所示，此类边坡滑动面一般没有能有效阻止边坡滑动的块石或块石群，或者由于块石空间位置分布导致滑面向边坡临空面发展，滑面一般较为光滑，与均质土体相近，稳定性也与均质土边坡相近。此类情况一般有两种可能性：①含石量低；②有一定含石量但滑动面位置没有良好的阻挡，或滑面形状被改变。

块石倾角对于块石土的物理力学性质影响较大^[21-22]。本文针对块石倾角不同，建立4种不同块石倾角的块石土边坡模型，含石量40%，长短轴S_ab取1.5～2.0 m，块石长轴与坡面夹角可定义为块石相对倾角θ，分别为0°、45°、90°、135°，旋转方向为逆时针。

图15为不同块石倾角的块石土边坡破坏时的水平位移分布，图16为不同块石倾角块石土边坡破坏时的塑性区分布。

图15(Fig. 15)

图15 不同块石倾角的块石土边坡破坏时的水平位移分布 Fig. 15 Failure horizontal displacement distribution of soil-rock mixture slopes with different inclinations of rock

图16(Fig. 16)

图16 不同块石倾角块石土边坡破坏时的塑性区分布 Fig. 16 Failure plastic zone distribution of soil-rock mixture slopes with different inclinations of rock

由图15、16可知：当θ=90°、135°时，滑动面与块石长轴方向相交，滑动面的发展受到块石的阻挡，此时塑性区呈发散状，滑面呈锯齿状，安全系数较大；当θ=0°、45°时，滑动面与长轴方向部分重叠，土–石界面强度低，为滑动面发展提供条件，塑性区与块石滑动面重叠，滑动面光滑，安全系数较小。

块石长短轴比对块石土结构的物理力学性质有很大影响。本文计算3种由不同长短轴比块石组成的块石土边坡，长短轴比S_ab取值范围分别为1.0～1.5、1.5～2.0、2.0～2.5 m。含石量设置为40%，块石的最大粒径取1 m，棱角指数AI为1.0～1.5，面积差异比S_cp的取值范围为0.10～0.15。在重力作用下沉积的颗粒，其长轴方向往往水平，具有方向排列性形成层面，故块石长轴设置为水平，边坡模型与网格的划分如图6所示。

图17为块石长短轴比不同时块石土边坡破坏点横向位移分布，图18为不同长短轴比块石土边坡破坏时的塑性区分布。

图17(Fig. 17)

图17 不同长短轴比块石土边坡破坏时的水平位移分布 Fig. 17 Failure horizontal displacement distribution of soil-rock mixture slope with different axial ratio of rock

图18(Fig. 18)

图18 不同长短轴比块石土边坡破坏时的塑性区分布 Fig. 18 Failure plastic zone distribution of soil-rock mixture slopes with different axial ratio of rock

由图17、18可知：随着块石长短轴比的增大，块石土边坡稳定性降低，滑动面越发平滑，塑性区范围也更加集中。这是因为块石土边坡最薄弱面为土–石接触面，而边坡在破坏时会于坡脚出现应力集中，即坡脚处的强度对于整体稳定性有较大影响。当长短轴比增大时，块石处于水平状态，坡脚处接触面面积增大，且接触面更为平滑，强度降低，故稳定性下降。

经过模型计算分析，棱角系数AI和面积差异比 S_cp对块石土边坡稳定性及破坏形式影响只有在粒径足够大时才能有所体现，限于篇幅，不再进行详细论述分析。

基于块石形状数据库，考虑块石多种几何参数，实现了块石土边坡模型随机生成，并利用ABAQUS有限元软件进行边坡稳定性折减计算，分析了块石粒径参数、空间参数、长短轴比对块石土边坡的特性影响，得出以下结论：

1）含石量对块石土边坡稳定性的影响。含石量增大，边坡安全系数上升，且含石量在50%～60%区间内上升速率最快。含石量增大，同粒径组边坡安全系数标准差增大，模型安全系数差异变大。

2）最大粒径d_max对块石土边坡稳定性的影响。最大粒径增大，块石分布对于边坡安全系数的影响增大，边坡稳定性先减小后增强，同粒径组边坡安全系数标准差增大。

3）块石空间分布对块石土边坡稳定性的影响。含石量和最大粒径增大，块石空间分布对于边坡稳定性影响增大。粒径大小不同，块石空间分布对于块石土边坡稳定性影响方式也不同，大粒径块石存在令滑动面向边坡内部发展，而密集小粒径块石大多是令边坡呈折线形，二者皆可增大边坡稳定性。

4）块石倾角对块石土边坡稳定性的影响。当块石长轴与坡面夹角为90°、135°时，边坡塑性区呈发散分布，滑面更靠内部，稳定性更好；当块石长轴与坡面夹角为0°、45°，与上述情况相反。

5）块石长短轴比对于块石土边坡的影响。在通常沉积条件下，边坡稳定性随长短轴比增加而降低，滑动面越发平滑，而塑性区范围也更加集中。