基于数字图像分析的冰水堆积体结构建模与力学参数研究

石 崇 ，王盛年 ，刘 琳 ，陈鸿杰

（1. 河海大学 岩土力学与堤坝工程教育部重点试验室，南京 210098；2. 河海大学 岩土工程研究所，南京 210098）

1 引 言

冰水堆积体是由第四纪冰川运动引起，经历了后期的斜坡再造运动堆积形成的一类特殊地质体。其主要特点[1－2]：①物质成分为混合介质，即由软弱的“土”和坚硬的“石”组成，一般土为基质，石为填充物；②不同尺度块石分布具有强烈的不均匀性和随机性；③在长久的地质作用过程中基质土有一定的黏结力。因此，冰水堆积体可认为是一种典型的土石混合介质，其力学特性受极不均匀的土石颗粒组成、胶结程度和颗粒物理力学特性的控制，具有参数确定困难、分析方法困难、加固治理困难的“三难”特点。特别在高烈度地震、强降雨等影响下很容易产生地质灾害，其稳定性问题是坝址选择、渠道或公路选线、工程治理及地质灾害评估的重点。

近年来，对冰水堆积体的研究主要集中在现场试验和数值模拟两方面，并形成了一些统一认识：冰水堆积体受物质组成、结构分布、含石量、颗粒级配、排列方式等因素的影响，其力学性质极其复杂；传统的物理试验和数值模拟都难以对介质的细观结构特征进行描述，使研究成果存在很大的制约；因此，采用新的思路解释介质的细观结构，才能从本质上揭示其力学特性和变形破坏规律。Holtz等[3]于 1956年就发现，只有满足尺寸与颗粒级配的要求，类似于冰水堆积体的土石混合介质试验才能够反映实际情况；后来 Chandler[4]进一步指出，当试样中含有异常大砾石时强度值会大幅提高，但该值并不是该介质的准确强度；此后，Lanaro等[5]通过使用激光扫描技术对砾石土构成的多元混合介质进行扫描，得到砾石的三维图像，并且对扫描结果分别用傅里叶以及几何分析方法进行了分析，从而得到了土石结构的分布特征；李晓等[6]则基于野外土石混合体结构统计获取土石分布特征的岩土力学试验表明，冰水堆积土石混合体是非均质、非连续体，其力学性能主要受控于内部结构。此外，一些现场及室内试验[7－9]也表明，土石混合介质的应力-应变关系为非线性硬化型，其关系比较符合邓肯-张模型的双曲线假设，而强度则与所含石块的组成、分布密切相关。

但限于试验条件，冰水堆积体室内试验一般剔除了大尺寸的块石，由此得到的介质力学参数一般小于含块石原位介质，由此进行的稳定性分析往往与现场条件存在偏差。如何考虑大尺寸块石对冰水堆积体的影响成为冰水堆积体性质认识的关键。近年来数字图像处理技术[10－11]在岩土工程领域的广泛应用丰富了堆积体细观结构的研究手段，利用数字成像辨别“土石”的方法[12－15]已有学者进行了探讨，并取得了卓有成效的成果。在进行力学性质研究时，已经摒弃了模拟连续介质的方法，转而采用非连续的颗粒离散元方法[16]、DDA方法等，但根据图像转化为地质力学模型时，一般是依据像素位置采用隔层错位进行构建细观模型，该方法导致纵轴方向被压缩，与事实不符，得到的力学规律误差较大。

本文基于土石结构特征，对冰水堆积体混合介质的细观建模提出了颗粒组合构建法，使得颗粒受力特征更为合理，并采用该方法与颗粒离散元方法针对云南古水水电站枢纽区坝前堆积体工程，对介质的宏观力学参数进行了探讨；可为全面认识冰水堆积体混合介质的力学特性和岩土力学参数提供理论依据。

2 数字图像结构特征

研究冰水堆积体介质，首先必须获得介质中土石结构分布特征，以建立介质分布模型。常规的方法是在现场采用统计窗的方法进行，耗时费力。借助图像处理技术可使得冰水堆积体细观结构统计简单易行。其基本步骤如图1所示，首先通过现场拍摄的图片作为对象，利用数字图像处理技术分辨出“土”、“石”结构特征；其次，根据该特征进行力学建模；最后，利用该模型进行力学分析。

