土体表面干缩裂隙的形态参数定量分析方法

，

(1.三峡地区地质灾害与生态环境湖北省协同创新中心， 湖北 宜昌 443002；2.上海交通大学 土木工程系，上海 200240)

1 研究背景

土体干缩裂隙的产生是一种常见的自然现象，其影响范围涉及到农业土壤学、岩土工程、环境保护等众多领域。由于外界气候交替作用，边坡工程，尤其是膨胀土边坡的裂隙不断发育，一方面导致土体的强度不断衰减，另一方面直接增加了土体的渗透性，两者相互影响，大大提高了裂隙发育范围内浅层性滑坡的发生机率。裂隙的长度、宽度，裂隙网络分布规律及分形特征等形态参数的变化都会对土体工程特性产生重要影响。研究土体干缩裂缝，对其形态参数做出定量分析、描述，对提高边坡稳定性、控制环境污染(尤其是土壤污染)、提高农业用水利用率等均具有重要的科学意义和实践价值[1-4]。

早期对于岩土体裂隙形态参数的获取[5]，主要依靠现场手工测量，效率较低、误差较大；Vogel 等[6-7]借助于拓扑学的概念和方法，建立了裂缝面密度、长度密度、欧拉数等用于描述裂缝大小和连通性的指标体系，对裂缝的定量化做了有益的贡献，但过于简单、抽象。CT法[8]也可以动态、定量地量测岩土体微细观结构的发展变化，但设备过于昂贵，不便推广。

新兴发展的数字图像识别分析技术为实现高效的裂隙测量提供了一种可能途径。数字图像分析通过对裂隙图像的计算机识别，进而得到各种裂隙形态参数。

为了实现对平面裂隙图像的无损、连续、动态识别和形态参数分析，本文借助数字图像处理技术，基于MATLAB软件，开发了裂隙网络图像分析系统。通过该系统对黏土裂隙的识别，可得到蒸发过程中裂隙的裂隙率、裂隙长度、平均宽度，以及隙宽分布规律和分维等形态参数。

2 平面裂隙图像的识别

本程序主要包括：图像二值化处理程序、二值图像桥接去杂处理程序、骨架提取程序，二值图像特征参数提取程序以及图形用户界面设计处理程序。

2.1 图像二值化处理程序

前处理程序主要内容为，首先将图片裁剪成统一大小，通过imread(‘.jpg’)函数读取裂隙图片，然后利用rgb2gray(RGB)函数将裂隙真彩图像转化灰度图像，再利用imhist(gray)函数得到灰度图像的像素直方图，最后根据直方图特点选择适当阈值，一般情况下直方图上有一个峰，将阈值设在峰值前的谷底处效果较好。如图像中裂隙较多且分散，裂隙与背景对比不明显时，裂隙和背景会各对应一个峰值，此时在两峰之间选择阈值效果较好。利用im2bw(gray，level)函数将灰度图像转换成二值图像。图像二值化处理程序流程如图1所示。

图1 裂隙图像二值化处理程序流程图

2.2 二值图像桥接去杂处理程序

由于受到拍摄过程中的光线变化，以及图像噪声的影响，经过二值化程序处理后的裂隙图像中包含大量的杂点，使得原本连续的裂隙可能会断成几段。在进行裂隙特征提取过程中，程序会将某些本是裂隙的视为非裂隙或将本不是裂隙的视为裂隙，从而造成系统误差，因此必须对裂隙的二值图像进行桥接和去杂处理。

桥接过程运用逻辑计算中的膨胀运算和腐蚀运算。这种运算给出一个裂隙图像集合A和一个结构元素集合B(俗称刷子，其大小比待处理的图像要小得多，一般可取为3×3矩阵)，利用结构元素对裂隙图像进行操作。

由图3(a)可知：在同一地震烈度下,从Ⅰ类场地到Ⅳ类场地输电塔地线挂点位移依次增大;在同一场地条件下,从6度到9度地震烈度输电塔地线挂点位移依次增大.因此,对于该钢管塔,最严重的工况组合是Ⅳ类场地和9度地震动.

集合论中，膨胀运算被定义为：A⊕B={x|(B^)x∩A≠Ø}，B^为B关于原点对称的集合，作为B的映射。膨胀运算的具体操作为：用膨胀结构元素来扫描原始图像中的各个像素，将两者做“与”操作，即如果像素都为0(非裂隙像素)，则原始图像对应的像素取为0，否则取为1(裂隙像素)。膨胀运算填补了图像中的一些空洞，即桥接那些本是连续而不连续的裂隙，但也使原始图像中的目标被加长、加粗，这时就需要用到腐蚀运算来修正。

