凍融巖石CT圖像處理及細觀特征分析

張慧梅，袁 超，慕娜娜，張 嬋，路亞妮

(1.西安科技大學 理學院，陜西 西安 710054；2.西安科技大學 建筑與土木工程學院，陜西 西安 710054；3.湖北工程學院 土木工程學院，湖北 孝感 432000)

0 引 言

“一帶一路”倡議的提出，為中國煤礦行業提供了新的發展契機，隨著開采水平的提高，開采深度不斷增加[1]。相應地，季節性凍融的影響深度和影響范圍也隨之增加，對巷道工程施工，尤其是砂巖巷道的支護提出了更高的要求。因此，開展凍融循環作用下煤礦砂巖的力學特性研究對指導西部寒區煤礦安全開采至關重要。

巖石凍融破壞為寒區工程建設中的常遇災害，針對此問題，國內外學者進行較為全面的研究，如徐光苗等通過進行力學特性試驗，分析不同宏觀力學參數與凍融循環次數之間的關系[2]；張慧梅等對巖石在凍融循環作用下的損傷劣化規律進行研究，提出凍融荷載耦合的損傷模型[3-4]。任曉龍等對四川芙蓉白膠煤礦的煤樣進行不同溫度下的單軸壓縮試驗，研究溫度對煤巖物理力學性質的影響[5]。LUO等通過以輝綠巖為研究對象，進行不同次數的凍融循環試驗，得到其質量和P波波速隨凍融作用的變化規律[6]。KHANLARI等以伊朗的典型紅砂巖為研究對象，進行室內以及現場的凍融循環試驗，分析其物理力學性質隨凍融次數的變化情況[7]。ZHANG等通過真實紅砂巖預制裂隙的方法，研究不同裂隙幾何特征對巖體強度及裂紋擴展演化規律的影響[8]。

但僅從宏觀角度出發研究，難以解釋試樣在凍融過程中產生的層狀破裂和片落破壞[9]。為揭示凍融破壞的細觀機理，細觀尺度的研究受到研究人員的重視，如翟成等利用核磁共振，分析凍融作用下煤體孔隙結構特征的變化規律[10]；李杰林等進行核磁共振試驗，研究了巖石在凍融作用下內部孔隙結構的損傷特性，分析核磁共振圖譜與凍融循環次數之間的關系[11]。聞名等基于SEM圖像，分析水分及凍融環境對巖石動態抗壓強度的影響[12]。在眾多細觀試驗技術中，計算機層析掃描技術(CT)，因其可視化程度高，使用最為廣泛，國內楊更社較早地開展相關試驗，并提出了一種采用CT數描述的凍融損傷模型，實現對凍融損傷的定量描述[13]；劉慧等運用CT圖像直方圖技術，定量分析巖石在凍結情況下的損傷信息[14]。國際上，ARGANDONA等利用CT檢測技術對凍融巖樣進行斷面掃描，研究不同凍融條件下巖樣的孔隙結構變化規律[15]；PARK和KOCK等則分別利用CT技術研究火成巖和沉積巖在凍融循環作用下孔隙擴展過程[16-17]。

隨著巖石細觀試驗技術的成熟和試驗裝置的普及，合理的圖像后處理成為決定試驗的準確性關鍵因素，王超等提出一種基于LBP和GLCM的煤巖圖像特征提取方法，可以實現巖石表面幾何特征的快速識別[18]。針對CT圖像，張嘉凡等基于改進的K均值聚類算法對砂巖CT圖像進行閾值分割，實現凍融巖石橫截面孔隙分布的有效識別[19]；張慧梅等則考慮孔隙結構的縱向空間分布，實現凍融巖石孔隙分布的三維重構[20]。

常用的圖像處理技術多數以傳統歐氏幾何為基礎，僅通過二維孔隙面積或者三維孔隙體積來描述巖石的凍融損傷程度，無法對巖石孔隙的結構變化特征進行描述，巖石分形理論的提出為砂巖孔隙特征定量描述提供新的思路[21]，史宏財對煤樣進行各個溫度預損傷和三軸聲發射試驗，探討高溫預損傷后煤巖的力學行為[22]。

