冻融岩石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)将煤矿砂岩细观损伤特征与冻融荷载联系，发现冻融损伤是由孔隙引起的有效面积损伤和分形维数表征的结构性损伤共同构成。从内部细观角度揭示岩石的破坏机制，为指导煤矿安全开采提供一定的理论基础。