南海钙质砂微观构造三维图像分析与表征

2023-03-27万佳怡孔德琼苏思杨彭宇朱斌

中南大学学报（自然科学版） 2023年1期

万佳怡，孔德琼，苏思杨，彭宇，朱斌,3

(1.浙江大学超重力研究中心，浙江杭州，310058；2.浙江大学建筑工程学院，浙江杭州，310058；3.浙江大学岩土工程研究所，浙江杭州，310058)

随着国家“一带一路”倡议的深入开展，在珊瑚礁区开展的海洋油气、渔业资源开发和国防建设过程规模不断增大。钙质砂作为南海岛礁的主要沉积物，主要矿物成分为碳酸钙，其沉积过程大多未经长途搬运因而保留有丰富的原生生物骨架细小孔隙[1]，表现出形状复杂、内孔隙丰富、易破碎、易胶结[2]的特征，与一般陆相、海相沉积物存在显著差异。

国内外学者针对钙质砂力学性质进行了大量试验研究。WU等[3]通过室内三轴试验定量分析了围压、颗粒级配、初始密度、排水条件和含水量等因素对颗粒破碎总量的影响。JAⅤDANIAN等[4-5]对钙质砂颗粒进行动三轴试验，分析了剪切模量、剪切刚度和阻尼比等动态力学性质。张季如等[6-7]开展室内压缩试验分析了钙质砂的压缩性能和破坏特征。WANG等[8]基于渗透试验结果，对颗粒级配、孔隙比和渗透性之间的关系进行了评估。丁选明等[9]结合振动台模型试验,对孔隙水压力、加速度、位移和动应变等动力响应进行了测试和分析。基于室内试验，相关学者在钙质砂本构模型的发展方面也取得了进展[10-11]，但由于目前在试验过程中缺少实时监测颗粒级配变化的方法，大都只关注了试样被加载至破坏或临界状态时的颗粒破碎量，而较少探讨具体应力路径下颗粒破碎的中间发展过程，因此目前所建立的颗粒破碎模型对于具体应力路径的适用性仍有待验证[11]。此外，当前对钙质砂力学特性的研究主要集中于宏观层面的室内试验与原位测试，在颗粒尺度对其运动、接触、破碎等微观构造演变行为方面的研究仍较为缺乏[12]，对颗粒级配分布、孔隙分布、颗粒形状、孔隙形状等影响渗蚀发生的内部因素研究较少[13]。

近年来，高分辨率成像设备和图像处理技术的迅速发展为岩土材料内部结构的高精度无损检测和荷载作用下颗粒尺度力学行为评估提供了新的可能。当前，高精度成像技术在颗粒材料的尺寸、方向和形态表征[14-15]以及接触关系[16]、颗粒运动[17]、局部应变[18]等颗粒尺度行为评估方面受到了广泛应用。相关学者[13,19-21]结合CT、SEM 扫描试验结果和图像分析算法对南海钙质砂颗粒三维图像形貌和内部孔隙结构分布规律进行了研究，而在钙质砂颗粒形状与尺寸参数相关性分析和颗粒间接触关系的研究等问题仍有待进一步探索。

本文作者将侧限压缩试验与CT断层扫描技术相结合，获取特定应力状态下样本中钙质砂颗粒三维图像。然后，运用自适应图像分割方法和主成分分析方法提取钙质砂颗粒尺寸、主轴方向等基本特征，进而量化钙质砂方向比、球度、凸度、棱角度等表征参数，并结合压缩过程探究受力状态下颗粒形态演化规律。同时利用膨胀算法获取颗粒间接触网络，并实现受力状态下接触对数量、方向演化规律量化和可视化的统一，以期为后续钙质砂本构模型开发与离散元数值模型验证与标定提供参考。

1 钙质砂CT扫描与图像处理

本研究选取南海钙质砂颗粒为研究对象。试验前，取适量砂样经清水洗掉盐分后烘干备用。对原始试样进行初筛后从(0.85，1.27]mm、(1.27，2.36]mm 和(2.36，3.35]mm 三个粒径级别中取等量砂样均匀混合，配置级配均匀的砂样进行试验。

1.1 “原位”侧限压缩试验

