基于形態(tài)學(xué)理論的土體裂隙邊緣分形維數(shù)計(jì)算

劉 興，王 媛，馮 迪

(1.河海大學(xué)巖土力學(xué)與堤壩工程教育部重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210098;2.河海大學(xué)土木與交通學(xué)院，江蘇南京 210098)

土體中廣泛發(fā)育著裂隙，裂隙是影響土體強(qiáng)度和滲透性的主要因素。定量描述裂隙形態(tài)是進(jìn)一步揭示土體裂隙開(kāi)展機(jī)理及其演化規(guī)律的前提。長(zhǎng)期以來(lái)，由于裂隙擴(kuò)展的隨機(jī)性和非均勻性，導(dǎo)致其定量化程度較低。

土體裂隙具有很強(qiáng)的非線(xiàn)性特征。已有研究表明［1-4］土體表面的裂隙符合分形理論所要求的“自相似性”，因此很多學(xué)者嘗試應(yīng)用分形理論對(duì)土體的力學(xué)效應(yīng)及工程性質(zhì)展開(kāi)相關(guān)研究［5-8］。有效利用分形理論進(jìn)行裂隙參數(shù)定量化統(tǒng)計(jì)的前提是能夠準(zhǔn)確提取裂隙的發(fā)育形態(tài)。研究初期，對(duì)土體裂隙形態(tài)的獲取多是以手工描繪為主［9］。隨著計(jì)算機(jī)技術(shù)的發(fā)展，例如AutoCAD這樣的繪圖輔助軟件逐漸被應(yīng)用到裂隙形態(tài)的描繪中［10-11］，但無(wú)論是手工繪制還是軟件繪圖均易受人為因素的影響，引起的隨機(jī)誤差較大，并且精度較低。數(shù)碼圖像信息處理技術(shù)的發(fā)展有望克服這些不足。已有學(xué)者應(yīng)用這種技術(shù)在相關(guān)領(lǐng)域進(jìn)行研究，譬如:彭瑞東等［12］借助數(shù)字圖像技術(shù)從SEM圖中提取巖石開(kāi)裂面裂紋的擴(kuò)展曲線(xiàn);劉春等［13］提出一套裂隙圖像計(jì)算機(jī)識(shí)別和定量分析方法，對(duì)裂隙圖像進(jìn)行形態(tài)學(xué)操作，獲取裂隙長(zhǎng)度、寬度、走向等裂隙形態(tài)參數(shù)。

以上研究多側(cè)重于裂隙單一基本參數(shù)的提取，僅憑某幾個(gè)孤立的指標(biāo)還不足以反映復(fù)雜裂隙的幾何形態(tài)，難以全面體現(xiàn)裂隙發(fā)育的隨機(jī)行為和細(xì)節(jié)特征，但目前對(duì)綜合反映裂隙整體形態(tài)的指標(biāo)鮮有研究。鑒于此，筆者基于數(shù)學(xué)形態(tài)學(xué)理論，對(duì)裂隙圖像進(jìn)行處理，獲取裂隙邊緣曲線(xiàn)，把裂隙邊緣的分形維數(shù)作為裂隙發(fā)育程度的表征指標(biāo)，建立了邊緣分形維數(shù)與裂隙寬度之間的關(guān)系，并嘗試將裂隙邊緣的分形維數(shù)應(yīng)用到立方定律的修正中。

1 裂隙邊緣分析方法

數(shù)學(xué)形態(tài)學(xué)不僅是一種理論，而且是一門(mén)強(qiáng)大的圖像分析技術(shù)。邊緣檢測(cè)是進(jìn)行圖像信息提取的基本手段，它可以剔除非相關(guān)信息，保留圖像的重要屬性。邊緣檢測(cè)通常借用微分算子利用卷積來(lái)實(shí)現(xiàn)，筆者采用的Sobel算子對(duì)噪聲具有平滑作用，可以提供較精確的邊緣方向信息，且計(jì)算簡(jiǎn)單、速度快。