数字图像处理是指将图像信号转换为数字信号，并利用计算机对数字信息进行处理的过程。数字图像在计算机中由一个个像素点构成，而数码图片颜色空间常采用 RGB格式，通常以一个三维矩阵对其信息进行存储，可用一离散函数f(i,j,k)来描述图像矩阵信息，其中，(i,j)代表像素点在该图像对应像素点矩阵中的行、列号，k代表通道信息，彩色图像k=1～3，分别代表R、G、B通道。因此，通过命令 R=f (:,:,1)即可获取图像的R通道信息。R通道信息以一个二维矩阵进行存储，可以用式（1）表示，其中任一值的范围为[0～255]。

式中：min(R,G,B)表示每个像素中R、G、B的最小值。

常用的颜色空间很多，如RGB、CMYK、HSI、HSV、HSL等。其中，HSI色彩空间是用色调(Hue)、色饱和度(Saturation)和亮度(Intensity)来描述色彩。HSI对于人的视觉特性较 RGB的优越性使之经常被用于人的视觉系统。在图像处理和计算机视觉中也经常采用HSI色彩空间，这大大简化了图像分析和处理的工作量。

HSI和 RGB色彩空间只是同一物理量的不同表示法，因而它们之间存在着转换关系。按照式（2）～（4）可将RGB色彩空间转化为HSI色彩空间。

图1 基于图像检测的冰水堆积体细观结构统计Fig.1 Structure statistics of outwash deposits microstructure based on image detection

在数字图像数字化后，根据色调分布与饱和度分布特征找到土、石分布的分界值。该方法得到的结果往往是带有损伤的，需要对照片进行前处理，如增加对比度、光照等措施才可以获得较好的效果。

图像识别是一种有损识别方法，由于光照、阴影、拍照角度的差异每一幅图片中土、石区分的阈值均不一样，其不可能完全识别出土、石。为了克服这一缺陷，本文首先在图像识别后，将有损的土石照片导入AUTOCAD，按照同等比例与原照片进行处理（见图2），然后对照原照片连接每一块石的边界。由于加入了人工识别的因素，对部分不合理的边界可进行人工修正。

图2 基于图像识别冰水堆积体结构建模方法Fig.2 Modeling of outwash deposits based on image recognition

3 模型实现及参数选取

3.1 颗粒组合细观模型

图像信息经数字图像处理系统处理后可获取堆积体“土”、“石”介质的相关信息，在将二元化的介质由像素转化为力学模型时，常见的处理方法如图3所示。

图3 常见的基于像素的颗粒组合方式Fig.3 Common combinations of particles based on pixel

如图 3(a)模型，像素直接对应颗粒位置，但由于各颗粒半径间的对应关系，在承受垂直向荷载时，水平向变形接近 0，因此，模型在双轴压缩试验时泊松比接近 0，难以模拟土、石的变形力学特性，一般不采用；图3(b)采用偶数排像素移位，使得颗粒间成60°角，该方法导致纵轴（垂直向）被压缩，使得数字信息中的石块尺寸失真，土石结构效应不准确；图 3(c)双球组合模型，该方法不会造成纵轴压缩，但像素代表的实体面积不相等，且边界上小颗粒处于应力奇点，容易造成结果不收敛，这与事实不符。

为了克服这些方法的缺陷，本文建立了颗粒组合构建方法：首先在选定的研究范围内，采用平均颗粒尺寸4.5 mm，最大、最小半径比2.0生成随机分布颗粒，如图2所示模型在0.5 m×1.0 m范围内共生成颗粒12 483个；通过PFC2D软件自带fish函数调整颗粒接近零应力状态；然后，在前述图片识别基础上，将识别出来的块石视作不同的多边形区域(可为凹多边形)，搜索所有颗粒，若某颗粒中心位于其中一个多边形内部，则判断该颗粒属于岩石颗粒，位于同一多边形域内的颗粒通过clump组装在一起模拟岩石介质；最后，为了能模拟块石间的接触，不同多边形域的颗粒赋予不同的编号，土、石分别赋予不同的参数以分别模拟“土”、“石”性质。

3.2 土石颗粒宏、细观力学参数

PFC2D采用圆盘(或圆球)形颗粒，从散体介质的细观力学特性出发，采用显式差分算法，在计算循环中交替运用力–位移定律和牛顿运动定律，根据力–位移定律更新接触颗粒间的接触力；通过运动定律更新颗粒的位置，构成颗粒的新接触，从而模拟颗粒介质的运动及其相互作用过程。