張志鎮等利用壓汞法研究高溫處理后花崗巖樣品的孔隙特征和分形結構特征[23]。王秀娟等利用分形理論，計算得到煤滲流孔的分形維數，探討分形維數與滲透率的關系[24]。通過引入分形維數來描述圖像復雜性特征的方法受到廣泛認可，該方法也被用于巖石孔隙幾何特征的研究，但目前已有研究多數分析特定CT圖像的較多，分析砂巖整個縱向高度分形特征的較少；定性分析的較多，定量描述的較少。

以新鮮采集的煤礦砂巖為研究對象，通過CT掃描試驗，獲取掃描試樣在不同凍融循環次數下的斷面掃描圖像，經預處理后，得到不同截面分形維數隨凍融循環次數的變化規律。并通過分形維數定義一種描述巖石孔隙結構損傷的新型損傷變量，定量地分析不同凍融循環次數下煤礦砂巖的損傷演化規律。

1 試驗概況

受地理環境的影響，中國礦區分布的煤礦砂巖受凍融災害影響較多，考慮到研究對象的典型性，筆者依托陜西彬縣和長武縣內的煤礦礦井建設項目，采集一批新鮮的煤礦砂巖樣品。試樣制備完成后，初篩剔除外表明顯缺陷的試件，然后利用聲波檢測儀篩選波速相近的試件進行相關試驗，以便于盡可能的排除自身差異性。

1.1 凍融循環試驗

依照《工程巖體試驗方法標準》(GB/T50266—2013)[25]，凍融試驗前先將巖樣置于(105±1)℃的恒溫箱中烘干24 h，然后采用真空飽和儀對巖樣進行0.1 MPa條件下的強制飽水。

凍融循環試驗由圖1所示的XMT605快速凍融試驗機完成，試驗過程中試樣處于無約束狀態，按溫度控制凍融循環時間，巖心溫度由+20 ℃降低到-20 ℃后再重新升溫到+20 ℃為一個凍融循環。為保證試樣中心溫度達到設定值，試驗中設置有控制試樣，內置有高精度熱敏電阻溫度傳感器。

圖1 XMT605凍融循環試驗機Fig.1 XMT605 freeze-thaw cycle testing machine

通過對試驗現象的觀測發現，凍融次數較少時，砂巖試樣整體完好，僅表面有零散的顆粒脫落；凍融次數達到5次時，大部分砂巖試樣仍能保持完整，有少量試樣表面出現淺裂紋；當凍融次數增加到10次時，試樣表面顆粒開始大量脫落，脫落體多為顆粒狀，表現為剝落破壞；當凍融次數達到20次時，試樣的破壞形式變化不大，部分試樣表面開始出現片狀脫落；當凍融次數到達40次時，巖樣的破壞主要仍集中于兩端，主要為嚴重的層裂破壞，小部分出現了貫穿的橫向裂紋，試件沿著裂紋面逐漸滑移成2個部分，典型的破壞形式如圖2所示。

圖2 砂巖試樣典型的破壞模式Fig.2 Typical failure modes of sandstone specimens

1.2 計算機層析掃描(CT)試驗

CT試驗試驗裝置為圖3的YTU225型工業X射線CT機，掃描精度為0.005 mm，掃描照片分辨率為96 DPI。

圖3 YTU225型工業CT Fig.3 YTU225 industrial CT machine

實驗掃描層數預設為900層，由于吸漲作用和端面脫落，實際得到不同凍融次數下的掃描圖像數量有輕微差異，凍融次數n=0，5，10，20，40時第450層的CT圖像，如圖4所示。

圖4 不同凍融次數下第450層的CT圖像Fig.4 CT image of layer 450 under different numbers of freeze-thaw cycle

CT圖像是與巖石密度相關的灰度圖像，它是由0～255個灰度值所構成，灰度值直接反映了物體的密度，圖4中高亮點代表的是密度大的巖石顆粒，暗些的點代表的則是密度較小的孔隙結構。

2 CT圖像預處理

CT圖像雖然能夠更加直觀地體現孔隙變化規律，但由于肉眼分辨能力較低，不能準確察覺灰度圖像中的細微變化，且獲取圖像的過程中會受到來自各個方面的噪聲影響，給后續凍融巖石的損傷劣化定量分析帶來了困難，因此還需進一步對圖像進行預處理，減少噪聲，豐富CT圖像中的力學信息。

2.1 中值濾波

中值濾波就是把圖像中的每個點都用其領域窗口內的中值代替，從而消除孤立的噪聲點，其好處是可以減少圖像中的非線性噪聲，有效地保護邊緣信息，方法簡便，算法也容易實現。圖5為n=40時第450層圖像進行處理前后的對比。