1.1 Sobel算子

Sobel算子是一階導(dǎo)數(shù)的邊緣檢測(cè)算子，該算子包含一對(duì)3×3矩陣，分別用來(lái)檢測(cè)水平邊緣和垂直邊緣，利用式(1)計(jì)算某一像素鄰近區(qū)域的梯度值來(lái)確定該像素的梯度，通過(guò)梯度值的大小與設(shè)定的閾值(T)進(jìn)行比較、取舍，結(jié)果大于設(shè)定閾值的便是邊緣部分。

式中:G(x，y)——(x，y)像素點(diǎn)的梯度值;S'x，S'y——檢測(cè)水平方向與垂直方向的有向算子;T——程序調(diào)試過(guò)程中確定的合理閾值。

1.2 裂隙邊緣提取

邊緣檢測(cè)是要計(jì)算圖像每一像素的梯度值，數(shù)字圖像中的每個(gè)像素都對(duì)應(yīng)一個(gè)灰度值，若將每一像素的灰度值預(yù)先處理成簡(jiǎn)單數(shù)值，便會(huì)給梯度的計(jì)算帶來(lái)方便。二值圖像中所有元素的灰度值只能從0和1中取值，其灰度值沒(méi)有中間過(guò)渡，圖像二值化在圖像處理中相當(dāng)于將被關(guān)注的區(qū)域與背景相分離，所以在圖像邊緣提取之前將圖像二值化往往會(huì)達(dá)到更好的效果。筆者在程序編制過(guò)程中亦遵循這種順序。

裂隙圖像處理結(jié)果如圖1所示。圖1(a)為待處理圖像，先處理成圖1(b)所示的二值圖像(黑色區(qū)域表示土體裂隙，純白色表示非裂隙區(qū)域)。圖1(c)中白色曲線(xiàn)即裂隙邊緣，與圖1(a)比較，基于Sobel算子提取的裂隙邊緣能準(zhǔn)確地反映原始裂隙的寬度、走向，以及邊緣曲折等幾何形態(tài)。

圖1 裂隙圖像處理示意圖Fig.1 Sketch map of cracking image processing

2 裂隙邊緣分形維數(shù)計(jì)算

2.1 盒維數(shù)

很多數(shù)學(xué)家雖然在分形維數(shù)的定義方法上做過(guò)卓有成效的工作，但由于存在于自然界中的物理分形往往表現(xiàn)出某種隨機(jī)性和尺度性，分形的確切定義仍未給出，有的定義在理論上比較科學(xué)，但實(shí)際計(jì)算中卻很困難，如Hausdorff維數(shù)、相似維數(shù)，因此很多等價(jià)的或近似的分形維數(shù)定義被提出來(lái)，盒維數(shù)是目前廣泛采用的一種分形維數(shù)計(jì)算方法。

設(shè)F是Rn上任意非空的有界子集，Nδ是邊長(zhǎng)最大為δ，可以覆蓋F集的最少盒子數(shù)，逐漸減小邊長(zhǎng)δ，相對(duì)應(yīng)地，統(tǒng)計(jì)能夠覆蓋F集的盒子數(shù)Nδ，即可用式(2)計(jì)算盒維數(shù):

根據(jù)盒維數(shù)的定義可知，通過(guò)逐漸細(xì)分盒子邊長(zhǎng)可對(duì)應(yīng)得到一系列覆蓋原圖像所需的盒子數(shù)，以獲取足夠多的數(shù)據(jù)點(diǎn)(δi，Nδi)。但在實(shí)際計(jì)算中受圖像的分辨率限制，只能劃分到一個(gè)像素點(diǎn)尺寸。原則上，圖像分辨率越高，計(jì)算精度也越高，δ的取值與圖像的大小以及形狀有很大關(guān)系，當(dāng)圖像較大并且為正方形時(shí)，為避免出現(xiàn)半個(gè)盒子的情況，δ可以取能被正方形邊長(zhǎng)所整除的數(shù)字，而當(dāng)圖像為矩形時(shí)δ可取矩形兩邊長(zhǎng)度的公約數(shù)。

