基于形态学理论的土体裂隙边缘分形维数计算

刘 兴，王 媛，冯 迪

(1.河海大学岩土力学与堤坝工程教育部重点实验室，江苏南京 210098;2.河海大学土木与交通学院，江苏南京 210098)

土体中广泛发育着裂隙，裂隙是影响土体强度和渗透性的主要因素。定量描述裂隙形态是进一步揭示土体裂隙开展机理及其演化规律的前提。长期以来，由于裂隙扩展的随机性和非均匀性，导致其定量化程度较低。

土体裂隙具有很强的非线性特征。已有研究表明［1-4］土体表面的裂隙符合分形理论所要求的“自相似性”，因此很多学者尝试应用分形理论对土体的力学效应及工程性质展开相关研究［5-8］。有效利用分形理论进行裂隙参数定量化统计的前提是能够准确提取裂隙的发育形态。研究初期，对土体裂隙形态的获取多是以手工描绘为主［9］。随着计算机技术的发展，例如AutoCAD这样的绘图辅助软件逐渐被应用到裂隙形态的描绘中［10-11］，但无论是手工绘制还是软件绘图均易受人为因素的影响，引起的随机误差较大，并且精度较低。数码图像信息处理技术的发展有望克服这些不足。已有学者应用这种技术在相关领域进行研究，譬如:彭瑞东等［12］借助数字图像技术从SEM图中提取岩石开裂面裂纹的扩展曲线;刘春等［13］提出一套裂隙图像计算机识别和定量分析方法，对裂隙图像进行形态学操作，获取裂隙长度、宽度、走向等裂隙形态参数。

以上研究多侧重于裂隙单一基本参数的提取，仅凭某几个孤立的指标还不足以反映复杂裂隙的几何形态，难以全面体现裂隙发育的随机行为和细节特征，但目前对综合反映裂隙整体形态的指标鲜有研究。鉴于此，笔者基于数学形态学理论，对裂隙图像进行处理，获取裂隙边缘曲线，把裂隙边缘的分形维数作为裂隙发育程度的表征指标，建立了边缘分形维数与裂隙宽度之间的关系，并尝试将裂隙边缘的分形维数应用到立方定律的修正中。

1 裂隙边缘分析方法

数学形态学不仅是一种理论，而且是一门强大的图像分析技术。边缘检测是进行图像信息提取的基本手段，它可以剔除非相关信息，保留图像的重要属性。边缘检测通常借用微分算子利用卷积来实现，笔者采用的Sobel算子对噪声具有平滑作用，可以提供较精确的边缘方向信息，且计算简单、速度快。

1.1 Sobel算子

Sobel算子是一阶导数的边缘检测算子，该算子包含一对3×3矩阵，分别用来检测水平边缘和垂直边缘，利用式(1)计算某一像素邻近区域的梯度值来确定该像素的梯度，通过梯度值的大小与设定的阈值(T)进行比较、取舍，结果大于设定阈值的便是边缘部分。

式中:G(x，y)——(x，y)像素点的梯度值;S'x，S'y——检测水平方向与垂直方向的有向算子;T——程序调试过程中确定的合理阈值。

1.2 裂隙边缘提取

边缘检测是要计算图像每一像素的梯度值，数字图像中的每个像素都对应一个灰度值，若将每一像素的灰度值预先处理成简单数值，便会给梯度的计算带来方便。二值图像中所有元素的灰度值只能从0和1中取值，其灰度值没有中间过渡，图像二值化在图像处理中相当于将被关注的区域与背景相分离，所以在图像边缘提取之前将图像二值化往往会达到更好的效果。笔者在程序编制过程中亦遵循这种顺序。

裂隙图像处理结果如图1所示。图1(a)为待处理图像，先处理成图1(b)所示的二值图像(黑色区域表示土体裂隙，纯白色表示非裂隙区域)。图1(c)中白色曲线即裂隙边缘，与图1(a)比较，基于Sobel算子提取的裂隙边缘能准确地反映原始裂隙的宽度、走向，以及边缘曲折等几何形态。

图1 裂隙图像处理示意图Fig.1 Sketch map of cracking image processing

2 裂隙边缘分形维数计算

2.1 盒维数