PFC2D的接触本构模型包括滑动模型、接触刚度模型和平行黏结模型。本文土体采用接触刚度模型，它适用于以颗粒接触为主的散体介质模拟。该模型有5个参数需要确定：法向刚度kn、切向刚度ks、摩擦系数f、法向黏结强度、切向黏结强度，这些参数需根据土、石室内试验，通过不断调整细观参数，使应力-应变曲线与室内试验曲线吻合得到。

表1 宏观岩土力学参数Table 1 Mechanical characteristics of macro medium

图4 土体数值试验与室内试验应力-应变曲线对比Fig.4 Comparison of stress-strain curves between simulation test and laboratory test

云南省某水电站坝前冰水堆积体由砂土构成，充填碎块石，碎块石为变质砂岩、板岩。宏观岩土力学参数如表1所示。由于土体的直接剪切试验、大三轴压缩试验仪器的限制，在进行岩土力学试验时剔除了粒径大于2 cm的碎块石。本文采用的图像识别方法正是为了评价这些大于2 cm石块的影响，因此，土、石是一个相对的概念，本文中定义“土”为颗粒小于2 cm的散体，“石”为大于2 cm堆积体块石。

通过不断调整土颗粒参数（含石率 0%）细观参数，使其逼近室内试验应力-应变曲线，反演出土细观力学参数如表2所示，试验及模拟试验曲线对比如图4所示。在模型加载过程中，采用应变控制加载方式，利用伺服控制程序保持围压恒定。模型的上、下边界作为加载板，加载速率为5 mm/s，当轴向应变达到峰值强度后再继续执行 5%应变时停止，综合峰值强度与围压对应关系定出强度包络线。

表2 土石颗粒细观力学参数Table 2 Mesoscopic parameters of soil and rock particles

根据现场数据共进行了24幅数字图像分析，结果表明，该冰水堆积体介质中超过2 cm的块石含量在20%～40%之间，定义块石边界任意点最远距离为长轴，与长轴垂直方向点最大距离为短轴。进行几何统计，发现块石长轴在4.1～15.4 cm，均值和方差分别为10.3、3.3 cm；轴比（短轴与长轴比值）的均值为0.33，方差为0.083；长轴方向与水平向夹角均值为94.5°，方差为 50.0°。同时，定义块石所占图片面积占总面积的比率为含石率，结果表明，24组图片中，含石率（粒径大于2 cm）有一组约10%，其余在 20%～40%之间，为了评价土石结构组成及块石含量对冰水堆积体的影响，对每个模型进行0.1、0.3、0.5 MPa 3个围压下的双轴压缩数值试验，研究其变形性能及强度特性变化。

4 试验结果分析

4.1 应力变形分析

从图5可以看出，冰水堆积体双轴压缩时应力-应变过程经历了4个阶段。第1阶段为初始压缩阶段，由于土、石的刚度差异巨大，此时的荷载主要由边界上的块石颗粒承担，导致其变形模量较高，变形曲线近似为线性；由于土颗粒强度低，很快进入屈服，压缩变形进入非线性段，土体被压缩，在块石与土体接触部位产生应力集中（见图6），土、石进入了协调变形阶段，此时的应力-应变曲线为非线性硬化型，更接近于土体；当堆积体达到峰值强度后，继续变形导致土体的黏结强度逐步丧失，承载力下降，进入应变软化段；当土体颗粒黏结力完全消失后，此时堆积体的强度主要受土颗粒的摩擦控制，但仍具有一定的残余强度。

在变形过程中，由于土体模量低、变形快，不均匀变形导致块石承受较强的拉应力集中（见图6），块石边界受力不均可引起几何旋转，使得土石咬合作用改变。在块石周边应力较为集中的部位土体首先破坏，因此，块石分布会影响土体的破坏路径。

图5 典型介质变形特征曲线Fig.5 Curves of typical deformation characteristic

图6 压缩截面I-I上正应力分布Fig.6 Normal stress distribution on section I-I

图7 压缩区域内轴向应力分布Fig.7 Axial stress distribution in compressive area

如图7所示不同截面处的轴向应力随块石分布变化剧烈，可导致如图6所示压、拉应力的急剧转变，因此，不同含石率、几何结构组成、块石空间分布的块石对冰水堆积体变形与强度特性的影响不同，对混合介质的破坏路径也有重要影响，这是造成冰水堆积体岩土力学参数离散性较大的主要原因。

