基于数字图像技术的冰水堆积体边坡稳定研究

杨 彬1，王潇弘2，吉 锋1

（1.地质灾害防治与地质环境保护国家重点实验室（成都理工大学），四川成都610059；2.华北水利水电大学，河南郑州450045）

摘 要：冰水堆积体在西部地区广泛存在，具有高度的非均匀性与非连续性。与一般均质土体相比，冰水堆积体在力学行为上存在很大差异，而通常计算中也经常将这一特征忽略，由此造成了较大误差。以野外调查为基础，运用数字图像处理技术对典型冰水堆积体边坡处理后得到细观物理模型，并导入FLAC3D中采用强度折减法计算其安全系数。计算结果表明，受块石存在的影响，边坡内部的应力场分布“异常”，埋深对应力集中有明显的“放大效应”，潜在滑动面呈现多级发展，具有显著“绕石”特征。对比计算表明，冰水堆积体边坡的安全系数较均质边坡有较大提高。

关键词：冰水堆积体边坡；数字图像处理；强度折减法；安全系数

冰水堆积体是冰水消融时冰下径流和冰水前缘水流的堆积体，是第四纪以来（大多属于第四纪更新世早中期）逐渐形成的。其经历了较长时期的自然固结，存在不同程度的泥质或钙质胶结，内部块石大小混杂，级配较为宽广。土石混合物的研究正处于起步阶段，虽然现有岩土体工程分类体系没有对其作出解释，但是在实际的工程中其引发的问题却较为广泛。本文所研究的冰水堆积体也属于土石混合体，在细观结构及力学研究中，可以借助土石混合体的研究成果帮助我们对冰水堆积体的破坏形态、破坏机制进行细致深入的研究。关于冰水堆积体的物理力学性质已有较多的研究，如冯俊德等［1］对云南某铁路冰碛土进行了大型直剪试验，认识到原状冰碛土具有一定结构强度，重塑冰碛土强度会随着含水率的增加而降低；范余敏等［2］对雅泸高速西冲特大桥22＃墩冰水堆积体边坡开展了研究，以大剪试验为基础，计算了不同工况下边坡的稳定性；袁广祥等［3］在研究川藏公路然乌—鲁朗段冰水堆积体的过程中，首先分析了堆积体的成因，进而论述了其结构特征、力学特征，最后采用圆弧法对边坡稳定性进行了分析。

现阶段计算冰水堆积体边坡稳定性大多是将整个边坡视为均质连续体，冰水堆积体为赋予一组等效物理力学参数，进而得出结论。然而，冰水堆积体的细观结构非常特殊，尤其是内部块石，对整个边坡稳定性及破坏方式有着重要影响。

本文以西藏某地区野外调查资料为基础，利用数字图像处理技术及有限差分数值模拟技术对冰水堆积体边坡稳定性进行分析研究。

1 研究区冰水堆积体概况

研究区位于西藏林芝地区，冰水堆积体分布非常广泛，内部块石、土体分选性较差，磨圆程度也较差，大小混杂，粒径大小悬殊。虽然冰水堆积体也属于堆积体，但在强度特征上，冰水堆积体与一般松散堆积体差别很大，其形成于第四纪冰期，岩土体都经过了长时期的压实固结，有些还存在不同程度的胶结，使得冰水堆积体与一般堆积体相比呈现出诸多特点。现场调查的某一冰水堆积体边坡如图1所示。

图1 某冰水堆积体边坡

为了得到此冰水堆积体边坡的粒度分布，并满足统计分析的需求，使测量结果更符合实际情况，在边坡的三个不同地点分别取样，进行颗粒筛分试验，得到三个试样的粒组频率分布，见图2～图4。

图2 1＃试样粒组频率分布

图3 2＃试样粒组频率分布

图4 3＃试样粒组频率分布

从图2～图4可以看出，研究区冰水堆积体内部土体和块石粒径分布具有不均匀性，各粒组占比差异较大。这反映了在整个冰水堆积体形成的过程中，冰水和冰融水受地形、气候等因素的影响，其携带堆积下来的物质出现空间上的分异，呈现局部某粒组富集或缺失，同时也表明冰水堆积体是一种高度非均质的物质。

2 冰水堆积体边坡模型建立

2.1 数字图像处理技术

