基于三维重构钙质砂的结构特征及相关性分析

王志兵，孙广，刘金明，邹永胜，于江涛

(1.桂林理工大学土木与建筑工程学院，桂林 541004；2.广西岩土力学与工程重点实验室，桂林 541004；3.广州协安建设工程有限公司，广州 510075)

南海是中国南部边防重要领地和对外贸易的重要通道，而且富含大量海洋资源[1]。在“一带一路”倡议、“海洋经济与海洋强国战略”大背景下，在南海地区建立了大量的岛礁工程，钙质砂广泛存在于南海地区，因此常被用来作为岛礁的填筑材料。钙质砂一般由海洋生物遗骸沉积形成，形成的钙质砂孔隙较多且颗粒形态复杂多变，这种多孔以及形貌复杂的特性极大地影响了钙质砂的力学性质[2]，因此研究钙质砂的结构特征及其相关性是有必要的。

已有大量学者对钙质砂结构特征进行了研究。吕海波等[3]通过扫描电镜试验描述了孔隙形状，运用压汞试验研究了钙质土的孔径分布。朱长岐等[4]应用分秒激光切割技术将钙质砂进行切片处理，使用图像处理软件分析了钙质砂横断面孔隙，研究得出了随着钙质砂颗粒尺寸的增大，其断面孔隙度也会增大。将明镜等[5]使用MATLAB开发的图像处理程序分析了大量钙质砂扫描电子显微镜(scanning electron microscope，SEM)图像，定量研究了不同粒径和形状钙质砂的孔隙特征。周博等[6]应用计算机断层扫描(computed tomography，CT)技术重构了钙质砂三维结构，并计算了钙质砂颗粒的孔隙率、孔径分布、分形维数。曹培等[7]使用压汞试验研究钙质砂的连通孔隙的孔径分布，并运用CT扫描试验提取钙质砂孔隙三维结构分析了其封闭内孔隙度，发现钙质砂的封闭孔隙度远低于连通孔隙度。除了以上对于钙质砂孔隙特征的研究以外，一些研究人员也进行了钙质砂形貌特征的研究。陈海洋等[8]通过二维钙质砂图像，通过颗粒的长宽比和分维数较好地描述了钙质砂的二维形貌特征。王步雪岩等[9]通过颗粒形貌观测系统从不同投影面研究了大量钙质砂的二维结构特征，可将钙质砂从形貌上快速进行分类。赵书辉[10]、张斌等[11]分别通过三维和二维方法描述了钙质砂的形貌特征，并探究了形状对于钙质砂力学性质的影响。虽然已有大量研究人员进行了钙质砂的结构特征研究，但大多研究是基于二维平面，而且不同方法对于参数的提取是有局限的。

为了较全面地研究钙质砂的结构特征，现基于高精度X射线CT扫描技术建立了钙质砂颗粒与孔隙的三维结构。通过Avizo软件提取分析钙质砂的孔隙度、封闭孔隙度、颗粒形状系数、孔隙分形维数等结构特征，并对这些参数进行相关性分析，以期为预测钙质砂的结构特征提供一定的理论依据。

1 材料及图像采集

1.1 材料

试验所用钙质砂取自南海某岛礁，试样为珊瑚屑沉积物且互不胶结。共选取了7个钙质砂作为研究对象，粒径为10～20 mm，其形状分为块状、片状、条状三种，钙质砂样如图1所示，其中1号、2号、3号钙质砂形状为块状，4号、5号钙质砂形状为片状，6号、7号钙质砂形状为条状。由于钙质砂的易破碎性，取样后放入防震泡棉盒中，在运输过程中尽量减少震动和挤压。送至实验室需用蒸馏水洗尽钙质砂的盐分及附着杂质，放入105～110 ℃烘箱中烘烤8 h以上，然后存于自封袋内备用。

图1 钙质砂样

1.2 图像采集

1.2.1 CT扫描原理

X射线是一种具有穿透性的电磁波，当其穿过物质时，由于物质的组成成分不同，对于X射线的衰减程度也不同[12]。基于X射线的性质，人们开发了可用于医学、工业、农业等领域的CT设备。高分辨率工业CT设备一般由X射线源、载物台、探测器三部分组成。采集图像时，X射线源及探测器固定不动，将样品放置在载物台上，然后精确控制载物台旋转角度，待累计旋转角度达到360°时完成数据采集，再通过计算机一系列处理生成最后的CT图像，计算机断层扫描系统示意图如图2所示。

图2 计算机断层扫描系统(CT)示意图

1.2.2 试验装置

采用桂林理工大学有色金属及材料加工新技术教育部重点实验室的Xradia 510 Versa高分辨率三维X射线显微镜进行钙质砂的CT图像采集，仪器内部构造如图3(a)所示。扫描得到930张1 000×1 024像素CT切片，由于钙质砂颗粒大小不一，扫描精度为10～20 μm，颗粒越大单个像素点大小所对应的值越大。

图3 仪器内部构造及试样

