土石混合体边坡是指由土夹碎石或卵石以及碎石或卵石夹土等土石混合体组成的边坡，常常发育在第四系松散堆积层中，主要由风化残积、滑坡崩堆积以及冲洪堆积等形成^[1]。土石混合体边坡的二元结构是指坡体中除有土石混合体外下部还有基岩分布。土石混合体边坡的二元结构常常越接近地表其风化作用越强，随着深度的增加风化作用逐渐减弱，块石含量逐渐增多，直至新鲜基岩^[2]。

目前为止，有关土石混合体的研究进展较大，研究的手段主要是室内外试验及数值模拟。室内外试验方面，徐文杰等^[3–4]通过图像处理获取的块石粒度组成为基础制备试样，采用现场大型直剪试验研究了含石率对其强度和变形的影响；油新华、李晓等^[5–7]采用现场水平推剪及压剪试验研究了多个地方土石混合体的强度及变形特性；徐东升等^[8]利用大型单剪仪进行了21组大型单剪试验，分析了含石率对土石混合体抗剪强度、剪胀及剪缩特性的影响。数值模拟方面，丁秀丽等^[9]采用APDL语言生成了土石混合体随机结构模型，建立了非饱和土石混合体数值模拟方法；杨冰等^[10]采用PFC3D进行了土石混合体侧限压缩试验，研究了含石率对其骨架效应等结构性的影响；郑宏等^[11]采用数值流形法（NMM）建立了土石混合体的计算模型，研究了围压及土石错动对土石混合体抗压强度的影响。

土石混合体边坡稳定性分析方面，目前的块石建模方法主要是数字图像建模法和随机生成法两种。数字图像建模法方面，徐文杰等^[12]运用数字图像处理技术建立了西南地区一土石混合体边坡的细观结构模型，并运用强度折减法对其稳定性进行了分析研究。随机生成法方面，邵帅等^[13]采用有限元强度折减法探讨了块石空间分布对土石混合体边坡剪切带及稳定性的影响；李亮等^[14]采用MATLAB编制相应程序生成一定级配的块石数据，并采用FLAC3D内嵌的有限差分强度折减法对土石混合体边坡进行了稳定性分析。到目前为止，不管是数字图像建模法还是随机生成法，都未系统考虑块石含量及随机分布、基覆岩层倾角对土石混合体边坡稳定性的影响。

作者提出了一种可考虑不同块石含量、块石随机分布、基覆岩层倾角的土石混合体边坡稳定性分析方法。其主要的实施过程包括：根据土石混合体边坡中块石含量、级配、基覆岩层倾角等参数，随机投放不同直径圆形基料来生成土石混合体边坡初始模型，通过多边形最大边数的设定，获得带有棱角的多边形文件；再通过VBA的多线函数绘制出对应的块石边界轮廓线；最后将数据导入到Optum G2中建立土石混合体边坡模型，采用有限元极限分析法对不同块石含量及随机分布、基覆岩层倾角的土石混合体边坡进行了稳定性分析，并将计算结果与采用等效参数模型的计算结果进行了对比验证。研究结果可为土石混合体边坡的设计及施工提供参考依据。

1 随机块石模型生成系统

随机块石模型生成的主要的流程为：根据各块石级配下的块石面积及含量，生成随机元基料，控制块石的畸变系数及最大边数，最终确定满足级配特点、含石量及畸变系数随机块石模型。其中畸变系数（ $\mu $ ）定义为块石中最小、最大外接圆直径的比值，即为 ${{{d_{\min }}} / {{d_{\max }}}}$ 。

1.1 多边形块石的生成

根据土石混合体边坡中含石率及级配，先随机生成一个对应级配的随机圆，控制块石的位置和最大直径，其对应的字符串为 ${x_0}$ 、 ${y_0}$ 、 $d$ ；设置最大和最小直径角为：

${\theta _{\max }} = 2{\text{π}} rand$

(1)

${\theta _{\min }} = 2{\text{π}} rand$

(2)

式中， $ rand$ 表示0~1之间随机值。则对应的极值点坐标为：

$({x_0} \pm 0.5\delta \cos\;{\theta _{\max }},{y_0} \pm 0.5\delta \sin\; {\theta _{\max }})$

(3)

$({x_0} \pm 0.5\delta \cos\;{\theta _{\min }},{y_0} \pm 0.5\delta \sin\;{\theta _{\min }})$