很多数学家虽然在分形维数的定义方法上做过卓有成效的工作，但由于存在于自然界中的物理分形往往表现出某种随机性和尺度性，分形的确切定义仍未给出，有的定义在理论上比较科学，但实际计算中却很困难，如Hausdorff维数、相似维数，因此很多等价的或近似的分形维数定义被提出来，盒维数是目前广泛采用的一种分形维数计算方法。

设F是Rn上任意非空的有界子集，Nδ是边长最大为δ，可以覆盖F集的最少盒子数，逐渐减小边长δ，相对应地，统计能够覆盖F集的盒子数Nδ，即可用式(2)计算盒维数:

根据盒维数的定义可知，通过逐渐细分盒子边长可对应得到一系列覆盖原图像所需的盒子数，以获取足够多的数据点(δi，Nδi)。但在实际计算中受图像的分辨率限制，只能划分到一个像素点尺寸。原则上，图像分辨率越高，计算精度也越高，δ的取值与图像的大小以及形状有很大关系，当图像较大并且为正方形时，为避免出现半个盒子的情况，δ可以取能被正方形边长所整除的数字，而当图像为矩形时δ可取矩形两边长度的公约数。

采集数据点(δi，Nδi)，并在双对数坐标平面内用最小二乘法线性拟合，所得直线斜率的绝对值即该图像的盒维数D。

2.2 分形维数计算软件设计

通过以上对盒维数的概念分析，基于VisualC++语言编制了盒维数的计算软件，该软件具有可实现图像读取、图像二值化、输出“盒子”边长、计算盒维数、得出最小二乘法拟合曲线等一系列功能。软件界面见图2。

软件可自动计算图像两边长度的公约数，并将其作为“盒子”的边长，有效避免了格子数无法除尽导致的统计误差，并且不再局限于处理正方形图像。

图2 软件界面Fig.2 Software interface

3 实例分析

为再现土体在自然条件下裂隙的发育过程，进一步研究裂隙边缘分形维数与裂隙宽度之间的关系，进行了土体开裂观测试验。土样取自镇江南徐大道黄山边坡原状土样，取土深度为地表下3 m，从取样表观分析可知，土样裂隙发育、颜色多为棕黄色，经刀片切削后呈现的土体表面呈蜡状光滑，具有典型的裂隙土特征。土样的基本物理力学指标如下:液限wL=41.3%，塑限wP=24.8%，自由膨胀率r=50% ～90%，密度ρ=1930 kg/m3，黏聚力 c=37.3 kPa，内摩擦角 θ=14.4°。

试验在四五月份完成，水分蒸发过程较均匀，试验开始后第2天出现裂隙，以后每24 h拍照记录1次，待裂隙宽度波动量持续3 d小于2%时认为裂隙不再开展。为避免边界效应对裂隙开展的影响，选取中间区域为观测研究区域，用上述方法提取每张图片的裂隙边缘曲线，部分结果如图3所示。

拍摄图片的同时，取观测区域外的土体进行含水率测定，并统计裂隙平均宽度、分形维数等参数，试验结果如图4、图5所示。

试验结果分析如下:如图4所示，含水率随着试验进行呈逐渐减少趋势，试验进行到第8天趋于稳定，第12天认为干缩完成;裂隙宽度在土体开裂初期增长较快，随着干缩时间的推移其增势减缓，最终稳定在65像素左右;初始裂隙边缘呈齿状形态，盒维数值开始较大，随后逐渐减小，说明裂隙边缘逐渐趋于平缓。对比图4(a)和图4(b)可知，含水率的微小波动会引起裂隙宽度的增减，说明将形态学理论应用于裂隙边缘提取及对裂隙宽度的统计是足够精确的。

图3 提取得到的不同干缩状态下的裂隙边缘Fig.3 Extraction of crack edge in different shrinking states

图4 含水率、裂隙宽度像素、盒维数随时间变化曲线Fig.4 Changes of water content，crack width pixel，and box dimension with time

在现有的研究成果中，普遍认为基质吸力是引起土体开裂的主要因素，由非饱和土的土－水特征曲线可知，基质吸力与含水率密切相关。本文试验的统计结果也验证了这一理论。如图5(a)所示，含水率随干缩时间的增长而逐渐降低，体现在非饱和土力学中即反映了基质吸力的增加，进而引起土体开裂，裂隙宽度不断增大。裂隙边缘盒维数值随含水率的降低呈现减小趋势，如图5(b)所示。需要指出的是，在试验过程中并未控制湿度的变化，期间含水率会有小幅度变动。

