风化花岗岩孤石边坡的稳定性和破坏机理研究

2020-01-18刘晓华韦彬唐皓戴智颖

铁道科学与工程学报 2019年12期

刘晓华，韦彬，唐皓，戴智颖

(深圳市综合交通设计研究院，广东深圳518003)

广东省南部地区的燕山期花岗岩分布十分广泛，由于花岗岩抗风化能力的差异，在构造运动及化学风化等因素作用下，在花岗岩风化带(残积层、全风化层和强风化层)中发育球状风化体(俗称风化球)是该地区一种较为常见的地质现象[1]。由于特殊的地质构造，花岗岩孤石在深圳地区分布较为丰富[2]。花岗岩孤石作为常见的地质灾害，具有分布不确定性、埋藏边界的隐蔽性等特点，在施工的过程中，会带来很多干扰和困难。且花岗岩孤石容易形成落石等自然灾害。会给我们带来极大的生命与财产损失。为了减少生命财产损失，对花岗岩孤石的研究很有必要。刘治军等[3]对花岗岩孤石的3 种存在形态进行了探讨，并提出了不同形态下的稳定性理论计算方法，能够为自然界花岗岩孤石的评估提供依据。杨智翔等[4]分别采用理论计算方法和Rockfall 数值模拟对孤石的稳定性计算进行了研究，能够反映真实的情况。杜甫志[5]对有花岗岩孤石存在的边坡稳定性进行研究，并提出了有效的孤石施工应对措施。本文基于某市洲石路K1+880 左侧边坡塌方的实际工程案例，通过对现场情况进行分析，建立有孤石存在的边坡模型，研究孤石对边坡稳定性的影响，并分别就2 种不同的孤石存在模式采用FLAC3D进行计算分析。

1 工程背景与孤石数值化模型获取

1.1 工程背景

深圳市洲石路K1+880 左侧边坡经过持续暴雨作用，2008年6月18日出现滑坡。本边坡原始地貌为由花岗岩组成的低丘陵。岩体风化极不均匀，山坡上发育大量大小不一的风化球(孤石)，山体坡顶标高为67.14 m，自然山体坡度约为45°。边坡中段K1+880～K1+895 浆砌片石挡土墙顶，至截水沟之间发生浅层小型滑坡，平面形态呈楔形，滑体长14 m，宽15 m，面积210 m2，体积约400 m3，后缘陡坎呈近似直线，陡坎高2～4 m，坡度约80°，为原生结构面，浆砌片石挡土墙未见明显变形。滑坡物质主要为坡残积砾质黏性土及全风化花岗岩，呈松散状，黄褐色，含风化球体(孤石)，大者达2×3 m。

1.2 数字图像处理获取孤石数字化模型

基于MATLAB 平台，通过编写程序对现场获取的花岗岩孤石块体进行轮廓提取，并读取孤石块体轮廓点的坐标。具体实现过程如下。

为了方便检测块石的边缘，我们使用灰度转换函数将自然彩色图像(图1(a))转换为灰度图像(图1(b))。

式中：g(x,y)是输出的灰度图像；f(x,y)是输入彩色图像；T是f(x,y)上的算子；R(x,y)，G(x,y)和B(x,y)是彩色图像的RGB分量。

然后对图像进行高斯模糊处理[6]，以减少图像噪声并降低细节层次。高斯模糊是一种图像模糊滤波器，它用正态分布计算图像中每个像素的变换。所谓“模糊”，可以理解成每一个像素都取周边像素的平均值，是一种简单的平滑处理方法。这种模糊技术生成的图像，其视觉效果就像是经过一个半透明屏幕在观察图像，产生高斯模糊效果(图2(c))。

一维高斯函数的方程为：

根据一维高斯函数，可以推导得到二维高斯模糊方程：

式中：x是水平轴上任意一点到坐标原点的距离；y是竖直轴上任意一点到坐标原点的距离；σ是高斯分布的标准差。

为了检测孤石的外轮廓，我们对图像进行二值化处理(图1(d))，将图像中的孤石和背景分离，孤石区域为黑色(像素值为0)，背景区域为白色(像素值为1)。其中的关键是确定一个合理的阈值进行二值化图像的划分[7]，因此我们选取分离系数η(k)用于评价阈值是否合适：

式中：N是像素点总数；L是灰度级数；Pi是密度函数；ni是第i级灰度内的像素点总数。我们首先计算k从1 到L内所有的σB2(k)，σB2(k)的最大值即为二值化划分的阈值。

通过领域搜索方法即可在二值化的像素矩阵中搜索出花岗岩孤石块体的轮廓坐标(图1(e))。

图1 数字图像处理提取块体轮廓Fig.1 Extraction of rock contour by digital image processing

2 花岗岩孤石赋存的边坡模型构建

2.1 边坡数值模型构建