(4)

式中， ${x_0}$ 、 ${y_0}$ 分别为随机圆的横坐标和纵坐标； ${\theta _{\max }}$ 、 ${\theta _{\min }}$ 分别为最大、最小直径角； $\delta $ 为畸变系数控制下的随机边长，计算公式为：

$\delta = (1 + k(2rand - 1))\mu d$

(5)

式中： $rand$ 为0～1之间随机值； $\mu $ 为畸变系数，取值范围为0～0.9； $d$ 为随机圆直径； $k$ 为畸变系数变化阀值，取值范围为0～0.1。

控制点生成后，完成块石辅助随机边界点生成，其直径随机函数和角度随机函数分别为：

${d_{rand}} = \delta + (d - \delta )rand$

(6)

${\theta _{rand}} = 2{\text{π}} rand$

(7)

式中， $\delta $ 为畸变系数控制下的随机边长， $d$ 为随机圆直径， $rand$ 为0～1之间随机值。

由此可以根据圆心点的坐标计算出不同的块石边界点坐标值，其计算次数应以设置的多边形最多边数 $(n - 4)$ 为控制标准；由于随机边界点角度随机，因此需对单个块石的所有边界点按照逆时针进行排序。

边界点有理化后的块石面积采用三角形求和法计算，对任意多边形A1、A2、···，以圆形 $({x_0},{y_0})$ 为三角形顶点，底边为 ${l_{ab}}$ ，则任一三角形面积计算公式为：

${s_i} = \frac{1}{2}\left| {\begin{array}{*{20}{c}} {{x_{0\; }}}&{{y_{0\; }}}&1\\ {{x_n}\; }&{{y_n}}&1\\ {{x_{n + 1}}}&{{y_{n + 1}}}&1 \end{array}} \right|$

(8)

则单个块石的面积值为 $ S=\displaystyle\sum S_{i}$ 。

随机块石分布程序的界面如图1所示。

1.2 多边形块石的投放 1.2.1 投放算法

关于多边形块石之间的入侵判断，文献[15]中提出了向量面积法，文献[16]主要采用夹角之和测试法，De等^[17]则采用了空间分割填充法。3种方法虽然都可以判定块石是否入侵，但计算均较为复杂，作者通过对投放流程的优化，采用直线距离法直接判定块石角点是否入侵随机圆内，简化了复杂的入侵判断流程；又通过投放级配由大到小的顺序，直接避免了随机圆在块石内部的可能性，极大的提高了块石投放效率。

1.2.2 投放流程

块石投放流程为：①输入初始参数，主要包括边坡轮廓函数，块石级配参数、含石量及畸变参数；②在规定级配内生成随机圆；③判断是否有块石入侵到随机圆中；④按照畸变系数生成多边形块石，并将单个块石节点数据有理化；⑤绘制对应的块石轮廓图，重复步骤；②③④⑤，直至系统生成满足含石率、级配、畸变系数的轮廓图。投放流程如图2所示。

1.2.3 不同含石率及基覆岩层倾角下土石混合体边坡块石的生成

土石混合体边坡坡高为10.0 m,坡度为45º，坡脚到左侧边界距离为10.0 m，坡顶到边坡右侧的边界距离为15.0 m，坡脚向下边界的距离为10.0 m。

参照文献[12, 18]，土石混合体边坡的工程特征尺度 ${L_c}$ 为10.0 m，相应的土–石阈值 ${d_{S/RT}} = $ 0.05 × 10 m = 0.5 m。取块石的最小粒径 ${d_{\min }}$ 为0.5 m、最大粒径 ${d_{\max }}$ 为2.0 m。土石混合体边坡中块石级配粒组如表1所示。

含石率（volumetric rock block proportion，VBP）定义为块石总面积与上部土石混合体边坡面积的比值。采用基覆岩层的倾角为10º，研究含石率对土石混合体边坡稳定性的影响，不同含石率下土石混合体边坡的示意图如图3所示。采用含石率为40%，研究基覆岩层倾角对土石混合体边坡稳定性的影响，不同基覆岩层倾角下土石混合体边坡的示意图如图4所示。

2 有限元极限分析计算模型的建立

将随机生成的块石模型导入到有限元极限分析软件Optum G2中形成土石混合体边坡数值计算模型。在模型的左右边界设置水平约束，在模型底部边界同时设置固定约束，有限元极限分析模型示意图如图5所示。