图5 含水率、裂隙宽度像素、盒维数之间的关系Fig.5 Relationships between water content，crack width pixel，and box dimension

含水率、裂隙宽度像素、盒维数三者之间均有较好的线性关系。裂隙边缘D随裂隙宽度像素b的增加而减小，如图5(c)所示，两者有如下负相关线性关系:

D和b两者相关系数为0.9663，裂隙边缘盒维数能够作为表征裂隙宽度的度量指标。一般来说，表面越粗糙，盒维数越大，式(3)揭示土体裂隙宽度越大，盒维数越小，也就是说随着裂隙开展，裂隙表面越平滑，这与文献［14］所述的有关压裂岩体裂隙分形规律有所出入，究其原因，可能与岩土体开裂机理不同有关。土体开展多属张拉破坏，且土体抗拉强度较小，一旦开裂易被拉伸，导致裂隙表面较平展;前述岩体为压裂破坏，随着压力不断增大，裂隙咬合面越粗糙。

4 立方定律的修正

立方定律是裂隙渗流研究的基础理论，根据光滑平行板裂隙水力学试验推导而来。但天然裂隙表面往往是凹凸起伏的，粗糙程度极大地影响了裂隙的水力特性。若将立方定律应用到实际中必须对它进行修正，修正需建立在裂隙面粗糙性的定量描述上。Barton等［15］建议的节理粗糙度系数JRC已经被国际岩石力学学会所推荐使用。

许多学者［16－20］已经注意到可以采用分形维数的大小来表征裂隙面的粗糙程度并试图将分形数维与JRC之间建立某种关系，以达到定量化研究裂隙面形态的目的，两者之间的关系在理论研究上已经较成熟，但如何准确获取现场实际裂隙的分形维数是目前面临的主要问题。以形态学理论提取的裂隙边缘曲线能较好地反映裂隙形态，其分形维数亦能准确表征裂隙的粗糙程度，并与裂隙宽度有较好的线性关系，这就为立方定律的修正提供了可靠参数。

必须指出，现有的研究成果均是描述裂隙表面的曲线形态，由于裂隙张开度很小，曲面形态的研究尚存在困难。如何将分形理论运用到曲面粗糙度的描述上将是下一步要进行的工作。

5 结 论

a.基于数学形态学理论，利用Sobel算子对裂隙边缘曲线进行提取，经过与原始图像对比，发现所提取的边缘符合裂隙开展形态，并能较好地反映裂隙边缘的粗糙程度。

b.编制了计算数字图像盒维数的软件，该软件具有良好的人机交互界面，只需读入256色原始BMP图像即可实现盒维数的计算过程。

c.试验表明，裂隙边缘分形维数与裂隙宽度有较好的负线性关系，因此分形维数可以被用来刻画裂隙的发育程度，进一步研究将有助于边缘分形维数与粗糙度系数JRC建立关系，从而达到对立方定律修正的目的。

［1］袁俊平，殷宗泽，包承纲.膨胀土裂隙的量化手段与度量指标研究［J］.长江科学院院报，2003，20(6):27-30.(YUAN Junping，YIN Zongze，BAO Chenggang.Quantitative description method and index for fissures in expansive soils［J］.Journal of Yangtse River Scientific Research Institute，2003，20(6):27-30.(in Chinese))

［2］袁俊平.非饱和膨胀土的裂隙概化模型与边坡稳定研究［D］.南京:河海大学，2003.

［3］BAER J U，KENT T F，ANDERSON S H.Image analysis and fractal geometry to characterize soil desiccation cracks［J］.Geoderma，2009，154:153-163.

［4］PERFECT E.Fractal models for the fragmentation of rocks and soils:a review［J］.Engineering Geology，1997，48:185-198.

［5］徐永福，孙婉莹，吴正根.我国膨胀土的分形结构的研究［J］.河海大学学报:自然科学版，1997，25(1):18-23.(XU Yongfu，SUN Wanying，WU Zhenggen.On fractal structure of expansive soils in China［J］.Journal of Hohai University:Natural Sciences，1997，25(1):18-23.(in Chinese))