腐蚀运算定义为：AΘB={x|(B)x⊆A}，即将B平移x, 平移后的B仍在A中。腐蚀运算与膨胀运算相反，其作用是使图像中的目标收缩、细化，可以消除裂隙图像中毫无意义的杂点。腐蚀运算的具体操作与膨胀运算一样，只是在最后做“与”操作时有些不同。在腐蚀运算中，如果结构元素与原始图像所对应的像素都为1时，原始图像中对应的像素取为1，否则取为0。

2.3 裂隙骨架提取

图2 图像骨架示意图

骨架广义的定义是用一组与原始形状连通性和拓扑结构分布性相一致的能够完整表达物体形状的曲线集合。通过骨架的定义，我们可以计算得到各种物体图像的骨架图，如图2所示。

裂隙的骨架我们可以定义为到裂隙两边缘距离相等点的集合，即裂隙的中轴。通过骨架化算法可以提取裂隙主干。主干只有一个像素宽，记录了裂隙的主体框架，如图3(a)所示为裂隙二值图，图3(b)为裂隙骨架二值图。

3 土体表面干缩裂隙形态参数

本程序主要对土体裂隙率、裂隙长度，平均宽度，裂隙分维等平面形态参数进行定量化统计。设计流程如图4所示。

3.1 裂隙度(裂隙率)

裂隙长度、裂隙宽度、裂隙深度只能代表裂隙的某单个方面的特征，并不能综合代表裂隙的整体特征，因此，需要一个综合的指标既能反映出裂隙特征，又能体现出裂隙的总体几何分布特征。于是裂隙度就被引入进来，袁俊平等[9]对裂隙度提出了几种概念，但较为复杂。唐朝生等[10]对裂隙度概念做了简化处理，在基于MATLAB图形图像处理中像素统计功能的基础上，提出了裂隙率的概念，即

(1)

式中：δf为裂隙率；Ai为第i条裂隙面积,可通过MATLAB中bwarea(‘裂隙二值图像’)函数得出裂隙图像中裂隙的面积；n为裂隙的条数；A为统计试样总面积；nb为黑色像素点数，即试样中开裂部分；nw为白色像素点数，即试样中未开裂部分；n为总像素点数。

3.2 裂隙长度、宽度、隙宽分布率

图5 裂隙长度计算示意图

裂隙宽度(又称为裂隙开度)反应裂隙的张开程度，其中，裂隙的平均宽度通常被作为衡量裂隙发育程度的一个重要参考指标。本文利用裂隙的面积及长度来推导计算得出裂隙的平均宽度，具体公式为

(2)

式中：d为裂隙的平均宽度；li为第i条裂隙的长度；其他参数含义与上文相同。

考虑到裂隙宽度分布的不均匀性，本文提出隙宽分布率的概念。

提取骨架二值图中骨架坐标以及桥接去杂二值图中裂隙边缘坐标，搜索骨架某点到裂隙边缘的最短距离。如前文所述，裂隙骨架被定义为到裂隙两边缘距离相等点的集合，因此该距离的2倍即可定义为该点处裂隙宽度。本文将介于一定隙宽范围内的裂隙称为隙宽组。隙宽分布率即可定义为各个隙宽组的相对含量(指各隙宽组的裂隙长度占裂隙总长度的百分数)，表达式为

(3)

式中：ci为第i组隙宽分布率；li为第i隙宽组的裂隙长度。

3.3 裂隙分维

分维，即分形维数，应用于量度复杂形体的不规则性，反映复杂形体占有空间的有效性。已有研究表明[11]，裂隙网络分布在一定尺度范围内具有体积分形特征，岩土体裂隙网络的分维值越大，表示风化作用越强烈，裂隙越发育，其空间分布也越复杂。

本文采用盒维数法计算裂隙网络的分维值。具体做法是使用边长为ε的矩形格子(盒子)网络去覆盖裂隙分布区域，存在裂隙几何体的格子数目N(ε)将随着ε而变化，两者关系如方程(4)所示：

lgN(ε)=B-Dlg(ε) 。

(4)

式中：D即为裂隙分维；B为常数。

