基于區(qū)域生長算法的土石混合介質(zhì)細(xì)觀識別與數(shù)值模擬研究

戴 薇 石 崇 阮懷寧 孔 洋 楊俊雄

(1. 河海大學(xué) 巖土力學(xué)與堤壩工程教育部重點實驗室， 南京 210098； 2. 河海大學(xué) 巖土工程科學(xué)研究所，南京 210098)

土石混合介質(zhì)[1-2]是指第四紀(jì)以來形成的，由具有一定工程尺度、強(qiáng)度較高的塊石、細(xì)粒土體及孔隙構(gòu)成且具有一定含石量的極端不均勻松散巖土介質(zhì)系統(tǒng)，是一種由作為骨料的石塊與作為充填料的土組成的地質(zhì)體．由于兩種組分在強(qiáng)度、變形特征等方面的差異，使得土石混合介質(zhì)表現(xiàn)出典型的非均質(zhì)性、非線性性以及非連續(xù)性，同時，在自然界中，不同場地下土石混合介質(zhì)分布情況各不相同，甚至有很大差別，故基于圖像數(shù)字識別技術(shù)建立土石混合介質(zhì)的數(shù)值模型，對于分析不同土石混合介質(zhì)的特性，具有更加普遍的意義．

目前對土石混合介質(zhì)的研究主要分為現(xiàn)場試驗和數(shù)值模擬兩方面．徐文杰等[3]對虎跳峽龍蟠右岸土石混合介質(zhì)分別進(jìn)行了天然狀態(tài)下和浸水條件下的原位大型推剪試驗，獲得了土石混合體強(qiáng)度參數(shù)的參考值并證明土石混合介質(zhì)的強(qiáng)度與含石量有關(guān)．胡峰[4]等考慮不同含石量、不同上覆壓力、不同塊石尺寸進(jìn)行剪切變形試驗，揭示了土石混合體滑坡剪切帶形成演化規(guī)律、破壞模式．但由于現(xiàn)場試驗受到場地以及試驗條件的約束，數(shù)值試驗成為更令人青睞的研究方式，廖秋林等[5]基于基于塊石與土體的顏色屬性差異識別土石混合介質(zhì)，并通過有限元方法研究了土石混合介質(zhì)的力學(xué)特性．徐文杰等[6]通過數(shù)字圖像方法及大型直剪試驗分析了塊石含量與抗剪強(qiáng)度的關(guān)系．丁秀麗等[7]對土石混合介質(zhì)進(jìn)行灰度處理和邊界重構(gòu)，進(jìn)而建立了土石混合介質(zhì)的顆粒流模型，并對模型進(jìn)行了雙軸壓縮試驗．金磊[8]等通過CT掃描不規(guī)則顆粒生成三維離散元模型，進(jìn)行大三軸試驗數(shù)值模擬研究．

目前對于土石混合介質(zhì)的研究多基于隨機(jī)重構(gòu)，不能反映現(xiàn)場的實際土石分布狀況，另外，基于數(shù)字圖像分析方法建立研究模型，可以避免由于室內(nèi)試驗采樣所造成的土樣擾動等影響，可以較真實再現(xiàn)原位試驗下土石混合介質(zhì)的力學(xué)特性．對于現(xiàn)場拍攝的數(shù)碼照片，在預(yù)處理的基礎(chǔ)上采用區(qū)域生長算法進(jìn)行數(shù)字圖像識別，實現(xiàn)土石混合介質(zhì)的二值化．在此基礎(chǔ)上生成土石混合介質(zhì)的數(shù)值模型，進(jìn)行數(shù)值試驗，進(jìn)一步研究其強(qiáng)度特性及顆粒間的接觸關(guān)系．

1 數(shù)字圖像原理

在計算機(jī)中，數(shù)字圖像[9-10]是由一個個的像素點構(gòu)成的，像素點的值記錄該像素點的色彩、亮度和對比度等圖像信息．因此可以根據(jù)像素點的信息了解整個圖像的信息，這些像素點對應(yīng)的各個離散數(shù)據(jù)便成了數(shù)字圖像處理的基礎(chǔ)．