文献[12, 18–19]表明，土石混合体边坡中“块石”与“土体”应具有力学特性的差异，需满足 ${{{E_{\rm block}}} / {{E_{\rm matrix}}}} > 2$ 以及 ${{\tan {\varphi _{\rm block}}} / {\tan {\varphi _{\rm matrix}} > 2}}$ ，具体的计算参数如表2所示。土石混合体边坡稳定性计算时“块石”、“土体”、基岩的本构模型采用摩尔–库仑模型。

有限元极限分析法在计算土石混合体边坡稳定性时，采用强度折减法，边坡的安全系数为：

${F_s} = \frac{{\tan \;\varphi }}{{\tan \;{\varphi _{\rm cr}}}} = \frac{c}{{{c_{\rm cr}}}}$

(9)

式（9）中， ${c_{\rm cr}}$ 和 ${\varphi _{\rm cr}}$ 为土石混合体边坡处于临界滑动状态时岩土体的黏聚力和内摩擦角， $c$ 和 $\varphi $ 为岩土体的黏聚力和内摩擦角。

图6中给出了采用有限元下限法计算土石混合边坡安全系数的分析流程。安全系数上限法与此类似，上限解输出改为 ${F_{s\max }}$ 。 $TOL$ 为收敛精度，论文中取0.001。

3 不同含石率下边坡稳定性结果分析 3.1 不同含石率下边坡的安全系数

采用有限元极限分析法分别进行了土石混合体边坡中含石率为0%、10%、20%、30%、40%、50%及60%的稳定性分析计算，获得各含石率下土石混合体边坡的下、上限安全系数。为分析块石分布位置对边坡安全系数影响，每种含石率生成8个随机试样，共进行49个土石混合体边坡的稳定性分析。不同含石率下土石混合体边坡的下、上限安全系数分别如图7、8所示，土石混合体边坡的下、上限安全系数的平均值、标准差及变异系数如表3所示。

从图7可以看出，当含石率低于40%时，随着含石率的增加，土石混合体边坡的下限安全系数增加较小；当含石率大于40%时，边坡的下限安全系数显著增长。与均质土边坡下限安全系数相比，当土石混合体边坡的含石率分别为10%、20%、30%、40%、50%、60%时，其下限安全系数的平均值分别增加0.4%、1.1%、1.8%、3.0%、9.5%、16.4%。这与Maria等^[18]的有限元计算结果较为一致。

从图8中可以看出，当含石率低于20%时，随着含石率的增加，土石混合体边坡的上限安全系数增长较为缓慢；当含石率大于20%时，边坡的上限安全系数显著增长。与均质土边坡上限安全系数相比，当土石混合体边坡的含石率分别为10%、20%、30%、40%、50%、60%时，其上限安全系数的平均值分别增加2.8%、6.1%、11.7%、16.9%、36.4%、68.5%。产生这个结果的原因是，随着含石率的增加，土石混合体边坡的剪切带产生了绕石及分叉绕石效应，而不再是圆弧形。

从表3可以看出，相同的含石率下，由于块石空间位置的不同，土石混合体边坡的下、上限安全系数呈现出较大的差异性。当含石率从10%增加到60%时，土石混合体边坡下限安全系数的变异系数从0.6%增加到2.8%；土石混合体边坡上限安全系数的变异系数从2.0%增加到6.7%。可见，对于分析土石混合体边坡这种各项异性材料边坡的稳定性时，应采用随机的方法而不应该采用确定性方法。

3.2 不同含石率下边坡的剪切耗散图

选取不同含石率下土石混合体边坡的上限剪切耗散图，每种含石率选取一种块石位置，如图9所示。从图9可以看出，随着土石混合体边坡含石率的增加，坡体中的剪切带不再是圆弧型，而是呈现出不同的剪切带扩展模式。根据坡体中块石空间分布位置的不同，土石混合体边坡的剪切带呈现出3种典型的扩展模式，分别是“绕石”模式、“分流”模式以及“包含”模式，如图10所示。对于图10（a）中的“绕石”模式，将会引起坡体剪切带的偏移；图10（b）中的“分流”模式将会导致剪切带范围的扩大，剪切带变长；图10（c）中的“包含”模式是导致土石混合体边坡产生多滑动面的主要原因。本文的研究结果与文献[14, 18, 20]中剪切带扩展趋势一致。