本次试验在浙江大学Nikon XTH 225/320 LC CT 扫描仪进行，并配合自主设计的微型侧限压缩装置对钙质砂颗粒进行初始表征与原位压缩试验。侧限压缩装置示意图如图1所示。装置制作材料为聚甲醛，该材料密度较低，有效减少了因高吸收率造成的X 射线穿透过程中的能量耗散，同时和其他非金属材料相比在抗压强度、耐磨损性能上更具优势。加载过程通过旋转顶部驱动螺杆，继而推动活塞向下运动控制，螺杆与装置外壳通过螺纹连接，螺距设置为0.5 mm，每旋转一圈螺杆向下移动0.5 mm。螺杆底部设计为光滑的半圆球体，从而减少因摩擦而引起的下方活塞旋转。装置顶部、底部均设有力传感器，量程为2 000 N。装置内样品的直径d=25 mm，高度H=25 mm，共有钙质砂颗粒2 000余个。

图1 CT“原位”侧限压缩装置Fig.1 In situ confining compression device

试验过程中，压缩装置固定于X 射线光源与探测器间的旋转平台上。本试验所选用的X 射线加速电压为160 keⅤ，在样品旋转180°范围内均匀采集2 500 帧投影图，其单帧曝光时间为500 ms。扫描结果分辨率设置为17.1 μm(单位体素长度)，组成颗粒的最小体素数量超过3 000，保证了砂粒真实形态表征的准确性。

试验步骤为：通过装置顶部的驱动旋钮对钙质砂颗粒压缩至目标荷载。考虑到本试验为位移控制加载，且在扫描过程中试样反力可能随时间发生松弛现象，故等待20 min 直至传感器示数基本稳定后启动CT进行拍摄，记录此目标荷载下样品三维影像。试验共设置有100，1 500 和3 000 kPa 三级荷载。

1.2 CT图像处理

为获取试样中单个颗粒的形态信息，本文采用自适应分水岭法[13]对扫描获得的样品三维CT图像进行分割。X 射线穿过样本后，受样本物质成分、密度差异影响其强度发生衰减变化，根据衰减系数生成CT扫描灰度矩阵(图2(a))；然后基于双重阈值算法确定合适的阈值对图像进行二进制化(图2(b))，从而实现实体材料与背景部分的分离。此时颗粒材料因相互接触而无法被识别为独立颗粒，通常采用基于分水岭法的图像处理技术进行分割。图像分割过程通常存在过度分割(oversegmentation)和欠分割(under-segmentation)的问题，需对二进制图像计算所得集水盆地深度进行调整[22-23]。因钙质砂形状极不规则且尺寸分布较广，设定全局控制参数调整集水盆地深度往往不适用。本文采用自适应分水岭分割法迭代进行图像分割，在每一次迭代过程中，针对新分割出的区域计算相应的局部调整参数，且在迭代过程中保留欠分割颗粒至下一次迭代继续处理，从而实现图像的精细分割。图2(c)所示为最终分割图像的二维切片，其中不同颜色区域即对应不同颗粒。可以看到该方法在分割钙质砂这一类不规则颗粒土样的图像中具有较好的效果。将扫描结果与石英砂扫描结果对比[24]，可以看出钙质砂颗粒形状复杂，内部保留有丰富的原生生物骨架细小孔隙。以上过程均在MATLAB 中进行，其他具体细节参考文献[13]。

图2 CT图像分割(仅展示二维切片)Fig.2 Illustration of image segmentation(showing a 2D slice)

2 钙质砂三维形态表征

本文基于图像分割结果对装置内全部钙质砂颗粒进行逐个分析，进一步提取钙质砂颗粒形心坐标、体积和三维尺寸等基本信息，进而计算出钙质砂颗粒的主轴方向、平坦度、伸长度、球度和凸度等典型形态参数，实现钙质砂颗粒三维形态的定量表征。

2.1 颗粒形态参数确定