將一副數(shù)字圖像看作由N×M的像素點組成的矩陣，每個像素點的灰度值對應(yīng)到其所在的坐標(biāo)，以二維函數(shù)f(i，j)表示，其中i，j為空間坐標(biāo)．

(1)

式中，N與M分別為圖像中水平與垂直方向的像素數(shù)量．

2 數(shù)字圖像識別過程

2.1 圖像預(yù)處理

如圖1(a)所示，將現(xiàn)場拍攝的彩色圖像進(jìn)行灰度化，分析其特征線上的灰度分布曲線．如圖1(b)所示，原始圖像中由于受光照條件、土石混合介質(zhì)內(nèi)部礦物質(zhì)色彩差異等因素的影響，存在大量噪聲，難以選取適當(dāng)?shù)淖R別準(zhǔn)則將內(nèi)部的土石區(qū)分．因此，必須對原始圖像進(jìn)行適當(dāng)?shù)念A(yù)處理，增大土石混合介質(zhì)內(nèi)部土和石塊的差異，便于土石閾值的選?。?/p>

圖1 土石混合介質(zhì)灰度圖像及其特征線的灰度分布曲線

將原始圖像首先經(jīng)過中值濾波[11]，去除圖像中的椒鹽噪聲，再通過聚類方法[12-13]將像素點進(jìn)行分類，將圖像中的目標(biāo)區(qū)域和背景區(qū)域互不重疊的區(qū)分開來，最后適當(dāng)調(diào)節(jié)圖像的對比度，使土石差異更加明顯，預(yù)處理后圖像如圖2(a)所示．

2.2 土石閾值的確定

在土石混合介質(zhì)圖像識別中，需要確定一個土石閾值來區(qū)分土和石塊，閾值T的選擇很重要，它是區(qū)分土石邊界的關(guān)鍵問題．閾值確定函數(shù)可以寫為一般形式：

T=T[x，y，p(x，y)，f(x，y)]

(2)

式中，f(x，y)表示灰度值；p(x，y)為某種圖像性質(zhì)，即閾值T一般是p(x，y)，f(x，y)，(x，y)的函數(shù)．圖像的灰度直方圖可以反映一幅圖像中的灰度分布信息，是閾值選取的重要參考依據(jù)．如圖2(b)可見預(yù)處理圖片中特征線上的灰度直方圖中出現(xiàn)明顯的峰段和谷段特征，分別對應(yīng)圖2(a)中的塊石和土，以圖像特征線的灰度分布直方圖為依據(jù)，選取區(qū)域生長算法中土石介質(zhì)的閾值．

圖2 預(yù)處理圖像及其特征線的灰度分布

2.3 區(qū)域生長算法

土石混合介質(zhì)是由具有不同特征的土和石塊構(gòu)成的，土和石塊有著相似的內(nèi)部特征，而相鄰的土與石塊具有不同的特征，即在一個區(qū)域內(nèi)部滿足某種特征的一致性，而相鄰的不同區(qū)域具有不同的特性．區(qū)域生長算法[14-16]即基于上述特點的區(qū)域分割方法，其基本思想是將某范圍內(nèi)具有某些相似性質(zhì)的像素合并起來作為一個“區(qū)域”，首先需要尋找一個種子點，如將石塊的中心點作為區(qū)域生長的種子點，根據(jù)一定的生長規(guī)則，將種子點周圍鄰域的與種子點具有相同性質(zhì)的點合并到種子點的區(qū)域．然后以新點作為種子點，繼續(xù)以上生長過程，直到該區(qū)域內(nèi)沒有滿足要求的點為止．區(qū)域生長算法是一種串行區(qū)域分割的圖像識別方法，其優(yōu)點是基本思想相對簡單，通常能將具有相同特征的聯(lián)通區(qū)域分割出來，并能提供很好的邊界信息和分割結(jié)果．