采集數(shù)據(jù)點(diǎn)(δi，Nδi)，并在雙對(duì)數(shù)坐標(biāo)平面內(nèi)用最小二乘法線(xiàn)性擬合，所得直線(xiàn)斜率的絕對(duì)值即該圖像的盒維數(shù)D。

2.2 分形維數(shù)計(jì)算軟件設(shè)計(jì)

通過(guò)以上對(duì)盒維數(shù)的概念分析，基于VisualC++語(yǔ)言編制了盒維數(shù)的計(jì)算軟件，該軟件具有可實(shí)現(xiàn)圖像讀取、圖像二值化、輸出“盒子”邊長(zhǎng)、計(jì)算盒維數(shù)、得出最小二乘法擬合曲線(xiàn)等一系列功能。軟件界面見(jiàn)圖2。

軟件可自動(dòng)計(jì)算圖像兩邊長(zhǎng)度的公約數(shù)，并將其作為“盒子”的邊長(zhǎng)，有效避免了格子數(shù)無(wú)法除盡導(dǎo)致的統(tǒng)計(jì)誤差，并且不再局限于處理正方形圖像。

圖2 軟件界面Fig.2 Software interface

3 實(shí)例分析

為再現(xiàn)土體在自然條件下裂隙的發(fā)育過(guò)程，進(jìn)一步研究裂隙邊緣分形維數(shù)與裂隙寬度之間的關(guān)系，進(jìn)行了土體開(kāi)裂觀(guān)測(cè)試驗(yàn)。土樣取自鎮(zhèn)江南徐大道黃山邊坡原狀土樣，取土深度為地表下3 m，從取樣表觀(guān)分析可知，土樣裂隙發(fā)育、顏色多為棕黃色，經(jīng)刀片切削后呈現(xiàn)的土體表面呈蠟狀光滑，具有典型的裂隙土特征。土樣的基本物理力學(xué)指標(biāo)如下:液限wL=41.3%，塑限wP=24.8%，自由膨脹率r=50% ～90%，密度ρ=1930 kg/m3，黏聚力 c=37.3 kPa，內(nèi)摩擦角 θ=14.4°。

試驗(yàn)在四五月份完成，水分蒸發(fā)過(guò)程較均勻，試驗(yàn)開(kāi)始后第2天出現(xiàn)裂隙，以后每24 h拍照記錄1次，待裂隙寬度波動(dòng)量持續(xù)3 d小于2%時(shí)認(rèn)為裂隙不再開(kāi)展。為避免邊界效應(yīng)對(duì)裂隙開(kāi)展的影響，選取中間區(qū)域?yàn)橛^(guān)測(cè)研究區(qū)域，用上述方法提取每張圖片的裂隙邊緣曲線(xiàn)，部分結(jié)果如圖3所示。

拍攝圖片的同時(shí)，取觀(guān)測(cè)區(qū)域外的土體進(jìn)行含水率測(cè)定，并統(tǒng)計(jì)裂隙平均寬度、分形維數(shù)等參數(shù)，試驗(yàn)結(jié)果如圖4、圖5所示。

試驗(yàn)結(jié)果分析如下:如圖4所示，含水率隨著試驗(yàn)進(jìn)行呈逐漸減少趨勢(shì)，試驗(yàn)進(jìn)行到第8天趨于穩(wěn)定，第12天認(rèn)為干縮完成;裂隙寬度在土體開(kāi)裂初期增長(zhǎng)較快，隨著干縮時(shí)間的推移其增勢(shì)減緩，最終穩(wěn)定在65像素左右;初始裂隙邊緣呈齒狀形態(tài)，盒維數(shù)值開(kāi)始較大，隨后逐漸減小，說(shuō)明裂隙邊緣逐漸趨于平緩。對(duì)比圖4(a)和圖4(b)可知，含水率的微小波動(dòng)會(huì)引起裂隙寬度的增減，說(shuō)明將形態(tài)學(xué)理論應(yīng)用于裂隙邊緣提取及對(duì)裂隙寬度的統(tǒng)計(jì)是足夠精確的。