根据第2 节获取的花岗岩数字模型，本文构建了3 个不同的边坡模型：无花岗岩存在的边坡(图2(a))，花岗岩完全嵌入型的边坡(图2(b))，花岗岩部分露出的边坡(图2(c))。3 个边坡的尺寸均如图2(a)所示。构建这3 个边坡的目的是研究花岗岩孤石对边坡稳定性的影响，以及不同孤石位置对边坡稳定性的影响。

图2 边坡模型构建Fig.2 Construction of slope model

2.2 网格划分

基于DISTMESH2D[8-9]来对3 个不同的边坡进行网格划分。DISTMESH2D 是一种针对隐式函数指定的几何体生成非结构化网格的新技术。通过求解单元边缘的力平衡，迭代地改进初始网格，并使用隐式几何定义投影边界节点。该方法通常可以生成高质量的三角形网格。

在使用DISTMESH 生成网格时，需要对网格的尺寸进行定义。在有限差分方法中，创建高质量的网格非常重要。数值模拟的计算速度通常对网格的数量十分敏感，在FLAC3D中，一般弹塑性问题的求解时间大致与N4/3成正比(N是单元的数量)[10]。网格尺寸太小，单元数量较为庞大，计算耗时较长；网格尺寸太大，计算精度又难以保证，因此，兼顾效率与精度，我们需要通过尝试定义一个合适的网格尺寸。

DISTMESH 允许我们定义某些点来作为固定点，这些固定点必定会是网格的节点。这样做的好处是：我们可以将模型的边界点以及块石的轮廓点固定起来作为网格的节点，以确保土石界面上的点是网格的节点，防止模型进行网格划分后出现变形。在进行网格生成时，为了确保图示界面上原有的两点之间不会有新的节点插入，我们需要保证土石界面上任意两点之间的距离小于网格尺寸。若土石界面上两点之间的距离大于网格尺寸，则将该线段等分成若干小于网格尺寸的线段。

采用DISTMESH 对3 种结构土石混合体模型进行网格划分后的效果如图3所示。

图3 边坡网格划分Fig.3 Slope meshing

网格划分后，通过编写接口程序将边坡模型导入FLAC3D。同时，为了能够在FLAC3D中进行计算，在接口程序中我们将边坡沿垂直纸面方向拉伸1 mm 的宽度，即形成一个单位厚度的三维模型。

3 含花岗岩孤石的边坡稳定性分析

采用FLAC3D内嵌的有限差分强度折减法对边坡进行稳定性分析。

3.1 有限差分强度折减计算理论

3.1.1 基本原理

强度折减法中边坡稳定的安全系数定义为：使边坡刚好达到临界破坏状态时，对岩、土体的抗剪强度进行折减的程度，即定义安全系数为岩土体的实际抗剪强度与临界破坏时的折减后剪切强度的比值。强度折减法的要点是利用式(6)～(7)来调整岩土体的强度指标c和φ然后对边坡稳定性进行数值分析，通过不断地增加折减系数，反复计算，直至其达到临界破坏，此时得到的折减系数即为安全系数Fs。

式中：cF为折减后的黏接力；φF为折减后的摩擦角；Ftrial为折减系数。

3.1.2 实现过程

随着数值计算理论的发展，强度折减法在岩土工程实践中的可行性逐渐被得到了证实。目前存在的困难为尚无统一的边坡失稳判据，即无统一安全系数数值求解过程的终止条件。现行的边坡失稳判据主要有以下几种：

1)以数值计算的收敛性作为失稳判据[11]；

2)以特征部位位移的突变性作为失稳判据[12]；

3)以塑性区的贯通性作为失稳判据[13]。

在FLAC3D中求解安全系数时，单次安全系数的计算过程主要采用的是第(1)种失稳判据。FLAC3D主要以力的不平衡比率小于10-3作为终止条件[13]。

为了控制整个强度折减法循环计算过程的求解时间，规定在计算过程中，只要满足上述3 个标准中的任何一个，便退出当前计算。

3.2 边坡稳定性分析

3 个不同边坡稳定性分析结果如下：

图4 边坡稳定性计算结果Fig.4 Calculation results of slope stability

计算所用的材料参数如表1所示[14-15]。

表1 计算参数Table 1 List of calculation parameter

图4可以看出，均质土体边坡最危险潜在滑动面是类圆弧滑动带；当有花岗岩孤石完全嵌入边坡时，边坡最危险潜在滑动面会绕开孤石呈现出明显绕石效应，且滑动发生在边坡后缘[14-15]；当花岗岩孤石部分露出边坡时，边坡的滑动带依然会绕开孤石，但由于孤石重力作用潜在滑动面靠近边坡前缘。

计算后3 个边坡的安全系数分别为：均质土体安全系数：1.08；孤石完全嵌入的边坡安全系数：1.18；孤石部分露出的边坡安全系数：0.83。对比破坏模式图我们发现，当孤石相对边坡的尺寸较大，且孤石完全嵌入边坡中时，受孤石的影响滑动发生在边坡后缘，相比均质边坡安全系数会增高。而当孤石部分露出边坡时，滑动发生在边坡前缘，且受孤石自身重力的影响，边坡更易发生滑动，因此边坡安全系数会降低。