数字图像处理（Digital Image Processing）是通过计算机对图像进行去除噪声、增强、复原、分割、提取特征等处理的方法和技术，近年来广泛应用于航空、材料、医学及岩土工程等领域。在岩土工程领域，数字图像处理技术也得到了很大的应用发展。数字图像在计算机中每个像素点对应一个整数值［4－9］，用以表示该点的亮度，即灰度。整个图像由相应的具有不同灰度值的像素点阵构成，不同的灰度值分别代表了图像所包含的不同信息，这些像素点对应的各个离散数据成为下一步数字图像处理的基础。在野外获取岩土体的断面照片后，首先对图像进行灰度处理，然后利用MAT⁃LAB去除噪声，并采用灰度拉伸、灰度值均衡化、自适应阈值分割等方法增强图像，最大程度消除外界各种不良因素（如光照、断面平整度）的影响。

将去噪和增强后的灰度图输入到数字图像处理系统并进行相应的二值化处理，这时就可以在新形成的图像中呈现粒径较大块石的形态。

2.2 模型的建立

岩土工程问题由定性到定量转变的重要标志就是物理模型的建立。不同于一般岩体和土体，冰水堆积体这种三元材料（土、石、接触面）具有典型非均质和非连续性的特征，常用的模型建立方法不能将这种差异表达出来，而借助数字图像处理技术有助于问题的对比分析。图5是图1经过数字图像处理后的二值化结构数字图像。

图5 边坡二值化图像

鉴于资金和时间原因，本研究未编制专门的数字处理系统，而是利用上述软件强大的图片处理功能突出所需的关键信息，然后将二值化图像导入COREL⁃DRAW，充分利用CORELDRAW的自由绘图功能实现图像关键信息的提取，即利用多义线对冰水堆积体中的砾石进行任意形状的逼近，以最大程度接近实际［10］。

将获得的线框模型在CAD中面域化，导入大型有限元程序ANSYS并划分网格，最终生成FLAC3D可以识别的文件格式。FLAC3D模型如图6所示，网格划分产生的节点总数为21 732个、单元总数为17 695个。

图6 冰水堆积体边坡计算模型

由图6可以看出，边坡中块石呈现随机分布的特征，与现场筛分试验得到的结果相吻合，保证了模型的非均匀性。根据油新华等［10］的研究成果，在土石交界面生成界面元（INTERFACE）可以模拟拉裂和滑移，保证模型的非连续性。

3 冰水堆积体边坡稳定性分析

将模型的边界条件设定为：底面边界固定约束，左侧边界法向约束，其余为自由表面。为了减小边界效应对数值计算的影响，对模型左侧边界进行了延拓［11］，延拓部分用均质土体代替。模型中块石和土体的计算参数见表1。

表1 块石、土体计算参数

成分堆积密度／（g·cm－3）弹性模量／MPa泊松比μ黏聚力／kPa内摩擦角／（°）块石1.7 2×1040.25 1 000 42土体2.8 40 0.30 40 24

3.1 临界失稳状态判据选择

本次边坡稳定性计算采用整体强度折减法，即针对岩土工程中应用最为广泛的Mohr－Coulomb准则进行计算，并将边坡的安全系数定义为使边坡刚好达到临界失稳状态时，对其抗剪强度进行折减的程度。此外，边坡临界失稳状态目前有3种判据：①塑性区贯通；②计算不收敛；③位移突变。经过分析可知，前两种判据存在较多的主观任意性，误差较大，本文不予采用。此次计算岩土体本构模型采用了Mohr－Coulomb理想弹塑性模型，因此当边坡进入极限状态时［12］，必然是其中一部分岩土体材料相对于另外一部分发生无限制的滑移，随着折减系数的增大，滑坡体内的岩土体有明显的位移增量，而稳定区的位移增量几乎为零。这样，依据位移突变判据就能确定边坡稳定形态与折减系数的定量关系，通过方程拟合快速准确地得到安全系数，且物理意义明确。位移突变判别时，选择坡顶作为监测点并分析位移—折减系数的关系曲线，可得到较合理的安全系数［13］。因此，本文选择边坡坡顶上某些点的位移突变作为临界失稳的判断依据。

3.2 边坡内部应力场特征

为了更好地研究冰水堆积体这种非均质、非连续体与一般均质土体的差别，为模型赋予一组等效物理力学参数，这样就得到了一个均质边坡模型，其边坡形态与非均质模型一致。冰水堆积体边坡和均质边坡在自重应力作用下最大主应力云图见图7、图8。

图7 冰水堆积体边坡最大主应力云图（单位：MPa）

图8 均质边坡最大主应力云图（单位：MPa）