CT 扫描重构可获得代表钙质砂颗粒的一系列相同标号三维顶点坐标。由于钙质砂颗粒具有形状复杂、不规则的特点，采用主成分分析(principal component analysis,PCA)方法对样本中每个单独颗粒的顶点坐标进行正交变换，从而将原始坐标转化为互不相关的新变量长轴、中轴和短轴，进而实现不规则颗粒尺寸、方向等基本信息的标准化测量。根据PCA 计算结果，将其从小到大依次定义为颗粒短轴模量c、中轴模量b和长轴模量a。参考KONG等[14]对钙质砂颗粒轴长和筛分试验对比结果，采用长轴模量a作为钙质砂颗粒等效粒径D。

基于长轴模量a、中轴模量b和短轴模量c计算颗粒形状比参数，从而对颗粒形状进行度量。形状比参数包括伸长度E和平坦度F，通常采用E=b/a和F=c/b分别进行计算[25]。图3所示为2个代表性颗粒模型及其伸长度、平坦度计算结果。可以看出伸长度E、平坦度F接近于1 时，颗粒各方向轴长模量近似相同。采用球度对颗粒球状性进行定义，通常采用(其中，S0为颗粒球度，V为颗粒体积，SA为颗粒等效球表面积)进行计算[24]，但当颗粒表面较为粗糙时使用该方法会使得计算结果产生较大误差。因此，由于钙质砂具有多棱角、形状复杂的特征，本文定义球度指标S=Vfill/VS(其中，Vfill为填充内孔隙后的颗粒体积，VS表示颗粒外接球体积)。图4所示为2个代表性颗粒模型及其球度计算结果。

图3 典型颗粒平坦度与伸长度Fig.3 Typical grains with different flatness and elongation values

图4 典型颗粒球度Fig.4 Typical grains with different sphericity values

采用凸度对颗粒表面凹凸程度进行定义，采用公式C=V/VCH(其中，C为颗粒凸度，VCH为包裹住颗粒所需的最小凸包体积)计算。图5所示为2个代表性颗粒模型及其凸度计算结果。从图5可以看出：凸度越大，颗粒表面越平整，凹凸程度越小。

图5 典型颗粒凸度Fig.5 Typical grains with different convexity values

为了对颗粒整体棱角特征进行评估，通常采用棱角度对颗粒表面形态进行量化。首先需对颗粒表面局部棱角特征进行定义，通过MATLAB 内置函数vertexNormal 和vertexAttachments 计算某一顶点与其相邻顶点矢量，获取该顶点最大曲率k1和最小曲率k2。之后，本文通过颗粒表面顶点平均曲率来判断该顶点是否为一个棱角，若该顶点处平均曲率km,j对应的平均曲率半径1/km,j小于颗粒内切球半径Rin，则判定顶点为一棱角。图6所示为2个典型颗粒表面平均曲率计算情况。从图6可以看出：棱角处对应的平均曲率较大，平缓处平均曲率较小，表明平均曲率能较好表征钙质砂颗粒局部棱角特征。在此基础上，定义棱角度A计算方法如下：

式中：kin为颗粒内接球半径的曲率；km,j为该颗粒第j个面元的平均曲率；Aj为第j个面元的面积。图6中典型颗粒棱角度计算结果显示，表面棱角丰富的钙质砂颗粒棱角度较大，而表面棱角较为平缓的颗粒棱角度较小，表明棱角度能够较好地表征钙质砂颗粒整体棱角情况。

图6 典型颗粒棱角度Fig.6 Typical grains with different angularity values

2.2 结果与讨论

图7所示为试验过程中3次扫描获得的钙质砂颗粒等效粒径的累积分布曲线演化情况。从图7可以看出：在加载过程中，钙质砂颗粒等效粒径中位数D50仅发生微弱减小，即在压力P为3 MPa 条件下，级配均匀的钙质砂样本并未表现出明显的颗粒破碎现象，与一般情况下大颗粒钙质砂在该压力下出现较大的颗粒破碎不符，这可能是良好的级配显著降低了钙质砂的破碎性所导致[26]。

图7 压缩过程颗粒等效粒径演化Fig.7 Grain size evolution during compression process

FONSECA 等[16]建议采用颗粒主轴a方向作为具体颗粒的方向向量。将每个颗粒形心空间坐标作为向量起点，主轴方向作为颗粒方向。为了探究压缩过程中钙质砂颗粒方向演化关系，本文分别求出颗粒主轴方向与X轴、Z轴的夹角，并绘制出压缩过程中颗粒方向的演化情况，如图8所示。从图8可以看出：在加载过程中，颗粒与X轴夹角逐渐减小，与Z轴夹角趋于垂直，即加载过程中样品中颗粒方向趋于向垂直于受力方向的平面演化。

