基于CT 成像的煤巖孔裂隙結(jié)構(gòu)重建及滲流模擬研究

劉維福，丁序海，周 光

（神華國能集團有限公司 煤炭管理部，北京100033）

隨著井工煤炭生產(chǎn)開采的需要[1]，地下礦井高效挖掘開采過程截割齒截割破碎過程總是伴隨著大量煤塵的產(chǎn)生[2]，長期暴露于工作面上高濃度粉塵條件下的煤礦工人及易罹患嚴(yán)重肺塵?。–WP）[3]，高濃度粉塵也是誘發(fā)粉塵爆炸災(zāi)害的主要原因，嚴(yán)重制約井下煤炭安全生產(chǎn)[4]。煤塵注水作為井下高效防塵防突、驅(qū)排瓦斯主要方法[5]，已被越來越多的應(yīng)用于生產(chǎn)過程中[6]。在煤層內(nèi)存在孔隙和裂隙2個系統(tǒng)[7]：孔隙系統(tǒng)發(fā)育在煤巖基塊中，天然裂隙是煤中流體滲透通道。煤的孔隙結(jié)構(gòu)和形態(tài)直接關(guān)系著孔隙率、滲透率。連通性是孔隙的重要特征，它直接關(guān)系到儲層的滲透作用。目前絕大多數(shù)的學(xué)者主要聚焦于煤的內(nèi)部孔隙裂隙分布、CT 掃描復(fù)原上，關(guān)于煤的細(xì)觀表征建模、水在煤體中滲流流動規(guī)律上的研究很少。藺麗君[8]等通過COMSOL 構(gòu)建一種煤體二維孔隙尺度流動與變形數(shù)值模型，分析了在應(yīng)力作用下巖石孔隙結(jié)構(gòu)演化趨勢及對滲流影響。郎穎嫻[9]等通過RFPA 建立了可以反映玄武巖內(nèi)部孔隙結(jié)構(gòu)的巴西盤數(shù)值模型，研究了孔隙對試樣的破壞機制和抗拉強度的影響。李相臣[10]運用FDK 算法重建煤巖的灰度圖像，對煤巖中的孔隙尺度以及空間分布進行綜合表征。白若男[11]等通過CT 掃描，重構(gòu)了煤體模型，分析了流場的流線、壓力、速度及滲流方向的演化規(guī)律。陳俊國[12]等通過圖像分析技術(shù)，分析了滲透率隨圍壓的變化規(guī)律。為了準(zhǔn)確描述煤層注水時煤樣滲流過程，通過對煤體進行CT 掃描，經(jīng)過MATLAB 對掃描圖像進行中值濾波處理，對孔隙變化準(zhǔn)確描述，并運用數(shù)值模擬軟件對煤體進行了重構(gòu)，分析研究水在煤體流動的機理及運動規(guī)律，為煤層注水降塵工藝提供指導(dǎo)意義。

1 實驗煤樣及實驗方法

研究煤樣取自阜新市雷家區(qū)4 號煤層，煤樣種類為長焰煤。按照GB/T 482—2008 國家標(biāo)準(zhǔn)進行采樣。在實驗室用標(biāo)準(zhǔn)煤樣切割機將煤樣切割成5 cm×5 cm×10 cm 的長方體，獲得3 個標(biāo)準(zhǔn)煤樣，選取1 個裂隙發(fā)育完好煤樣進行探究。

工業(yè)CT 采用先進的高頻恒壓X 射線源、數(shù)字圖像探測器以及高精度機械檢測平臺。工業(yè)CT 原理是利用X 射線穿過被檢測物體的截面，提取衰減后的信息。被掃描物體的密度由X 射線衰減系數(shù)來體現(xiàn)。物體的密度越大，對X 射線吸收能力越強。通過衰減系數(shù)轉(zhuǎn)換成CT 數(shù)，CT 數(shù)越大相當(dāng)于物體密度越高。通過對比CT 圖像中的黑白分布展現(xiàn)出被測物體的密度分布情況。通過CT 掃描后的物體，可以真實無損地將煤體內(nèi)部孔隙結(jié)構(gòu)表征出。

采用丹東奧龍射線儀器有限公司AL-CT-225工業(yè)X 射線CT 檢測系統(tǒng)，采用錐束掃描和DR 實時成像多種檢測方式，可掃描具有缺陷、孔隙分析和被檢測工件區(qū)掃功能。層析識別系統(tǒng)的CT 空間最小分辨率為6 μm，密度分辨率為0.5%，最大成像素為1 024×1 024，實驗共掃描出800 個切片圖像。

2 CT 圖像數(shù)據(jù)處理

2.1 圖像中值濾波處理