理论上对任意的ε值，都有一个和其对应的N(ε)，多组(lgε，lgN)即可得到裂隙分维D。但在程序实际操作中需注意到：二值图像由像素矩阵构成，在微观尺度上裂隙边缘表现为锯齿状，对像素进行切分将导致分形维数偏大[12]。因此，在分形维数计算过程中，ε应取整数。另外，在图片裁剪时，所选正方形区域的边长最好控制为一个具有较多约数的数值，如1 260×1 260，以保证能够获得足够多的N-ε组合。

3.4 图形用户界面GUI程序

本部分主要内容为图形用户界面(Graphical User Interface，简称GUI)设计， GUI是用户与计算机程序之间的交互方式，它是包含图形对象，如窗口、工具栏、菜单、文本、图标的用户界面。用户可以以某种方式选择或激活图形对象而运行与其相对应的程序命令。裂隙图像二值化处理、形态参数统计程序界面如图6，图7所示。

图6 二值化处理程序GUI设计

图7 裂隙特征参数提取程序GUI设计

考虑到裂隙演化是一个随时间增长的动态过程，程序可实现对若干裂隙演化图片的同时处理，最终得到的裂隙形态参数将会自动生成为EXCEL文件，便于后期数据处理。

4 应用实例

4.1 裂隙生成

本文设计了一套膨胀土裂隙室内模拟试验，试验装置如图8所示。本次试验所用试样盒的内径为29.5 cm×39.5 cm×10 cm。根据设计的土样干密度和土样厚度计算所需干土(过2 mm筛)的重量，再将土样配成含水均匀的土膏，平铺到试样盒，上覆铁板，利用千斤顶将铁板压到指定的位置，使土膏密度达到试验设计密度后撤去铁板。将试样盒、相机的相对位置调整好，固定在自制支撑架上，土膏在自然风干+风扇轻微加速风干条件下，含水量逐渐降低，裂隙逐渐生成、发展。数码相机与计算机连接，可在指定的任意时间，拍摄土膏的开裂状态。

图8 裂隙演化试验装置简图

试验过程中还需测量土膏的其它参数，为防止在试验称量过程中移动铁盒位置而照成裂隙图片的不统一性，在试验时准备2个大小、含水率等条件完全一致的试样，1个固定于相机下，用于裂隙图片的提取，1个用于其它数据的测量。

4.2 图像提取与形态参数分析

试验所取土样来南京某高速绕城公路路基，属弱膨胀土。本文选取土样初始条件为土样含水量ω=30%，干密度ρd=1.54 g/cm3的裂隙开展试验过程中部分照片来说明裂隙图像的识别分析过程。

在整个试验中，针对裂隙的发育演化过程定时拍照。所得裂隙数码照片在处理过程中可一次统一裁剪、识别、分析。

图9至图11分别为干缩裂隙原始图像、裂隙桥接去杂后的二值图像，以及裂隙骨架提取图。根据裂隙二值图像和裂隙骨架图像，系统对裂隙特征参数(裂隙率长度l、平均宽度d、隙宽分布率等)进行提取，其结果如表1、表2所示。

图9 试样表面裂隙原始图像(ω: 土样含水率)

图10 试样表面裂隙二值化图像(δ: 裂隙率)

图11 t=2 460 min时试样表面裂隙骨架图像

表1裂隙形态参数统计

Table1Statisticsofcrack’smorphologicalparameters

时间t/minω/%δ/%l/cmd/cmD030.000.1030.003024.580.8190.40.351.019018.083.22152.00.471.1218013.355.10269.90.871.203909.558.52371.91.441.382 4608.2812.80436.81.731.42

将系统所得到的各裂隙特征参数与土体含水率、基质吸力、强度及渗透系数等土性参数建立联系，可进一步揭示土体表面干缩裂隙开展规律及其对土体物理性质的影响，为后期裂隙动态分析提供准确依据。

表2 隙宽分布率统计

5 结 论

(1) 借助数字图像处理技术，基于MATLAB软件开发了裂隙网络图像分析系统，通过图像二值化、桥接去杂、裂隙骨架提取等处理程序，实现了土体平面裂隙图像的连续、无损、自动化识别。

(2) 在裂隙二值化图像基础上，对平面裂隙形态参数进行统计，可得到裂隙率、裂隙长度、平均宽度、隙宽分布率、裂隙分维等基本形态参数。

(3) 对不同时间点的裂隙图片进行统一分析，所得到的形态参数为土体干缩裂隙的动态演化分析提供准确依据。

