劉維福，丁序海，周 光

（神華國(guó)能集團(tuán)有限公司 煤炭管理部，北京100033）

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

1 實(shí)驗(yàn)煤樣及實(shí)驗(yàn)方法

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

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

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

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

2.1 圖像中值濾波處理

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

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

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

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 中值濾波對(duì)原圖像進(jìn)行濾除噪聲處理，通過(guò)中值濾波法保護(hù)了圖像的原有信息，濾除了原始圖像的噪聲點(diǎn)位，使得處理后的圖像信息更加準(zhǔn)確，并作為后續(xù)二值化處理的圖像基礎(chǔ)。

2.2 圖像二值化處理

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

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

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

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

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

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

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

2.4 CT 圖像的孔隙率計(jì)算

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

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

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

由此可以得到圖像中每1 個(gè)像素點(diǎn)的孔隙大小，煤樣每層孔隙率如圖4。

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

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

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

3.1 孔隙度的分布

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

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

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

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

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

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

在圖6（a）中水流速度從0.017 5 m/s 初始速度衰減到出口速度0.017 2 m/s。這主要是因?yàn)樵诖嗣簶又髁严吨?，裂隙發(fā)育較為平緩，水在流動(dòng)過(guò)程中受到的阻力較小。在圖6（b）中由于主裂隙發(fā)育的較為完好且彎折不大，水在流動(dòng)過(guò)程中受到阻力變大，在這一過(guò)程中水的流速會(huì)出現(xiàn)劇減的狀態(tài)，之后速度趨于穩(wěn)定減小。水在經(jīng)過(guò)次生裂隙的過(guò)程中水流動(dòng)速度會(huì)出現(xiàn)局部增大的趨勢(shì)，這主要是因?yàn)樗诮?jīng)過(guò)次生裂隙的過(guò)程中，因?yàn)榇紊严兜目卓诎霃姜M小且彎曲曲折方式急劇，水會(huì)出現(xiàn)流動(dòng)速度激增的狀態(tài)。雖然煤體中的次生裂隙很多，但水的流向不會(huì)通過(guò)所有的次生裂隙。這主要是因?yàn)樗牧鲃?dòng)具有選擇性，水會(huì)選擇從孔隙結(jié)構(gòu)發(fā)育完好的地方滲透。在次生裂隙和基質(zhì)夾雜部分流動(dòng)不是很明顯。

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

4 結(jié) 論

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

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

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