4 不同基覆岩层倾角下边坡稳定性结果分析 4.1 不同基覆岩层倾角下边坡的安全系数

采用有限元极限分析法分别进行土石混合体边坡中基覆岩层倾角为0º、5º、10º、15º及20º的稳定性分析计算，获得各基覆岩层倾角下土石混合体边坡的下、上限安全系数。为分析块石分布位置对边坡安全系数影响，每种含石率生成8个随机试样，共进行40个土石混合体边坡的稳定性分析。不同基覆岩层倾角下土石混合体边坡的下、上限安全系数分别如图11、12所示。土石混合体边坡的下、上限安全系数的平均值、标准差及变异系数如表4所示。

从图11、12中可以看出，相同含石率下，随着基覆岩层倾角的增加，土石混合体边坡的下、上限安全系数的平均值变化相对较小。与基覆岩层倾角为0º的边坡下、上限安全系数相比，当土石混合体边坡的基覆岩层倾角分别为5º、10º、15º、20º时，其下限安全系数的平均值分别减小1.0%、3.4%、3.1%、1.7%，其上限安全系数除基覆岩层倾角为5º时略有增加外，其余的分别减小2.2%、3.3%、2.7%。从以上分析可见，基覆岩层倾角的不同，对土石混合体边坡安全系数平均值的影响相对较小。

从表4中可以看出，相同的含石率及基覆岩层倾角下，由于块石空间位置的不同，土石混合体边坡的下、上限安全系数呈现出较大的差异性。当含石率为40%时，基覆岩层的倾角从0º增加到20º时，土石混合体边坡下限安全系数的变异系数在1.9%～3.1%之间变化；土石混合体边坡上限安全系数的变异系数在1.3%～4.9%之间变化。可见，对于分析土石混合体边坡的稳定性时，应采用随机的方法而不应该采用确定性方法。

4.2 不同基覆岩层倾角下边坡的剪切耗散图

限于文章的篇幅，选取不同基覆岩层倾角下土石混合体边坡的上限剪切耗散图，每种基覆岩层倾角选取一种块石位置，如图13所示。从图13中可以看出，当土石混合体边坡的基覆岩层倾角在0º～20º变化时，坡体中的剪切带不再是圆弧型，而是呈现出与不同含石率下相同的剪切带扩展模式，分别是“绕石”模式、“分流”模式以及“包含”模式。对于“绕石”模式，剪切带可能沿着块石的两侧发生扩展，剪切带具体的路径决定于块石与剪切带的相对位置。对于“分流”及“包含”模式，由于剪切带上分布较多块石，局部区已形成块石骨架，所以剪切带在扩展过程中将绕过块石骨架，在块石周边和土体相互贯通，最终可形成多个剪切带。本文的研究结果与文献[20]中剪切带扩展趋势的试验结果一致。

5 不同计算方法的对比验证

为验证本文提出的土石混合体边坡稳定性分析方法的可行性，以不同含石率为例，将本文方法的计算结果与采用等效强度参数模型的计算结果进行了对比验证。

5.1 土石混合体等效强度参数模型

目前土石混合体等效强度参数模型主要是Lindquist模型^[21]及Kalender模型^[22]。采用土石混合体等效强度参数模型分析土石混合边坡的稳定性时，将整个坡体看作由等效强度材料构成的均质边坡。因Lindquist模型本身较为粗糙，本文采用Kalender模型进行对比验证，Kalender模型的具体描述为：

${\varphi _{{\rm{SRM}}}} = {\varphi _{{\rm{matrix}}}}\left[ {1 + \frac{{1\;000\left[ {(\alpha /{\varphi _{{\rm{matrix}}}}) - 1} \right]}}{{1\;000 + {5^{((100 - VBP)/15)}}}}} \right](\frac{{VBP}}{{VBP + 1}})$

(10)

$UC{S_{{\rm{SRM}}}} = \frac{{\left(A - A\left(\dfrac{{VBP}}{{100}}\right)\right)}}{{(A - 1)}}UC{S_{{\rm{matrix}}}} ,\; \; \;0.1 \le A \le 500$

(11)

${c_{{\rm{SRM}}}} = \frac{{UC{S_{{\rm{SRM}}}} \times (1 - \sin\;{\varphi _{{\rm{SRM}}}})}}{{2\cos\; {\varphi _{{\rm{SRM}}}}}}$