［6］胡卸文，王治平.黏土体裂隙效应中的分形几何现象［J］.四川水力发电，1998，17(1):22-26.(HUXiewen，WANGZhiping.Fracture geometry properties in fissure effect of clay mass［J］.Sichuan Water Power，1998，17(1):22-26.(in Chinese))

［7］易顺民，黎志恒，张延中.膨胀土裂隙结构的分形特征及其意义［J］.岩土工程学报，1999，21(3):294-298.(YI Shunmin，LI Zhiheng，ZHANG Yanzhong.The fractal characteristics of fractures in expansion soil and its significance［J］.Chinese Journal of Geotechnical Engineering，1999，21(3):294-298.(in Chinese))

［8］PRESTON S，GRIFFITHSB S，YOUNG I M.An investigation into sources of soil crack heterogeneity using fractal geometry［J］.European Journal of Soil Science，1997，48:31-37.

［9］PRIESTSD，HUDSON JA.Estimation of discontinuity spacing and trace length using scanline surveys［J］.International Journal of Rock Mechanics and Mining Sciences ＆ Geomechanics Abstracts，1981，18(3):183-197.

［10］冯迪.裂隙土的裂隙几何描述和渗透特性研究［D］.南京:河海大学，2009.

［11］LI Jinhui.Field experimental study and numerical simulation of seepage in saturated/unsaturated cracked soil［D］.Hong Kong:The Hong Kong University of Science and Technology，2009.

［12］彭瑞东，谢和平，鞠杨.二维数字图像分形维数的计算方法［J］.中国矿业大学学报，2004，33(1):19-24.(PENG Ruidong，XIE Heping，JU Yang.Computation method of fractal dimension for 2-D digital image［J］.Journal of China University of Mining and Technology，2004，33(1):19-24.(in Chinese))

［13］刘春，王宝军，施斌，等.基于数字图像识别的岩土体裂隙形态参数分析方法［J］.岩土工程学报，2008，30(9):1383-1388.(LIU Chun，WANG Baojun，SHI Bin，et al.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):1383-1388.(in Chinese))

［14］李玮，闫铁，毕雪亮.基于分形方法的水力压裂裂缝起裂扩展机理［J］.中国石油大学学报:自然科学版，2008，32(5):87-91.(LI Wei，YAN Tie，BI Xueliang.Mechanism of hydraulically created fracture breakdown and propagation based on fractal method［J］.Journal of China University of Petroleum:Edition of Natural Science，2008，32(5):87-91.(in Chinese))

［15］BARTON N，CHOUBEY V.The shear strength of rock joints in theory and practice［J］.Rock Mechanics，1977，10(1/2):1-54.

［16］王媛，速宝玉.单裂隙面渗流特性及等效水力隙宽［J］.水科学进展，2002，13(1):61-68.(WANG Yuan，SU Baoyu.Research on the behavior of fluid flow in a single fracture and its equivalent hydraulic aperture［J］.Advances in Water Science，2002，13(1):61-68.(in Chinese))

［17］CARR J R.Fractal characterization of joint surface roughness in welded tuff at Yucca Mountain Nevada［EB/OL］.［2013-04-02］.http://www.onepetro.org/mslib/servlet/onepetropreview?id=ARMA-89-0193.

［18］LEE Y H.The fractal dimension as a measure of the roughness of rock discontinuity profiles［J］.International Journal of Rock Mechanics and Mining Sciences ＆ Geomechanics Abstracts，1990，27(6):453-464.

［19］谢和平，WILLIAN G.节理粗糙度系数的分形估算［J］.地质科学译丛，1992，9(1):85-90.(XIE Heping，WILLIAN G.Fractal estimate of the joint roughness coefficient［J］.Journal of Geoscience Translations，1992，9(1):85-90.(in Chinese))

［20］周创兵，熊文林.不连续面的分形维数及其在渗流分析中的应用［J］.水文地质工程地质，1996(6):1-5.(ZHOU Chuangbing，XIONG Wenlin.Discontinuous surface fractal dimension and its application in the seepage analysis ［J］.Hydrogeology and Engineering Geology，1996(6):1-5.(in Chinese))