图8 压缩过程颗粒方向演化Fig.8 Grain orientation evolution in compression process

图9所示为南海钙质砂与英国近海石英砂[24]三维形态参数的累积分布曲线，包括凸度和伸长度。从图9可以看出：南海钙质砂的凸度小于英国近海石英砂的凸度，而伸长度更大，这定量地说明了相较于石英砂颗粒，南海钙质砂外形更加复杂，棱角更加明显。此外，可以注意到相较于石英砂颗粒，南海钙质砂颗粒的所有三维形态参数的累积分布曲线的分布范围更大。上述研究结果表明，南海钙质砂颗粒受风化和剥蚀作用影响，形态更加复杂，棱角显著，形态特征变异性较大；而石英砂颗粒在地质搬运过程中趋于规则演化，形态特征变异性较小[27]。

图9 钙质砂与石英砂颗粒形态参数分布对比[24]Fig.9 Cumulative distribution curves of various morphological parameters of calcareous sand and quartz sand[24]

图10所示为南海钙质砂不同形态参数与颗粒尺寸的分布关系，包括球度和棱角度。从图10 可以看出：随着颗粒尺寸的减小，颗粒球度不断增加，而棱角度出现明显减小。该结果表明，细颗粒砂由于受地质搬运过程历时更长，颗粒形态趋于规则演化，棱角更加平缓。同时，受顶端压缩影响，钙质砂颗粒球度分布范围略有所提高，而棱角度分布范围则出现整体下移，表明压缩过程使钙质砂颗粒形态发生改变，且变化方向趋于规则。此外，分析结果表明伸长度等其他形态参数与颗粒尺寸间不存在明显的相关关系，因此本文仅给出球度、棱角度与颗粒尺寸分布关系。

图10 颗粒形态参数与粒径关系Fig.10 Relationship between size and shape parameters

3 颗粒间接触演变

颗粒间接触关系是影响砂土力学性质的重要因素之一。本研究采用基于MATLAB 膨胀算法的颗粒接触关系捕捉方法，将目标颗粒进行膨胀运算后与图像分割后的周边颗粒叠加，标记出现接触的部分从而确定配位数(coordination number)。此外，本方法仅对连接在同一分水线上的颗粒执行接触捕捉算法，有效避免了将非接触颗粒识别为接触对的错误发生，且与文献[16]采用的算法相比极大地提高了运算效率。图11所示为4 组典型颗粒接触对识别结果，图中蓝色箭头表示接触对方向向量，其中向量起点表示颗粒间接触点，向量方向指向两颗粒接触面法线方向。

图11 典型钙质砂配位数和接触关系Fig.11 Typical grains with different coordination number and contact network

图12(a)所示为加载前颗粒配位数与粒径分布情况。从图12(a)可见：颗粒配位数与颗粒尺寸呈明显的正相关关系，且颗粒最大可能配位数与等效粒径呈明显的线性关系。为研究压缩过程中配位数与粒径分布关系演化情况，分别计算特定粒径范围内配位数均值，得到平均粒径与平均配位数的关系(图12(b))。由于粒径在6～8 mm范围内的颗粒数量极少，为减小计算误差仅对粒径在1～6 mm范围颗粒进行分析。可见，伴随压缩过程土体密实程度增加，大颗粒钙质砂配位数出现明显增加。

图12 颗粒配位数与粒径关系演化Fig.12 Evolution of relationship between grain coordination number and size

将样品中所有颗粒的接触向量识别结果绘制于同一坐标系下，得到整个样品接触网络如图13(a)所示。本文将接触面法向量方向定义为接触方向，计算样品中各接触方向与X轴夹角，并绘制出压缩过程接触方向的演化情况，如图13(b)所示。可见：样品中颗粒间接触方向明显垂直于加载平面，且随着加载过程进行，接触方向继续向垂直于加载平面方向演化，这与文献[16，28]中认为接触方向趋向于平行于主应力方向从而形成稳定力链的结论一致。