由图7、图8可以看出，两类边坡最大主应力云图在形态上存在很大差异。均质边坡最大主应力分布较为规律，埋深较小时与地面线近似平行，埋深较大时则有变为水平分布的趋势；冰水堆积体边坡因块石的存在而使最大主应力分布“异常”，应力等值线极不规则，块石应力集中现象明显，应力场的结构效应非常显著，边坡内应力的分布（尤其是应力集中）基本受块石的分布与形状控制，而且埋深越大应力集中越显著。笔者认为原因是块石与土体的弹性模量悬殊，在自重荷载下，冰水堆积体内部块石有平移、转动的趋势，而块石与土体弹性不匹配又造成变形差异，埋深越大自重应力越大，对应力异常造成一种“放大效应”，最终在地质历史发展过程中，冰水堆积体边坡经过不断的应力调整形成现今的应力场。冰水堆积体边坡内部这种独特的细观力学特征［14－15］直接影响其宏观变形机制、稳定状态。

3.3 稳定性分析

通过对两类边坡采用整体强度折减法计算，采用FLAC3D内置的FISH语言，自行开发记录工具，提取边坡坡顶3个有代表性点的总位移，绘出折减系数与总位移的关系图，见图9。根据散点图的特征，假定拟合曲线形式为双曲线，且设总位移ζ＝（b＋cK）／（1＋aK），其中：K为折减系数，a、b、c为待定系数。当1＋aK→0且b＋cK≠0时，ζ→∞，边坡产生失稳滑动，从而K＝－1／a，此时的折减系数即为边坡的安全系数。利用MATLAB强大灵活的数据拟合和处理功能，将原散点图拟合后就可以得到拟合曲线的相关参数，见表2。

图9 总位移与折减系数的关系

表2 拟合曲线参数

边坡类别代表点a b c R2冰水堆积物边坡1＃－0.369 2－2.438 1.759 0.980 9 2＃－0.371 6－1.995 1.459 0.977 3 3＃－0.374 6 1.6560－0.226 0.991 4均质边坡1＃－0.470 3－75.220 54.700 0.950 3 2＃－0.487 5－38.210 28.650 0.994 8 3＃－0.482 9－59.800 44.680 0.983 5

从表2可以看出，拟合复相关系数（R2）均在0.95以上，表明双曲线拟合效果好，可信度高。通过计算得到，冰水堆积体边坡的安全系数F＝2.721 6，均质边坡的安全系数F＝2.132 8。可见，在考虑块石存在后，计算得到的边坡安全系数有了较大的提高（27.6%）。

由摩尔库仑强度理论可知，岩土体破坏原因是某一面上的剪应力达到了岩土体的剪切强度，此时，剪切面上必然发生较大的剪切应变［16］，因此可通过搜索最大剪应变增量的位置确定边坡滑动面。基于上述思路，在FLAC3D中绘出冰水堆积体边坡和均质边坡折减系数分别为2.7和2.1时的剪切应变增量云图，可以看出，冰水堆积体边坡的潜在滑动面呈现多级发展，贯穿带表现出“欺软怕硬”的特征，具有明显的绕石现象，这与现场的推剪试验现象（图10）一致。而均质边坡中潜在滑动面则与经典的极限平衡法的一般结论一致，呈圆弧滑动破坏。冰水堆积体边坡这种多滑面、绕石特征增加了抗剪强度，从而提高了其稳定性。

图10 野外原位推剪试验

4 结 论

土石混合体在自然环境中广泛存在，冰水堆积体作为土石混合体中的一种，随着西南地区工程建设的不断推进，已越来越频繁地涉及。本文以野外调查为基础，利用FLAC3D对冰水堆积体边坡稳定性特征进行了探讨，得出以下结论：

（1）对野外冰水堆积体进行三组颗粒筛分试验，结果表明同一边坡不同位置各粒组占比差异较大，反映了冰水堆积体形成过程中，冰水和冰融水受地形、气候影响，所携带物质出现空间分异，局部某粒组富集甚至缺失。

（2）冰水堆积体细观结构的复杂性决定了其力学性质的复杂性。冰水堆积体边坡块石的存在使得最大主应力等值线分布“异常”，边坡内应力的分布（尤其是应力集中）基本受块石的分布与形状控制，而且埋深越大应力集中越显著，呈现“放大效应”。

（3）采用边坡整体强度折减法，以位移突变作为边坡临界失稳判据，通过拟合曲线计算确定了两类边坡的安全系数，冰水堆积体边坡安全系数较均质土边坡有较大提高（27.6%）。从剪切应变增量云图可以看出，冰水堆积体边坡的潜在滑动面呈现多级发展，贯穿带表现出“欺软怕硬”的特征，具有明显的绕石现象。这与现场的推剪试验现象一致，而这恰恰提高了冰水堆积体边坡的稳定性。