4.2 土体黏结强度影响分析

随土颗粒黏结强度的提高，冰水堆积体混合介质的变形模量升高，相应峰值抗压强度增大（见图8）。表明土体颗粒间的黏结强度与冰水堆积体的峰值强度和变形模量呈正相关变化。

图8 不同颗粒黏结力下双轴压缩曲线Fig.8 Biaxial compression curves under different adhesion forces

但土颗粒黏结力的提高，对冰水堆积体的抗剪切强度指标并非同比增加，随着黏结力的提高，土石黏结后的抗剪内摩擦角有下降的趋势（见图9），而黏聚力上升（见图 10）。这表明由于块石分布引起冰水堆积体的强度变化主要由黏聚力引起，摩擦角的下降表明，黏结在一起的冰水堆积体更容易产生脆性破坏。

图9 内摩擦角随土体黏结力变化Fig.9 Friction angle changes with soil cohesive force

图10 黏聚力随土体黏结力变化Fig.10 Cohesion changes with soil adhesion

由以上不同细观参数下的应力-变形规律可知，冰水堆积体的力学特性受基质土颗粒属性影响明显。在漫长的地质作用下，若土颗粒间胶结程度较好，则堆积体峰值强度远高于残余强度，其属性更接近于弹塑性介质；若胶结程度较差，则堆积体属性更接近于非线性硬化土体，大尺寸块石存在可改变变形破坏路径，使得混合介质抗剪强度提高。

4.3 土石混合介质宏观力学参数分析

在不同块石体积含量情况下，冰水堆积体在稳定变形阶段受基质材料的影响较为明显，在块石体积含量少于25%时，此时基质为土、块石含量的变化对堆积体的抗剪切强度影响不明显，均接近于土体强度（见图11）；当块石超过25%时，随着块石含量的提高，堆积体内摩擦角提高，其幅度为块石含量每提高 1%摩擦角上升 0.118°，黏聚力提高6.1 kPa。

根据这一关系可近似拟合出该处冰水堆积体宏观强度与含石率的关系（见图12、13）：

图11 不同块石体积含量下典型特征曲线Fig.11 Biaxial compression curves under different stone contents

图12 块石体积含量与内摩擦角关系Fig.12 Relation between stone content and friction angle

式中：φ为内摩擦角（°）；x为块石百分含量（%）；c为黏聚力（kPa）。式（5）、（6）可为不同含石率的堆积体工程抗剪强度预测提供依据。

图13 含石率与黏聚力关系Fig.13 Relation between stone content and cohesion force

相关性研究发现，内摩擦角与含石率相关系数为0.46；黏聚力与块石体积含量的相关系数为0.79。这表明黏聚力与含石率的相关性要大于内摩擦角，即使含石率接近，内摩擦角也可能差距甚大，这是因为冰水堆积体介质变形过程中，块石的分布结构影响也非常明显造成的。

通过式(5)、（6）判断，考虑块石影响下的冰水堆积体介质抗剪强度要比纯土体试验获得的参数高0°～8°，平均为 3.7°。

对该工程堆积体边坡进行稳定性计算发现，采用纯土体抗剪参数得到的安全系数小于1.0，与工程判断边坡基本稳定结论不符，而通过式(5)、（6）提高摩擦角 3.7°后稳定性略大于 1.0，与工程判断吻合，表明借助颗粒流分析方法对介质含石影响研究具有一定可信性。

5 结 论

（1）基于数字图像分析确定边界，将随机生成的颗粒分别界定“土”、“石”，可克服常规像素对应生成颗粒导致的缺陷，为土石混合介质的细观建模提供了一种更为合理的途径。

（2）针对云南古水水电站坝前冰水堆积体统计发现，块石长轴尺寸，短轴/长轴比值统计规律较好，而长轴方向角方差大，可认为是随机分布的。由于土石混合介质的特性，其变形曲线在峰前为非线性硬化型，峰后则表现为塑性特性。

（3）冰水堆积体的强度受基质土体性质的影响，尤其是在块石含量小于25%时，其强度接近基质土；而当土颗粒的黏结力较大时，混合介质内摩擦角随黏聚力增加有下降趋势，表明由土体与块石黏结在一起的堆积体更显脆性。

（4）冰水堆积体的宏观抗剪切强度受块石含量、结构分布影响明显，24组数值试验表明，工程堆积体介质内摩擦角平均为29.85°，比均质土体的26.15°高3.70°，对强度指标提高明显，该数据对合理估计天然介质的抗剪强度提供了依据。