采用MATLAB 數(shù)字圖像技術(shù)對CT 掃描圖像中的孔隙結(jié)構(gòu)信息進行識別處理，將煤巖的孔隙裂隙和基質(zhì)進行分割。利用X 射線CT 掃描得到的圖像存在系統(tǒng)噪聲，噪聲定義為不可預(yù)測，只能用概率統(tǒng)計方法來認(rèn)識的隨機誤差。即噪聲在圖像上是隨機分布的，這些噪聲點位在圖像上會增加或減小圖像像素點位的真實像素值，常表現(xiàn)為孤立像素點或像素塊。所以降低了圖像的質(zhì)量，通過相關(guān)的濾波函數(shù)處理噪聲，選取中值濾波進行處理。

中值濾波對脈沖噪聲有良好的濾除作用，在濾除噪聲的同時還可以保護信號的邊緣，使之不被模糊。因此，為了將每個點的信息準(zhǔn)確提取，將圖片中噪聲點使用中值濾波函數(shù)進行處理，并利用MATLAB 將圖像上存在的每個像素點的值導(dǎo)出為矩陣的模式。中值濾波的基本原理為把數(shù)字圖像或數(shù)字序列中的1 點值，用其鄰域中各點值進行代替。對1 個數(shù)字信號序列進行濾波處理時，首先定義1 個長度L 為奇數(shù)的窗口，L＝2n＋1。窗口信號樣本為（x1，x2，x3，…，xn），將n 個數(shù)據(jù)大小進行排列：x1＜x2＜x3＜…＜xn。則有：

其中y 為窗口序列（x1，x2，x3，…，xn）的中值，將窗口在數(shù)據(jù)上滑動，將窗口正中所對的像素值用窗口各個像素的中值代替，即：

y（i）＝Med［x（i－N），…，x（i），…，x（i＋N）］

MATLAB 中值濾波算法前后原始圖像和中值濾波濾除噪聲比較如圖1。

圖1 原始圖像和中值濾波濾除噪聲比較Fig.1 Comparison of noise removal between original images and median filtering

利用MATLAB 中值濾波對原圖像進行濾除噪聲處理，通過中值濾波法保護了圖像的原有信息，濾除了原始圖像的噪聲點位，使得處理后的圖像信息更加準(zhǔn)確，并作為后續(xù)二值化處理的圖像基礎(chǔ)。

2.2 圖像二值化處理

二值圖像是1 個取值只有0 和1 的邏輯類型的二維數(shù)組。所以二值圖像可以看成是一種特殊的亮度圖像。1 種只有2 個亮度級（最暗和最亮）的亮度圖像。CT 掃描后的圖片為灰度圖像，其中的微小裂隙較難提取。選取OStu 法選確定最佳閾值，通過處理的圖像可以提取圖像的孔隙結(jié)構(gòu)信息特征。二值化處理的圖像如圖2。

圖2 二值化處理的圖像Fig.2 Binary images processing

原始灰度圖像經(jīng)MATLAB 的中值濾波處理后，去除了圖像的噪聲點位。清晰地分割了煤體中的孔隙和基質(zhì)。有效的保留了煤體的孔隙特征和邊界信息。而通過二值化處理后的原始圖像清晰地分割了煤體孔隙和基質(zhì)。為下一步的二維模型建模提供了基礎(chǔ)，并以此作為有限元方法構(gòu)建數(shù)值模型。

2.3 煤體微觀精細(xì)模型建立

選取了煤體試樣2 個孔隙發(fā)育完好煤體圖片，煤體微觀孔隙結(jié)構(gòu)模型如圖3。由圖3 可以看出，主裂隙貫穿整個煤樣，其中次生裂隙在煤體中分布。部分次生裂隙和主干裂隙相互連接，孔隙在試樣內(nèi)的分布均勻。

圖3 煤體微觀孔隙結(jié)構(gòu)模型Fig.3 Microscopic pore structure models of coal

將中值濾波、二值化處理后的圖像信息進行提取，運用COMSOL 中的插值函數(shù)將煤樣切片建模，以此為數(shù)值模擬的模型基礎(chǔ)。

2.4 CT 圖像的孔隙率計算

孔隙率指散粒狀材料表現(xiàn)體積中材料內(nèi)部的孔隙占總體積的比例。CT 掃描圖像中包括了煤體的骨架孔隙信息特征，其中包含了主裂隙和次裂隙可以將煤體內(nèi)部的孔隙結(jié)構(gòu)展示出來。

Taud H[13]在CT 掃描圖像內(nèi)，把CT 灰度級類比成為三維地形圖的高低，灰度水平等效于三維地形圖的地形高低?；叶戎翟酱蟮狞c位地形圖所在位置越低，結(jié)合數(shù)理統(tǒng)計方法獲得孔隙率的分布函數(shù)，通過灰度值的大小和統(tǒng)計灰度值出現(xiàn)的次數(shù)可以計算出試樣的孔隙率，孔隙率φ 公式為：

式中：ri為各像素點的灰度值，取值范圍［rmin，rmax］；H（ri）為［rmin，rmax］范圍內(nèi)的灰度直方圖；rmin為圖像灰度最小值；rmax為圖像灰度最大值。

由此可以得到圖像中每1 個像素點的孔隙大小，煤樣每層孔隙率如圖4。

圖4 煤樣每層孔隙率Fig.4 Porosity of each layer of sample coal