圖3 提取得到的不同干縮狀態(tài)下的裂隙邊緣Fig.3 Extraction of crack edge in different shrinking states

圖4 含水率、裂隙寬度像素、盒維數(shù)隨時(shí)間變化曲線(xiàn)Fig.4 Changes of water content，crack width pixel，and box dimension with time

在現(xiàn)有的研究成果中，普遍認(rèn)為基質(zhì)吸力是引起土體開(kāi)裂的主要因素，由非飽和土的土－水特征曲線(xiàn)可知，基質(zhì)吸力與含水率密切相關(guān)。本文試驗(yàn)的統(tǒng)計(jì)結(jié)果也驗(yàn)證了這一理論。如圖5(a)所示，含水率隨干縮時(shí)間的增長(zhǎng)而逐漸降低，體現(xiàn)在非飽和土力學(xué)中即反映了基質(zhì)吸力的增加，進(jìn)而引起土體開(kāi)裂，裂隙寬度不斷增大。裂隙邊緣盒維數(shù)值隨含水率的降低呈現(xiàn)減小趨勢(shì)，如圖5(b)所示。需要指出的是，在試驗(yàn)過(guò)程中并未控制濕度的變化，期間含水率會(huì)有小幅度變動(dòng)。

圖5 含水率、裂隙寬度像素、盒維數(shù)之間的關(guān)系Fig.5 Relationships between water content，crack width pixel，and box dimension

含水率、裂隙寬度像素、盒維數(shù)三者之間均有較好的線(xiàn)性關(guān)系。裂隙邊緣D隨裂隙寬度像素b的增加而減小，如圖5(c)所示，兩者有如下負(fù)相關(guān)線(xiàn)性關(guān)系:

D和b兩者相關(guān)系數(shù)為0.9663，裂隙邊緣盒維數(shù)能夠作為表征裂隙寬度的度量指標(biāo)。一般來(lái)說(shuō)，表面越粗糙，盒維數(shù)越大，式(3)揭示土體裂隙寬度越大，盒維數(shù)越小，也就是說(shuō)隨著裂隙開(kāi)展，裂隙表面越平滑，這與文獻(xiàn)［14］所述的有關(guān)壓裂巖體裂隙分形規(guī)律有所出入，究其原因，可能與巖土體開(kāi)裂機(jī)理不同有關(guān)。土體開(kāi)展多屬?gòu)埨茐?，且土體抗拉強(qiáng)度較小，一旦開(kāi)裂易被拉伸，導(dǎo)致裂隙表面較平展;前述巖體為壓裂破壞，隨著壓力不斷增大，裂隙咬合面越粗糙。

4 立方定律的修正

立方定律是裂隙滲流研究的基礎(chǔ)理論，根據(jù)光滑平行板裂隙水力學(xué)試驗(yàn)推導(dǎo)而來(lái)。但天然裂隙表面往往是凹凸起伏的，粗糙程度極大地影響了裂隙的水力特性。若將立方定律應(yīng)用到實(shí)際中必須對(duì)它進(jìn)行修正，修正需建立在裂隙面粗糙性的定量描述上。Barton等［15］建議的節(jié)理粗糙度系數(shù)JRC已經(jīng)被國(guó)際巖石力學(xué)學(xué)會(huì)所推薦使用。