（5）基于现场数字图像识别的冰水堆积体“土”、“石”识别方法可快速有效地将现场图像与力学模型联系起来，并可弥补室内试验不能考虑大尺寸块石的缺陷，是冰水堆积体介质岩土力学参数获取途径的有效辅助。

[1]韩金良, 吴树仁, 李东林, 等. 秦巴地区地质灾害的分布规律与成因[J]. 地质科技情报, 2007, 26(1): 101－108.HAN Jin-liang, WU Shu-ren, LI Dong-lin, et al.Distribution regularities and contributing factor of geological hazards in Qinling-Daba mountains[J].Geological Science and Technology Information, 2007,26(1): 101－108.

[2]于永贵. 三峡库区万州-开县段堆积层滑坡发育规律及破坏模式[D]. 北京: 中国地质大学, 2008.

[3]HOLTZ W G, GIBBS H J. Triaxial shear tests on pervious gravelly soils[J]. Journal of the Soil Mechanics and Foundations Division, 1956, 82: 1－19.

[4]CHANDLER R J. The inclination of talus, arctic talus terraces, and other slopes composed of granular materials[J]. Journal of Geology, 1973, 81(1): 1－14.

[5]LANARO F, TOLPPANEN P. 3D characterization of coarse aggregates[J]. Engineering Geology, 2002, 65(1):17－30.

[6]LI Xiao, LIAO Qiu-lin, HE Jian-ming. In-situ tests and stochastic structural model of rock and soil aggregate in the Three Gorges reservoir area[J]. International Journal of Rock Mechanics and Mining Sciences, 2004,41(3): 494－450.

[7]伍明. 土石混合非均质填料力学特性试验研究[J]. 公路, 1997, (1): 40－49.WU Ming. Research of mechanical characteristic test on nonhomogenous filling mixed soil and rock[J]. Highway,1997, (1): 40－49.

[8]赫建明. 三峡库区土石混合体的变形与破坏机理研究[D]. 北京: 中国矿业大学, 2004.

[9]徐文杰, 胡瑞林, 曾如意. 水下土石混合体的原位大型水平推剪试验研究[J]. 岩土工程学报, 2006, 28(7): 814－818.XU Wen-jie, HU Rui-lin, ZENG Ru-yi. Research on horizontal push-shear in-situ test of subwater soil-rock mixture[J]. Chinese Journal of Geotechnical Engineering, 2006, 28(7): 814－818.

[10]吴国平. 数字图像处理原理[M]. 北京: 中国地质大学出版社, 2007.

[11]廖秋林, 李晓, 朱万成, 等. 基于数码图像土石混合体结构建模及其力学结构效应的数值分析[J]. 岩石力学与工程学报, 2010, 29(1): 155－162.LIAO Qiu-lin, LI Xiao, ZHU Wan-cheng, et al. Structure model construction of rock and soil aggregate based on digital image technology and its numerical stimulation on mechanical structure effects[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(1): 155－162.

[12]YUE Z Q, BEKKING W, MORIN I. Application of digital image processing to quantitative study of asphalt concrete microstructure[C]//Transportation Research Record 1492. Washington D C: National Research Council, 1995.

[13]LEBOURG T, RISS J, PIRARD E. Influence of morphological characteristics of heterogeneous moraine formations on their mechanical behaviour using image and statistical analysis[J]. Engineering Geology, 2004,73(1/2): 37－50.

[14]REID T R, HARRISON J P. A semi-automated methodology for discontinuity trace detection in digital images of rock mass exposures[J]. International Journal of Rock Mechanics and Mining Sciences, 2000, 37(7):1073－1089.

[15]徐文杰, 胡瑞林, 岳中琦, 等. 土石混合体细观结构及力学特性数值模拟研究[J]. 岩石力学与工程学报, 2007,26(2): 301－311.XU Wen-jie, HU Rui-lin, YUE Zhong-qi, et al.Mesostructural character and numerical simulation of mechanical properties of soil-rock mixtures[J]. Chinese Journal of Rock Mechanics and Engineering, 2007,26(2): 301－311.

[16]CUNDALL P A. A computer model for simulating progressive, large scale movements in blocky rock systems[C]//Proceedings of the ISRM Symposium on Rock Fracture. Rotterdam: A. A. Balkama, 1971: 129－136.