(12)

式（10）、（11）、（12）中， ${\varphi _{{\rm{SRM}}}}$ 为土石混合体的内摩擦角， ${\varphi _{{\rm{matrix}}}}$ 为土石混合体中土体的内摩擦角， $\alpha $ 为块石的休止角， $UC{S_{{\rm{SRM}}}}$ 为土石混合体的单轴压缩强度， $UC{S_{{\rm{matrix}}}}$ 为土石混合体中土体的单轴压缩强度， $A$ 为由土体单轴压缩强度及块石的休止角共同确定的系数， ${c_{{\rm{SRM}}}}$ 为土石混合体的黏聚力。

根据文献[18, 22]，式（10）～（12）中的 $\alpha $ 取40º， $A$ 取0.25， $UC{S_{{\rm{matrix}}}}$ 取33 kPa。计算完成的等效强度参数如表5所示。土石混合体边坡稳定性分析时重度、弹性模量、泊松比及基岩参数的取值同表2。

5.2 采用等效强度参数模型时的结果分析

采用有限元极限分析法对采用等效强度参数模型下土石混合体边坡的稳定性进行了分析计算，获得了各含石率下土石混合体边坡的下、上限安全系数，计算结果如表6所示。

从表6中可以看出，当含石率低于30%时，此时的块石悬浮于“土体”中，采用Kalender模型计算得到的土石混合体边坡稳定安全系数略小于均质边坡的安全系数。当含石率大于30%时，Kalender模型计算得到的土石混合体边坡稳定安全系数均随着含石率的增加而显著提高，可见采用Kalender模型计算得到的安全系数的变化规律与本文随机块石模型得到的结果变化趋势一致。但是采用等效强度参数计算土石混合体边坡稳定性的这种等效均匀化的方法，不能考虑边坡滑动面的“绕石”等效应。因此，采用Kalender模型分析土石混合体边坡的稳定性时，不能正确地得到边坡实际的滑动面。

5.3 不同计算方法结果的对比验证

本文基于随机块石模型及Kalender模型的土石混合体边坡在不同含石率下的下、上限安全系数结果如图14、15所示。

从图14、15中可以看出，当土石混合体边坡中的含石率大于30%时，随着含石率的增加，两种方法计算所得安全系数的变化趋势一致，都是随着含石率的增加而增加；当含石率小于30%时，两种方法得到的安全系数存在较大差异。为分析不同含石率下3种方法计算结果的一致程度，本文以 $\alpha=5{\rm{\% }}$ 的显著性水平进行T检验，检验各含石率下两种方法计算所得下、上限安全系数的是否存在差异，计算结果如表7、8所示。

从表7、8中可以看出，采用本文的方法与Kalender模型所得的安全系数存在较大差异。因土石混合体边坡中块石空间分布的多样性，采用本文随机块石模型的方法能考虑边坡滑动面的“绕石”等效应，从侧面反映了本文提出的土石混合体边坡稳定性的分析方法是可行性。

6 结　论

根据土石混合体边坡中块石含量、级配、基覆岩层倾角等参数，建立了随机块石模型，最后将模型导入到Optum G2中采用有限元极限分析法对不同块石含量、块石分布及基覆岩层倾角的土石混合体边坡进行了稳定性分析，并将计算结果与采用Kalender模型的计算结果进行了对比验证，主要结论如下：

1）相同基覆岩层倾角下，由于块石空间位置的不同，相同含石率下土石混合体边坡的最大、最小下限安全系数的差值从0.022（ $ VBP$ =10%）变化到0.109（ $ VBP$ =60%）；最大、最小上限安全系数的差值从0.060（ $ VBP$ =10%）变化到0.396（ $ VBP$ =60%）。

2）相同含石率下，随着基覆岩层倾角的增加，土石混合体边坡的下、上限安全系数的平均值变化相对较小。

3）随着土石混合体边坡中含石率、基覆岩层倾角的增加以及块石空间分布位置的不同，坡体中的剪切带呈现出“绕石”、“分流”以及“包含”这3种典型的扩展模式。

4）将本文方法与Kalender模型所得的安全系数进行了对比分析，验证了本文方法分析土石混合体边坡的可行性。

最后需要指出的是：本文仅考虑了土石混合体边坡在自重荷载作用下的稳定性变化，其在降雨及地震作用下的稳定性变化规律需后续进一步研究。