許多學(xué)者［16－20］已經(jīng)注意到可以采用分形維數(shù)的大小來(lái)表征裂隙面的粗糙程度并試圖將分形數(shù)維與JRC之間建立某種關(guān)系，以達(dá)到定量化研究裂隙面形態(tài)的目的，兩者之間的關(guān)系在理論研究上已經(jīng)較成熟，但如何準(zhǔn)確獲取現(xiàn)場(chǎng)實(shí)際裂隙的分形維數(shù)是目前面臨的主要問(wèn)題。以形態(tài)學(xué)理論提取的裂隙邊緣曲線(xiàn)能較好地反映裂隙形態(tài)，其分形維數(shù)亦能準(zhǔn)確表征裂隙的粗糙程度，并與裂隙寬度有較好的線(xiàn)性關(guān)系，這就為立方定律的修正提供了可靠參數(shù)。

必須指出，現(xiàn)有的研究成果均是描述裂隙表面的曲線(xiàn)形態(tài)，由于裂隙張開(kāi)度很小，曲面形態(tài)的研究尚存在困難。如何將分形理論運(yùn)用到曲面粗糙度的描述上將是下一步要進(jìn)行的工作。

5 結(jié) 論

a.基于數(shù)學(xué)形態(tài)學(xué)理論，利用Sobel算子對(duì)裂隙邊緣曲線(xiàn)進(jìn)行提取，經(jīng)過(guò)與原始圖像對(duì)比，發(fā)現(xiàn)所提取的邊緣符合裂隙開(kāi)展形態(tài)，并能較好地反映裂隙邊緣的粗糙程度。

b.編制了計(jì)算數(shù)字圖像盒維數(shù)的軟件，該軟件具有良好的人機(jī)交互界面，只需讀入256色原始BMP圖像即可實(shí)現(xiàn)盒維數(shù)的計(jì)算過(guò)程。

c.試驗(yàn)表明，裂隙邊緣分形維數(shù)與裂隙寬度有較好的負(fù)線(xiàn)性關(guān)系，因此分形維數(shù)可以被用來(lái)刻畫(huà)裂隙的發(fā)育程度，進(jìn)一步研究將有助于邊緣分形維數(shù)與粗糙度系數(shù)JRC建立關(guān)系，從而達(dá)到對(duì)立方定律修正的目的。

［1］袁俊平，殷宗澤，包承綱.膨脹土裂隙的量化手段與度量指標(biāo)研究［J］.長(zhǎng)江科學(xué)院院報(bào)，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］袁俊平.非飽和膨脹土的裂隙概化模型與邊坡穩(wěn)定研究［D］.南京:河海大學(xué)，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］徐永福，孫婉瑩，吳正根.我國(guó)膨脹土的分形結(jié)構(gòu)的研究［J］.河海大學(xué)學(xué)報(bào):自然科學(xué)版，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］胡卸文，王治平.黏土體裂隙效應(yīng)中的分形幾何現(xiàn)象［J］.四川水力發(fā)電，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］易順民，黎志恒，張延中.膨脹土裂隙結(jié)構(gòu)的分形特征及其意義［J］.巖土工程學(xué)報(bào)，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］.南京:河海大學(xué)，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］彭瑞東，謝和平，鞠楊.二維數(shù)字圖像分形維數(shù)的計(jì)算方法［J］.中國(guó)礦業(yè)大學(xué)學(xué)報(bào)，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］劉春，王寶軍，施斌，等.基于數(shù)字圖像識(shí)別的巖土體裂隙形態(tài)參數(shù)分析方法［J］.巖土工程學(xué)報(bào)，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］李瑋，閆鐵，畢雪亮.基于分形方法的水力壓裂裂縫起裂擴(kuò)展機(jī)理［J］.中國(guó)石油大學(xué)學(xué)報(bào):自然科學(xué)版，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］.水科學(xué)進(jìn)展，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.節(jié)理粗糙度系數(shù)的分形估算［J］.地質(zhì)科學(xué)譯叢，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］周創(chuàng)兵，熊文林.不連續(xù)面的分形維數(shù)及其在滲流分析中的應(yīng)用［J］.水文地質(zhì)工程地質(zhì)，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))