参考文献：

［1］ 冯俊德，李建国，汪稔，等.云南某铁路冰碛土大型直剪强度特性试验研究［J］.岩土力学，2008，29（12）：3206－3210.

［2］ 范余敏，黄润秋，裴向军.雅泸高速冰水堆积体室内试验及典型边坡稳定性研究［M］.成都：四川科学技术出版社，2009：35－56.

［3］ 袁广祥，曾庆利，尚彦军，等.川藏公路然乌—鲁朗段冰碛高边坡稳定性分析［J］.地质灾害与环境保护，2007，18（4）：48－51.

［4］ 徐文杰，胡瑞林，岳中琦.基于数字图像处理的土石混合体细观结构研究［J］.辽宁工程技术大学学报，2008，27（1）：51－53.

［5］ 杨传俊.数字图像处理技术在土石混合体颗分中的应用［J］.云南水力发电，2009，25（3）：24－25.

［6］ LEBOURG I，RISS J，PIRAR E.Influence of Mor⁃Phological Characteristics of Heterogeneous Moraine Formations on Their Mechanical Behaviour Using Image and Statistical Analysis［J］.Engineering Geology，2004，73：37－50.

［7］ 徐安权，徐卫亚，石崇.基于数字图像的大型堆积体细观力学特性及力学参数研究［J］.岩土工程学报，2012，34（1）：58－64.

［8］ 廖秋林，李晓，朱万成，等.基于数码图像土石混合体结构建模及其力学结构效应的数值分析［J］.岩石力学与工程学报，2010，29（1）：156－162.

［9］ 冈萨雷斯·伍兹.数字图像处理［M］.3版.北京：电子工业出版社，2010：22.

［10］ 油新华，李晓，贺长俊.土石混合体实测结构模型的自动生成技术［J］.岩土工程界，2003，6（8）：60－62.

［11］ 徐文杰，王玉杰，陈祖煜，等.基于数字图像技术的土石混合体边坡稳定性分析［J］.岩土力学，2008，29（增刊1）：341－346.

［12］ 林杭，曹平，李江腾，等.边坡临界失稳状态的判定标准［J］.煤炭学报，2008，33（6）：643－647.

［13］ 林杭，曹平，宫凤强.位移突变判据中监测点的位置和位移方式分析［J］.岩土工程学报，2007，29（9）：1433－1438.

［14］ 徐鼎平，汪斌，江龙剑，等.冰碛土三轴数值模拟试验方法探讨［J］.岩土力学，2008，29（12）：3466－3470.

［15］ 廖秋林.土石混合体地质成因、结构模型及其力学特性、固流耦合特性研究［D］.北京：中国科学院研究生院，2007：56.

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

【责任编辑 张华岩】

Study on the Stability of Glacial Deposits Slope Based on Digital Image Technology

YANG Bin1，WANG Xiaohong2，JIFeng1

（1.State Key Laboratory of Geohazard Prevention and Geoenvironment Protection，Chengdu University of Technology，Chengdu 610059，China；2.North China University ofWater Resources and Electric Power，Zhengzhou 450045，China）

Abstract：The glacial deposits exists extensively in the west，with high degree of nonuniformity and discontinuity.Contrastwith the general homogeneous soil，there aremany differences inmechanical behavior of glacial deposits，which are frequently ignored in in the general calcu⁃lation and produced a large error.Based on the field investigation，the mesoscopic physical model of glacial deposits slope was obtained through digital image processing，which could be imported to FLAC3Dand safety factor was figured out with strength reduction method.The results indicate that rock blocksmake stress field of slope abnormal，and the burial depth has obviousmagnifying effect on stress concentra⁃tion.Meanwhile，there aremulti⁃level sliding surfaces with obvious zigzag.By the way，the safety factor of glacial deposits slope has been much improved in the contrast.

Key words：glacial deposits slope；digital image processing；strength reduction method；safety factor

中图分类号：TU413

文献标志码：A

doi：10.3969／j.issn.1000－1379.2017.09.022

收稿日期：2016－09－27

基金项目：国家自然科学基金资助项目（51308082）。

作者简介：杨彬（1989—），男，江苏徐州人，硕士研究生，研究方向为岩土体稳定性及工程环境效应。

E⁃mail：1356024040＠qq.com