設(shè)有一個數(shù)字圖像，以圖3(a)中灰度最大值點作為種子，該點的灰度為9．假設(shè)采用的生長相似準(zhǔn)則為鄰點的灰度與已生成的區(qū)域的平均灰度的差小于2，第一次區(qū)域生長得到3個灰度為8的鄰點，如圖3(b)所示．此時4個點已接受點的平均灰度(8+8+8+9)/4=8.25，故第二次區(qū)域生長只得到灰度為7的一個鄰點，如圖3(c)．此時，這5個已接受點的平均灰度為(8+8+8+9+7)/5=8，已不存在灰度大于6的鄰點，所以生長過程終止．

圖3 區(qū)域生長法實現(xiàn)過程

在生長過程中，為了識別結(jié)果更為準(zhǔn)確，采用人工選點方法，作為區(qū)域生長的種子點，即生長起點．根據(jù)第二節(jié)中確定的土石閾值，將區(qū)分土石的灰度值設(shè)置為160．在生長過程中，首先判斷種子點是否為圖片邊界點，然后判斷種子點周圍8個點的灰度值，將符合生長規(guī)則的點設(shè)置為黑場，并且作為下一步生長的起點．最后將所有滿足要求的點納入同一區(qū)域，生成土石混合介質(zhì)二值圖像．

土石混合介質(zhì)原始圖像如圖4所示，首先利用2.1節(jié)的方法對其進(jìn)行預(yù)處理，然后利用區(qū)域生長算法進(jìn)行土石介質(zhì)識別，識別結(jié)果如圖5所示．可見，預(yù)處理在很大程度解決了數(shù)碼圖像中的噪點、土石邊界模糊以及塊石內(nèi)部礦物色彩差異等問題，直接使用區(qū)域生長算法對數(shù)字圖像進(jìn)行識別，得到的二值圖像中已有明顯的土石分布特征．最后，對于圖中仍存在由于石塊邊緣處陰影造成的石塊周圍的孤點、偽點進(jìn)行局部處理，如圖6所示．

圖4 原始圖像 圖5 二值圖像

圖6 局部處理后的二值圖像

3 土石混合介質(zhì)離散元模型

前文已根據(jù)區(qū)域生長算法對土石混合介質(zhì)的原始圖像進(jìn)行了分割處理，在此基礎(chǔ)上可建立離散元模型[17-18]，通過離散元方法對土石混合介質(zhì)在直剪試驗下的力學(xué)特性進(jìn)行研究．

3.1 離散元模型的建立

現(xiàn)場拍攝的數(shù)字圖像是由一個個像素點組成的，若要生成數(shù)值模型，還需要確定實際尺寸與像素之間的比例關(guān)系，若S為數(shù)字圖像中每個像素單位所代表的實際尺寸，L為數(shù)字圖像在橫向或縱向的實際尺寸，N為該方向上的像素點個數(shù)，則轉(zhuǎn)換比例為S=L/N．現(xiàn)場拍攝的土石混合介質(zhì)數(shù)字圖像為0.3 m×0.3 m的正方形圖像，其像素尺寸為1 000Pixel×1 000 Pixel，即轉(zhuǎn)換比例為S=L/N=0.000 3 m/Pixel．

為分析土石混合介質(zhì)的力學(xué)效應(yīng)，基于離散元方法對上述模型進(jìn)行直剪試驗．如圖7所示，模型尺寸為0.3 m×0.3 m，像素尺寸1 000 Pixel×1 000 Pixel，像素單元數(shù)為1 000 000．在該模型范圍內(nèi)生成孔隙率16%，半徑在0.000 5 m到0.001 m之間的顆粒共41 244個．試樣直接剪切試驗的模擬邊界條件與加載條件如圖，上下剪切盒的高度均為0.15 m，根據(jù)室內(nèi)直剪試驗的規(guī)格，在試驗過程中保持上剪切盒不動，推動下剪切盒，利用伺服加載機(jī)制使其正應(yīng)力保持不變．建立墻體作為剪切盒的外墻模型，在試驗中認(rèn)為墻體是剛體，即將墻體的剛度設(shè)置得遠(yuǎn)大于土石介質(zhì)的剛度．當(dāng)剪切位移達(dá)到0.012 m時，終止試驗．試驗中每組試樣法向應(yīng)力分別取200、400以及800 kPa．