2 三维重构及参数定义

2.1 图像处理

2.1.1 滤波处理

通过CT扫描得到的灰度图像会存在系统噪声，这种噪声反应在灰度图像上会呈现出无数的孤立像素点，如图4(a)所示。这些噪声无疑降低了图像的质量，在后续的阈值分割中会造成干扰(如出现大量微像素点)，因此需要对原始CT图像进行滤波处理来消除噪声的影响。

常用的滤波方法有：中值滤波、高斯滤波、非局部均值滤波。本文研究采用了三种方法对钙质砂原始CT图像进行滤波处理，如图4所示，通过对比这三种滤波的效果，发现中值滤波和高斯滤波对图像轮廓的平滑及对噪声的过滤效果比非局部均值滤波的效果差，而且非局部均值滤波较完整地保留了原始灰度图像的边缘与细节特征，这与耿冲等[13]的结论一致。因而选用非局部均值滤波进行滤波处理。

图4 滤波处理

2.1.2 阈值分割

在钙质砂的CT图像中，存在着钙质砂颗粒、钙质砂内孔隙以及颗粒外的“空气”三种组成成分，如何有效的区分它们是决定三维重构精度的重要因素。对于钙质砂颗粒的重构，只需确定钙质砂颗粒与孔隙的阈值，而最大类间方差法(Otsu法)一直被认为是最优的方法[14]，许多学者应用此方法进行阈值分割取得了不错的效果，本文研究也采用Otsu法确定阈值，阈值分割后得到钙质砂颗粒的二值图像，如图5(a)所示。然而通过常规的阈值分割方法，在相同的灰度范围内难以区分内孔隙和颗粒外的“空气”，而采用形态学顶帽算法(Top-hat)可以较好地提取出由颗粒包裹形成的孔隙，钙质砂孔隙的二值图像如图5(b)所示。

阈值分割完成后，通过连续切片的堆叠，如图5(c)所示，应用Avizo软件中的“Volume Rendering”即可展示钙质砂颗粒与孔隙的三维重构图像，如图6和图7所示。

图5 二值图像

图6 钙质砂颗粒三维图像

图7 钙质砂孔隙三维图像

2.2 参数定义

2.2.1 孔隙度

孔隙度是指孔隙体积与总体积的比值，即

(1)

式(1)中：n为孔隙度；Vpore为所有孔隙的体积；Vtotal为所有孔隙和颗粒的体积和。

2.2.2 孔隙曲折度

孔隙曲折度是指孔隙的路径长度与切片高度的比值。该模块首先计算图像每个平面的质心，然后计算通过质心的路径长度，最后除以切片高度。计算原理如图8所示，计算公式为

(2)

式(2)中：T为曲折度；dj为相邻切片间孔隙质心路径长；i为切片数；H为切片高。

2.2.3 颗粒形状系数

对于颗粒的三维形状评价一般采用球度、扁度、棱角度等，但考虑到钙质砂的内部多孔、表面多沟壑的特性，在此采用吴野等[15]改进Wadell形状系数φ来作为钙质砂形状系数。从表达式可以看出，R的范围为0～1，R越大，表示颗粒投影面越接近圆形，表面越光滑。φ的范围为0～1，φ越大，表示颗粒形状越接近圆球形，轮廓越规则。计算公式为

Z(i)、Z(i-1)、…、Z(0)分别为第i+1、i、…、1张切片图8 曲折度模块计算原理

(3)

式(3)中：φ为形状系数；As为与颗粒具有相同体积的球体的表面积；Ap为颗粒表面积；R为颗粒最大投影面二维的圆度，计算公式为

(4)

式(4)中：R为圆度；S为颗粒与内孔隙之和等面积的圆周长；P为钙质砂颗粒外轮廓周长。

2.2.4 孔隙三维分形维数

分形维数是度量研究对象的自相似性和定量性质的指标，也是测量和比较表面粗糙度的有效指标[16-17]。三维的分形维数在2～3，对于标准几何曲面(立方体、平面、椭球体等)分形维数为2，表面越粗糙，对应分形维数值越大。

计盒维数法广泛应用于材料分形维数的研究，其含义为分形体对空间的覆盖极限。具体方法先用边长为δ的“盒”把需要分形的三维结构覆盖起来，三维结构中包含了实体与孔隙，所以一部分“盒”覆盖了实体，另一部分“盒”覆盖了孔隙。然后，把不是空的盒子数统计出来，记为N(δ)，再缩小“盒”的尺寸δ，随之N(δ)的数目将增多，计盒维数的计算公式为

(5)

式(5)中：Dbox为分形维数；δ为正方体“盒”的边长；N(δ)为非空“盒”的数量。

3 结果与分析

3.1 参数结果