圖5 中值濾波前后對比Fig.5 Comparison before and after median filtering

從圖5可以看出，由于巖石邊緣被其鄰域的中值代替，所以巖石邊緣和顆粒間隙略顯模糊，但前后并沒有太大的灰度差異，仍能保持圖像中重要細節部分。

2.2 圖像直方圖均衡化

圖5(b)的灰度分布直方圖如圖6所示，可以看出圖像的灰分布范圍較窄且極不均勻，直方圖多密集靠近在一起且中間突出一個高峰，在50～150的像素點個數幾乎為0。

圖6 CT圖像灰度直方圖Fig.6 CT Image gray histogram

如若將直方圖分布在整個水平方向上壓縮，展開成一個新直方圖，其清晰程度會提高，所需目標信息會更突出。也就是說，若圖像存在噪聲，直方圖均衡化會放大圖像中的噪聲。因此，若使用直方圖均衡化對圖像進行增強，應先對圖像進行降噪處理。圖7為依次進行中值濾波和直方圖均衡化處理后n=40時第450層CT圖像。

圖7 均衡化后的CT圖像Fig.7 CT image after equalization

對比圖5及圖7可知，經中值濾波和均衡化處理后的CT圖像整體效果更加鮮明，對比度也更加明顯，巖石內部細節部分處理的更好，使得巖石與孔隙對比更明顯，巖石孔隙更加清晰。

2.3 圖像二值化處理

為了進一步觀察巖石顆粒與孔隙的分布特點，本文采用閾值法對CT圖像進行二值化處理，其原理為：先由用戶指定或通過算法生成一個閾值，若圖像中某像素點的灰度值小于該閾值，則將該像素灰度值設置為0或255，否則設置為255或0。

二值化處理的變換函數表達式如下

(1)

式中T為閾值。閾值選擇是灰度圖二值化的關鍵步驟，常見的閾值選取方法有人工選取法、最大類間方差法等。人工選取分割閾值容易受主觀因素影響，而迭代法受圖像目標和背景比例差異影響較多，因此本文采用最大類間方差法對圖像進行二值化處理。圖8為n=40時第450層圖像二值化處理后的圖像。

圖8 二值化處理后的CT圖像Fig.8 CT image after binary operation

二值化處理后，圖中白色代表巖石顆粒，黑色代表巖石孔隙和背景，從圖8可以看出，巖石顆粒與孔隙黑白分明，極大地提高了原始CT圖像的直觀性。預處理完成后凍融次數n=0,5,10,20,40時第450層的二值化CT圖像如圖9所示。

圖9 預處理完成后的CT圖像Fig.9 CT image after all pre-operations

3 基于CT圖像的細觀特征分析

3.1 巖石孔隙率的分析

為了定量分析巖石孔隙大小的變化規律，基于最大類間方差法對不同凍融循環次數下不同截面的孔隙率進行了計算，并給出其不同凍融次數下的孔隙率均值，如圖10所示。

圖10 截面計算孔隙率Fig.10 Computational porosity of different cross-sections

從圖10可以看出，凍融0次時巖石的孔隙率為13.97%，與試驗測得的孔隙率14.3%較為接近，說明了處理方法的可靠性。紅砂巖孔隙率均值隨凍融次數的增加整體呈上升趨勢，但不同凍融次數下試樣孔隙率在軸向的分布比較穩定。然而在凍融10次與20次時，孔隙率增加的并不明顯，且曲線有較長的交錯段。然而同時進行的力學試驗表明[4]，凍融10次與20次試樣的彈性模量及抗壓強度有明顯差異，因此僅對孔隙率大小變化規律進行研究無法準確描述凍融巖石細觀損傷劣化規律。

3.2 巖石分形維數的計算

分形理論作為一種新方法新理論，正在許多領域應用探索，它跳出一維線、二維面、三維體的傳統藩籬，采用連續維度的思維來描述物體的屬性與狀態，進而表達出物體的多樣性與復雜性。分形維數是分形理論及其應用中最重要的概念之一，是描述物體復雜度的重要參數[26]。根據其不同的定義方法可以分為Hausdorff維數DH、相似維數DS、信息維數Di等等。其中，許多等價Hausdorff維數的維數被提出來，例如文中使用的計盒維數。

3.3 計盒維數的計算方法