圖7 土石混合介質(zhì)離散元模型

3.2 土石介質(zhì)細(xì)觀參數(shù)的確定

根據(jù)顆粒流方法進(jìn)行數(shù)值模擬首先需要進(jìn)行細(xì)觀參數(shù)的標(biāo)定，即探究模型宏觀參數(shù)與細(xì)觀參數(shù)的關(guān)系來確定細(xì)觀參數(shù)．目前對巖石材料宏細(xì)觀參數(shù)關(guān)系基本認(rèn)識如下[19]：1)細(xì)觀粘結(jié)模量控制宏觀楊氏模量(E)；2)k*影響彈性變形階段的泊松比；3)法向與切向粘結(jié)強(qiáng)度比影響試樣的破壞模式．

首先，采用clump來模擬碎石，在每次循環(huán)都不考慮構(gòu)成clump顆粒之間的相互作用，同時若一個接觸中兩個顆粒的細(xì)觀參數(shù)不同時，在計算時選取兩個參數(shù)中最小的一個，但顆粒的接觸剛度除外．根據(jù)石崇[20]提出的參數(shù)標(biāo)定方法，采用雙軸試驗對土和石塊的細(xì)觀參數(shù)進(jìn)行標(biāo)定，最終確定土石混合介質(zhì)離散元模擬參數(shù)見表2．

表1 土石混合介質(zhì)力學(xué)參數(shù)

表2 離散元模型參數(shù)

4 直剪試驗下土石混合介質(zhì)力學(xué)特性分析

4.1 剪切應(yīng)力-位移曲線

如圖8所示的不同法向應(yīng)力下土石混合介質(zhì)的剪切應(yīng)力-位移曲線，土石混合介質(zhì)的剪切應(yīng)力-位移曲線整體分為峰前峰后兩個階段．

圖8 不同法向應(yīng)力下土石混合介質(zhì)剪切應(yīng)力-位移關(guān)系曲線

峰前階段，隨著剪切位移的增大，剪切應(yīng)力也在逐漸增大，剪切應(yīng)力和剪切位移呈現(xiàn)出非線性的變化趨勢．隨著法向應(yīng)力的增加，剪切應(yīng)力峰值也在逐漸增加．峰后階段，隨著剪切位移的增大，剪切應(yīng)力先減小，經(jīng)歷一定剪切位移后又繼續(xù)增大，這是由于在峰值處，土石混合介質(zhì)中的顆粒旋轉(zhuǎn)到一定位置，顆粒間的初始咬合斷裂，造成剪切應(yīng)力下降，隨著剪切位移的增加，土石介質(zhì)間形成新的咬合，導(dǎo)致剪應(yīng)力升高．對于同一試樣，法向應(yīng)力越大，這種新的咬合作用表現(xiàn)得越明顯，因此，隨著法向應(yīng)力的增加，剪切應(yīng)力峰值也在不斷增加．另外，隨著剪切位移的增大，模型中逐漸形成不規(guī)則的剪切帶．相同剪切位移下，法向應(yīng)力越大，對應(yīng)的剪切帶也越大，即塊石間的相互作用隨著法向應(yīng)力的增大而增大．

4.2 顆粒間接觸力的演化過程

在直剪過程中，為了抵抗剪應(yīng)力，試樣內(nèi)部顆粒間接觸力的大小和方向會發(fā)生一定的調(diào)整，通過分析剪切過程中試樣內(nèi)部力鏈網(wǎng)絡(luò)和接觸組構(gòu)的演化過程，可以清楚地了解剪切過程中土石顆粒間接觸力的相互傳遞演化過程．