圖4 中試樣煤樣的孔隙率平均值為14.82%，試樣煤樣每個切片層數(shù)的孔隙率波動變化在0.02 范圍內(nèi)。說明試樣內(nèi)孔隙數(shù)據(jù)具有一定的可靠性。

3 煤體內(nèi)部水滲流通道的數(shù)值模擬

3.1 孔隙度的分布

通過對比經(jīng)過中值濾波、二值化、COMSOL 建模的模型，可以看出煤體內(nèi)部的主裂隙分布明顯，在中值濾波和二值化處理后的圖像，煤的裂隙分布在黑色線條內(nèi)。通過COMSOL 建模后得到的孔隙度分布圖如圖5?？梢钥闯觯后w的主裂隙呈現(xiàn)紅色區(qū)域，而藍色區(qū)域代表的是煤基質(zhì)。其他區(qū)域代表的是次生裂隙和煤基質(zhì)夾雜部分。圖5（a）中平均孔隙度為14.37%，圖5（b）中平均孔隙度為15.26%。

圖5 孔隙度分布圖Fig.5 Porosity distribution map

3.2 水在煤體中流動的數(shù)值模擬

煤體內(nèi)部的孔隙結(jié)構(gòu)錯綜復(fù)雜，這也就導(dǎo)致了水在煤體中運動過程的復(fù)雜性。流體介質(zhì)在多孔介質(zhì)流動中流動呈向的是多孔流動流出，當(dāng)水在煤體中運動過程中會通過主孔隙流動，分裂隙連同基質(zhì)處水部分流動?？紤]到水流動在煤體中的復(fù)雜性，現(xiàn)簡化水在煤體的流動過程。煤體的上邊界設(shè)置20 mm 的水頭，在下邊界設(shè)置為0 的水頭。這樣產(chǎn)生了20 mm 的水頭差，驅(qū)動流體流動,左右邊界不可滲透，流體密度為1 000 kg/m3,動力黏度0.001 Pa·s。滲流水模擬結(jié)果圖如圖6。

圖6 滲流水模擬結(jié)果圖Fig.6 Diagram of seepage water simulation results

從圖6 可以看出，煤體含有主裂隙和次生裂隙，水在經(jīng)過煤體的過程中，在主裂隙流動過程明顯，水也會在次生裂隙中流動。水在流動過程中，速度會呈現(xiàn)逐漸減小后的趨勢。水在主裂隙流動過程中，速度變化的趨勢不大。水流在主裂隙水流流速如圖7。

在圖6（a）中水流速度從0.017 5 m/s 初始速度衰減到出口速度0.017 2 m/s。這主要是因為在此煤樣主裂隙中，裂隙發(fā)育較為平緩，水在流動過程中受到的阻力較小。在圖6（b）中由于主裂隙發(fā)育的較為完好且彎折不大，水在流動過程中受到阻力變大，在這一過程中水的流速會出現(xiàn)劇減的狀態(tài)，之后速度趨于穩(wěn)定減小。水在經(jīng)過次生裂隙的過程中水流動速度會出現(xiàn)局部增大的趨勢，這主要是因為水在經(jīng)過次生裂隙的過程中，因為次生裂隙的孔口半徑狹小且彎曲曲折方式急劇，水會出現(xiàn)流動速度激增的狀態(tài)。雖然煤體中的次生裂隙很多，但水的流向不會通過所有的次生裂隙。這主要是因為水的流動具有選擇性，水會選擇從孔隙結(jié)構(gòu)發(fā)育完好的地方滲透。在次生裂隙和基質(zhì)夾雜部分流動不是很明顯。

圖7 主裂隙水流流速圖Fig.7 Velocity diagram of main fractures

4 結(jié) 論

1）通過CT 掃描煤體圖像，通過MATLAB 中值濾波對圖像進行處理進而得到圖像的每一點灰度值，并以二值化處理得到了煤體的孔隙結(jié)構(gòu)。運用COMSOL 插值函數(shù)將煤體進行二維模型重構(gòu)，將煤體上的孔裂隙分布特征進行了精細(xì)表征。

2）運用Taud H 的數(shù)理統(tǒng)計方法獲得了孔隙率的分布函數(shù)，對比COMSOL 模擬做出了孔隙度的分布情況，兩相對比反演了煤的孔隙率。運用數(shù)理統(tǒng)計方法得到煤體的平均孔隙率為14.82%，在COMSOL中模擬數(shù)值的結(jié)果為14.37%、15.26%。證明了這種建模的可靠性和精準(zhǔn)性。

3）模擬結(jié)果表明，水在煤體流動過程中，裂隙的發(fā)育程度會影響水在煤體中的運動。整體的趨勢為水流速呈現(xiàn)減小的趨勢。沿著主裂隙流動的流速衰減趨勢很小，而在次裂隙流動的水跟裂隙的發(fā)育完好程度以及孔裂隙的彎曲程度有關(guān)。在煤體基質(zhì)的流動速度基本為0。