采用不同半徑的盒子將圖像中的巖石部分全部覆蓋，當選取的盒子尺寸不同，其數量也不相同。設小盒子的半徑為ε，盒子的數量為Nε，選取不同的ε值，假設兩者滿足

Nε=AεB

(2)

式中A,B為參數。

對(2)式變形，可以得到

B=DB+k

(3)

圖像二值化處理完成后，選取不同的ε對圖像分區，在每一種網格劃分下，計算出所有覆蓋圖像中含巖石區域的網格的數目，記為網格覆蓋數，根據ε取值的不同，可得到一系列尺寸不同的“網格”和相應“網格覆蓋數”的數據對，即Nε-ε數據對，然后在雙對數坐標系下畫出數據對，利用最小二乘法擬合為一次函數，得到一次函數的斜率k，其相反數即為圖像的計盒維數。

3.4 計算結果分析

在CT掃描試驗中，經歷0次、5次、10次、20次、40次凍融循環后得到的掃描圖像分別為 884張、899張、775張、925張、905張，由于試件兩端掃描圖像不完整，故選取100～800層的圖像進行計算，在MATLAB軟件中通過FOR循環語句依次讀取文件夾中所保存的CT圖像，按照盒維數的計算原理編寫程序，分別計算出每一張CT圖像的分形維數，最后將其保存到統一的XLSX文件中，進行ORIGIN作圖，結果如圖11所示。

圖11 截面計盒維數Fig.11 Box dimension of different cross-sections

從圖11可以看出，計盒維數在試樣軸向的分布具比較穩定，但存在有小幅震蕩，為了挖掘計盒維數變化的統計學規律，應用移動平均法對圖11中的數據進行平滑處理，移動平均期數取值為50時的結果如圖12所示。

圖12 移動平均處理后不同截面分形維數Fig.12 Box dimensions of different cross sections after moving average operation

從圖11和圖12可以看出，雖然各截面的計盒維數不同，但整體上與凍融循環次數負相關，說明巖石內部不同位置具有不同的孔隙特征，且隨著凍融循環次數的增加，巖石內部復雜的微小孔隙結構逐漸擴展并聚合為分布簡單的較大孔隙。不同凍融次數下的計盒維數平均值如圖13所示。

圖13 孔隙分形維數均值隨凍融次數的變化Fig.13 Variation of mean value of pore fractal dimension with freeze-thaw times

從圖13可以看出，0～5次凍融階段，孔隙分形維數下降速度較慢，僅下降0.18%，孔隙率增加也較小，說明凍融初期巖石內部孔隙發展緩慢，凍融損傷增長不明顯；隨著凍融循環次數的增加，巖石內部損傷開始加速演化，表現為孔隙分形維數的大幅度降低和孔隙率的快速增加，即內部孔隙增長的同時也不斷發生連通和擴展，當凍融次數由5次增長為10次時巖石內部萌生許多新的微孔隙，與此同時原有的微孔隙也不斷地聚集聯通，損傷快速演化；隨后凍融次數由10次增長到20次的階段，孔隙率的增長放緩，而孔隙分形仍保持較為高速的降低，說明這一階段試樣的凍融損傷主要為孔隙重組；當凍融超過20次以后，全截面孔隙分維均值下降趨勢也開始變緩，僅下降原來的0.6%，是因為巖石孔隙空間已發展到足夠容納水相變時的體積膨脹，無法使孔隙進一步擴展，凍融循環作用對紅砂巖孔隙結構的影響逐漸減弱。

4 結 論

1)運用中值濾波法和直方圖均衡化手段，對煤礦砂巖的CT圖像進行預處理，極大的降低圖像噪聲，提高圖像的分辨率，為后續計盒維數的計算奠定基礎。

2)通過計算圖像分形維數，發現不同凍融次數下煤礦砂巖的截面分形維數具有相似性，說明初始損傷對巖石凍融損傷具有持續性的影響；且隨著凍融次數的增加，孔隙率逐漸增加，孔隙分形維數降低，說明巖石孔隙發育的同時其復雜度逐漸降低。

3)將煤礦砂巖細觀損傷特征與凍融荷載聯系，發現凍融損傷是由孔隙引起的有效面積損傷和分形維數表征的結構性損傷共同構成。從內部細觀角度揭示巖石的破壞機制，為指導煤礦安全開采提供一定的理論基礎。