(4) 由于拍摄手段的限制，目前还只能对平面裂隙形态参数进行分析。沿裂隙深度方向，亦即三维几何形态的发展演化分析方法还需进一步探讨。

参考文献：

[1] 袁俊平，殷宗泽. 考虑裂隙非饱和膨胀土边坡入渗模型与数值模拟[J]. 岩土力学, 2004，25(10): 1581-1586.(YUAN Jun-ping, YIN Zong-ze. Numerical Model and Simulation of Expansive Soils Slope Infiltration Considering Fissures[J]. Rock and Soil Mechanics，2004，25(10): 1581-1586.(in Chinese))

[2] 薛 强，梁 冰，刘建军. 裂隙岩体中污染物运移过程的数值模拟[J].岩土力学，2003, 24(增刊):35-38(XUE Qiang, LIANG Bing, LIU Jian-jun. Numerical Simulation of Contaminant Transport in Jointed Rockmass[J]. Rock and Soil Mechanics, 2003, 24(Sup.): 35-38. (in Chinese))

[3] SMALING E M A, BOUMA J. Bypass Flow and Leaching of Nitrogen in a Kenyan Vertisol at the Onset of Growing Season [J]. Soil Use and Management, 1992, 8(1): 44-47．

[4] 姚海林, 郑少河, 陈守义. 考虑裂隙及雨水渗入影响的膨胀土坡稳定性分析[J]. 岩土工程学报, 2001, 23(5): 606-609. (YAO Hai-lin, ZHENG Shao-he, CHEN Shou-yi. Analysis on the Slope Stability of Expansive Soils Considering Cracks and Infiltration of Rain [J]. Chinese Journal of Geotechnical Engineering, 2001, 23(5): 606-609. (in Chinese))

[5] INOUE H. Lateral Water Flow in a Clayey Agricultural Field with Cracks[J].Geoderma,1993,59(1-4):311-325.

[6] VOGEL H, HOFFMANN H, ROTH K. Studies of Crack Dynamics in Clay Soil. I. Experimental Methods, Results, and Morphological Quantification[J]. Geoderma, 2005, 125(3/4): 203-211.

[7] 裴银鸽,徐则民,张家明. 土体干缩裂缝研究进展[J]. 应用生态学, 2012 , 23(4): 1137-1144. (PEI Yin-ge, XU Ze-min, ZHANG Jia-ming. Desiccation Cracking of Soil Body: A Review [J]. Chinese Journal of Applied Ecology, 2012, 23(4): 1137-1144. (in Chinese))

[8] 杨更社. 岩体损伤及检测 [M]. 西安 : 陕西科学出版社 , 1998.(YANG Geng-she. Damage and Detection of Rock Mass [M]. Xi’an: Shaanxi Science Press, 1998. (in Chinese))

[9] 袁俊平，殷宗泽.膨胀土裂隙的量化指标与强度研究[J]，水利学报，2004，(6)：108-113.(YUAN Jun-ping, YIN Zong-ze. Quantitative Index of Fissure and Strength Characteristics of Fissured Expansive Soils [J]. Journal of Hydraulic Engineering, 2004，(6)：108-113. (in Chinese))

[10] 唐朝生，崔玉军，TANG A M，等. 膨胀土收缩开裂过程及其温度效应[J]，岩土工程学报，2012, 34(12):2181-2187.(TANG Chao-sheng, CUI Yu-jun, TANG A M,etal. Shrinkage and Desiccation Cracking Process of Expansive Soil and Its Temperature-dependent Behavior [J]. Chinese Journal of Geotechnical Engineering, 2012, 34(12): 2181-2187.(in Chinese))

[11] 王金安, 谢和平, 田晓燕, 等. 岩石断裂表面分形测量的尺度效应[J]. 岩石力学与工程学报, 2000, 19(1): 11-17. (WANG Jin-an, XIE He-ping, TIAN Xiao-yan,etal. Scale Effect on Fractal Measurement of Rock Fracture Surfaces [J]. Journal of Rock Mechanics and Engineering, 2000, 19(1): 11-17.(in Chinese))

[12] 刘 春，王宝军，施 斌，等.基于数字图像识别的岩土体裂隙形态参数分析方法[J]，岩土工程学报，2008,30(9)：1384-1388.(LIU Chun, WANG Bao-jun, SHI Bin,etal. Analytic Method of Morphological Parameters of Cracks for Rock and Soil Based on Image Processing and Recognition [J]. Chinese Journal of Geotechnical Engineering, 2008, 30(9): 1384-1388.(in Chinese))