基本参数如孔隙总体积、颗粒体积、孔隙度、颗粒表面积可以通过Avizo软件中的“Label Analysis”分析计算得到。而钙质砂孔隙包含了连通孔隙和封闭(孤立)孔隙，需要应用“Axis Connectivity”进行连通孔隙的提取，再通过减法运算即可得到封闭孔隙。提取到连通孔隙后，通过“Centroid Path Tortuosity”可计连通孔隙的曲折度。孔隙的分形维数通过以上计算方法得到一系列的双对数散点图，通过最小二乘法的线性拟合得到斜率，即分形维数，如图9所示。钙质砂的结构特征参数如表1所示。

表1 钙质砂结构特征参数

图9 计盒维数拟合结果

依据颗粒外形形态，初步将钙质砂颗粒分为两类，第一类钙质砂颗粒为1号、2号、4号和5号，表面孔隙较多，沟壑纵横起伏。第二类钙质砂颗粒为3号、6号和7号，表面孔隙较少，形状也较圆润。这种差异性可能是由于不同生物环境造成的。而且第一类颗粒的形状系数φ值均小于第二类颗粒的形状系数φ，因此用形状系数φ来描述钙质砂颗粒形状是较合理的。

由CT重构得到的钙质砂连通孔隙度在27%～38%，孔隙曲折度在1.7以上，封闭孔隙度一般在1%以内，再次验证了钙质砂是一种内部多孔且孔隙结构较为复杂的岩土工程材料。研究人员常用压汞试验、氮气吸附试验和CT扫描试验研究材料的孔隙特征，但每种方法都有自身的局限性。例如，压汞试验是一种破坏性试验，在侵入汞的过程中会破坏材料，从而会高估材料的孔隙体积，存在“墨水瓶”效应，而且对于封闭孔隙的检测无能为力。氮气吸附试验不适用与大孔的检测，与压汞试验相同，无法检测封闭孔隙。CT扫描试验作为一种无损的试验方法，在试验后还可以进行其他相关试验，而且对颗粒或孔隙可以进行定量分析，但这种方法受限于扫描的精度，无法提取到小于该精度的颗粒或孔隙。

3.2 相关性分析

两个变量的关联程度被称为相关性，相关性按关联的密切程度一般可分为完全相关、不完全相关、无关。相关性系数r取值在-1～1。当r为正值时，说明为正相关，当r为负值时，说明为负相关。|r|越接近于1时，说明其相关性越好，|r|=1时则为完全相关，|r|=0时为无关[18]。

为了研究钙质砂结构特征之间的关系，对表1中的结构特征参数做相关性分析，如表2所示。参照潘宗源[19]对相关程度的分类方法，可将相关性分为：极强相关(0.8≤|r|<1)、强相关(0.6≤|r|<0.8)、中度相关(0.4≤|r|<0.6)、弱相关(0.2≤|r|<0.4)、极弱相关(0<|r|<0.2)。

表2 钙质砂结构特征的相关系数矩阵

从表2可以看出，除连通孔隙曲折度、颗粒形状系数外，其余6项参数都与颗粒体积有中度相关以上的正相关性。孔隙度作为钙质砂重要的且难以大量测量的重要指标，与钙质砂颗粒体积(颗粒大小)存在正相关的关系，绘制钙质砂颗粒体积与孔隙度的散点图，如图10所示。从图10中可以发现3号、6号点是异常点，这是由于3号、6号钙质砂颗粒形状系数均大于0.3，属于上文所划分的第二类钙质砂，表面孔隙较少，形状也较圆润。除去异常点的影响，可以发现随着钙质砂孔隙体积的增大，钙质砂的孔隙度也在增大，增大到一定程度时增长速率减缓。孔隙分形维数与颗粒形状系数呈极强负相关，说明颗粒形状越粗糙，对应的孔隙形状也越复杂。与形状系数类似，第一类钙质砂孔隙分形维数均大于第二类钙质砂孔隙分形维数。

图10 不同体积钙质砂的孔隙度

4 结论

通过CT扫描技术重构了钙质砂的颗粒与孔隙的三维结构，通过Avizo软件提取了颗粒体积、孔隙体积、孔隙度、曲折度、颗粒形状系数等结构特征参数，计算了钙质砂颗粒的分形维数。最后建立了钙质砂结构特征的相关性系数矩阵，对部分指标进行了相关性分析，主要得出以下结论。

(1)三维重构钙质砂的孔隙度在27%～38%，孔隙曲折度在1.7以上，孔隙连通率一般在98%以上，通过三维重构的方法从三维空间上定量的验证了钙质砂是一种内部多孔且孔隙结构较复杂材料。

(2)钙质砂可以通过形状系数来进行分类，形状系数越小，颗粒表面越粗糙。按形状系数0.3为界，分为两种类型钙质砂，第二类钙质砂随颗粒体积的增大，其孔隙度也在增大，增大到一定程度时增长速率减缓。

(3)颗粒体积与除连通孔隙曲折度、颗粒形状系数外的其余6项参数有中度相关以上的正相关性。孔隙分形维数与颗粒形状系数呈极强负相关，说明颗粒形状越粗糙，对应的孔隙形状也越复杂。