圖9為法向應(yīng)力200 kPa下不同剪切位移下的力鏈網(wǎng)絡(luò)分布規(guī)律，其中，圖9(a)～(d)剪切位移分別為0，1，2，12 mm．可見，伺服完成后，剪切盒內(nèi)整體穩(wěn)定，力鏈在剪切盒內(nèi)均勻分布，表明此時剪應(yīng)力是由土體和石塊兩者共同來承擔(dān)的，且石塊與土體顆粒所承擔(dān)的剪應(yīng)力大小基本相同．隨著剪切的開始，由于受剪切帶處剪切作用，該處產(chǎn)生大量強(qiáng)力鏈而其他位置的力鏈多為弱力鏈．在峰值應(yīng)力狀態(tài)前，剪切盒中右上和左下部分的力鏈分布發(fā)生明顯變化，而試樣中部的力鏈仍保持豎直．在峰值應(yīng)力狀態(tài)處，試樣中部力鏈的分布特征發(fā)生完全改變，強(qiáng)力鏈分布在右上、左下以及剪切帶位置，其中剪切帶位置的力鏈逐漸沿30°至45°方向傾斜．隨著剪切的繼續(xù)進(jìn)行，強(qiáng)力鏈進(jìn)一步增強(qiáng)，但其方向基本保持30°至45°方向．

圖9 法向應(yīng)力200 kPa下不同剪切位移時的力鏈網(wǎng)絡(luò)分布規(guī)律

圖10為法向應(yīng)力200 kPa不同剪切位移下的接觸力分布規(guī)律，其中圖10(a)～(d)剪切位移分別為0，1，2，12 mm．可見，初始階段的各個方向接觸力的大小分布比較均勻，分布形狀類似于一個圓，隨著剪切位移的增加，為了抵抗剪切變形，顆粒間接觸力發(fā)生明顯變化，分布形狀由圓形逐漸變成橢圓形，且橢圓的長軸方向與剪切面形成一定角度．在到達(dá)峰值應(yīng)力時，橢圓形狀變得扁平，試樣內(nèi)部顆粒間的接觸力沿著抵抗剪切方向逐漸增大．當(dāng)剪切應(yīng)力到達(dá)峰值之后，隨著剪切位移的增加，試樣內(nèi)部顆粒間的接觸力沿著抵抗剪切方向也將逐漸減小，導(dǎo)致橢圓的長短軸之比減小，橢圓形狀明顯變得粗壯．隨著剪切位移的增加，剪切應(yīng)力大小基本不變，試樣內(nèi)部顆粒間的接觸力也將保持不變．

圖10 法向應(yīng)力200 kPa下不同剪切位移時的接觸分布規(guī)律

5 結(jié) 論

采用土石混合介質(zhì)數(shù)字圖像，建立了土石混合介質(zhì)的離散元模型，并對土石混合介質(zhì)進(jìn)行了二維直剪試驗，得到了以下結(jié)論：

1)基于數(shù)字圖像識別的土石混合介質(zhì)離散元模型開展直剪試驗數(shù)值模擬，可保證土石混合介質(zhì)中的塊石分布情況與實際情況一致．

2)應(yīng)用區(qū)域生長算法可以實現(xiàn)土石混合介質(zhì)的分割，可將具有相同特征的土和石塊識別出來，并且保留清晰的邊界輪廓信息和分割效果，克服了簡單灰度分割方法中可能出現(xiàn)的欠分割和過度分割等缺陷．

3)通過顆粒間的組構(gòu)分布和力鏈傳遞關(guān)系直接反映試樣內(nèi)部的顆粒間接觸以及接觸力的變化情況，可作為現(xiàn)場試驗的補(bǔ)充，從微觀角度反映土石混合介質(zhì)在直剪試驗下的力學(xué)性質(zhì)．