碎石土大孔隙结构定量化表征

吕澄浒，樊秀峰,2，吴振祥,2，王晓摇

(1.福州大学 紫金地质与矿业学院， 福建 福州 350108; 2.地质工程福建省高校工程研究中心， 福建 福州 350108)

碎石土是一种土石混合体，颗粒大小混杂，结构多元，土体内分布着许多大孔隙，颗粒粗大的大孔隙域与颗粒细小的基质域并存，不同于其他均质边坡相对均匀的渗流特性，此类边坡易形成局部集中、流速快的渗流网[1]，在地下水和降雨条件下，水流绕过大部分土壤基质域，沿着大孔隙通道快速进入深层[2]，形成连通大孔隙优先渗流通道，引发碎石土边坡失稳等自然灾害。目前已有研究表明碎石土大孔隙空间分布特征对优先流的产生密切相关。因此，展开碎石土大孔隙结构特征定量化研究，对预防此类边坡灾害发生十分必要。

目前，国内外学者对大孔隙界定方法主要有当量孔径法、毛细势界定法、导水率法及功能界定法等[3]，其中当量孔径法是最直接的方法。国外方面，Beven等[4]曾将土壤中大孔隙当量孔径界定为0.03 mm～3.00 mm；Carey等[5]将大于1 mm的孔隙定义为大孔隙。国内方面，吴华山等[6]将孔径大于0.3 mm孔隙看做大孔隙；王伟等[7]在研究林地得到大孔隙半径大多在0.3 mm～3.0 mm之间。

已往对大孔隙结构研究的方法主要有：水分渗透曲线法[8]、张力入渗仪法[9]、切片法[10]、染色示踪法[11]。在众多学者们既往研究试验过程中，难以避免会对土样造成不同程度的扰动，影响试验结果的准确性。而CT扫描法[12-13]是一种非破坏性探测技术，它可以在不影响土体结构的前提下，快速精确获取土壤的大孔隙结构，且CT技术对扫描碎石土稍大尺寸土样有一定的优势。随着高精度CT扫描仪的应用，极大地提高了扫描的准确度，从而实现对小于1 mm的大孔隙结构进行定性分析。此外，已有对碎石土大孔隙结构的研究只是局限于定性，大孔隙空间分布规律如何定量描述方面还没有涉及，其次对大孔隙参数定量分析不够深入，仍仅是大孔隙度、数量和直径等参数的初步统计分析，研究成果未成体系。

综上所述，本文在已有丰富的研究成果基础上，采用图像0.32 mm分辨率的高精度CT扫描机，故将识别到等效直径大于0.32 mm的孔隙定义为大孔隙。通过对CT图像进行二值化和大孔隙参数提取，以细观角度定量分析大孔隙特征和分布规律，从而为水分沿其快速下渗研究提供可靠理论依据，这对于强降雨型碎石土边坡失稳的研究、揭示发生机理具有重要的意义。

1 取样与研究方法

1.1 原状土取样

选取福建安溪某地碎石土边坡作为取样地点，在该试验区选择两个距离5 m左右的区域分别挖取相同深度的土样作为A和B，图1为原状土取样图。为保护原状土的完整孔隙结构，本文采取从四周竖直向内开挖的方式进行取样。具体取样过程如下：用量尺划定一个大约半径为50 cm的圆形区域，由周边向圆心方向依次进行垂直人工开挖，开挖过程中确保中间半径20 cm的土体不受扰动，再将土柱修剪成尺寸半径为20 cm，高为30 cm的圆柱体。最后，将在现场烧融的石蜡浇注在土柱周围，对其进行蜡封，再用保鲜薄膜将蜡封的土柱进行密封，小心装入特制好的土柱模型中，用法兰盘固定好，示意图见图2。

(a) 土样A (b) 土样B

图2 土柱模型装置示意图

1.2 土样的基本物理力学性质

试验所取土样为崩坡积的典型碎石土，其天然状态下的物理指标及颗粒级配曲线见图3和表1。由图3可知现场取回土样的不均匀系数Cu分别为26.47和25.01，说明土样中的颗粒级配悬殊、差异大，颗粒级配曲线的曲率系数Cc为0.60和0.52，级配较差。

表1 天然状态碎石土的基本物理指标

图3 天然状态碎石土的颗粒级配曲线

1.3 CT扫描过程

本文选用Revolution 256排512层超高端螺旋CT扫描机是美国GE公司耗时7年研制的尖端产品，如图4所示。将土柱试样居中放置扫描架，扫描仪器参数调试完成后，将相邻两CT切片的实际距离设置为0.625 mm，对整个土柱“自上而下”进行360度全方位扫描，每个土样扫描共得到480个切片图像，得到CT扫描值范围-1 024 HU～3 071 HU。

图4 CT扫描仪

1.4 CT图像二值化处理

本文将在已有二值化处理研究的基础上对操作步骤进行优化：先将切片的CT值分割为-1 024 HU～-500 HU，-500 HU～1 300 HU，1 300 HU～3 071 HU三个等级范围，并分别赋予其不同的颜色，绿色部分代表碎石，黄色部分代表土壤基质和根系，黑色部分代表孔隙，见图5(a)，通过图像处理软件对赋予颜色的CT图像进行统一裁剪，裁去渗流区域以外的CT图像，见图5(b)，再通过Python软件中执行灰度图像处理命令自动生成最佳阈值的灰度图，继续执行二值化命令，输出二值化图像，图像中黑色为土体中的大孔隙，白色为碎石、土壤基质等结构，见图5(c)。CT扫描成像过程中，由于土体结构对X射线存在吸收作用，使得CT图像出现不可避免的噪点，最后通过图像处理软件Image-Pro-Plus 6.0对二值化的图片进行滤镜降噪处理，得到完整的二值化图像见图5(d)。

图5 CT图像二值化过程

1.5 大孔隙参数的提取

通过Image J软件灰度识别功能，精确提取二值化处理后CT图像的大孔隙形状特征，再运用颗粒分析功能对大孔隙面积、周长、尺寸等形状参数进行定量化分析；其中大孔隙尺寸包括费雷特直径(Feret)、最大费雷特直径、最小费雷特直径、等效面积圆直径等。由于以上几何参数是通过软件直接从图像分析获取的，因此被称之为基本几何尺寸参数[14-16]，相关参数示意图如图6所示，已有研究表明这几个参数可以客观表征大孔隙的空间分布规律[14-18]，本文将这种方法迁移到碎石土大孔隙研究中，进一步分析了大孔隙的形状参数，如扁平度、成圆率、整体轮廓系数和圆度。

图6 大孔隙几何参数示意图

扁平度(Aspect Ratio，AR)表征大孔隙的狭长程度，其定义为最大Feret直径与最小Feret直径的比值，计算公式如下：

(1)

式中：AR为扁平度；L为该大孔隙单元的最大Feret直径，mm；B为该大孔隙单元的最小Feret直径，mm。

成圆率(Form Factor，FF)表征大孔隙形状接近于圆的程度，即大孔隙形状越接近于圆形，其数值越趋近于1，反之，表示大孔隙形状越不规则，计算公式如下：

(2)

式中：FF为成圆率；A为大孔隙面积，mm2；P为大孔隙周长，mm。

整体轮廓系数(Convexity)表征大孔隙棱角突出程度和轮廓粗糙变化程度，定义为等效面积圆的周长与孔隙周长的比值，计算公式如下：

(3)

式中：Con为整体轮廓系数；D为等效直径，mm；P为大孔隙周长，mm。

2 结果分析

2.1 大孔隙分布特征

2.1.1 大孔隙度随深度变化趋势分析

大孔隙度是CT图像中各大孔隙面积占CT图像总面积的百分数，其确定方法是计算二值化处理后的黑色像素区域面积和完整图像像素面积，二者之比可得到该CT切片的大孔隙度，其中黑色像素代表大孔隙。由于CT扫描过程时将切片间距设置为0.625 mm，共得到480张切片图像。若统一处理分析，工作量十分冗杂，为了提高工作效率，每隔10张提取一张切片图像进行统计分析，即层间距为6.25 mm。对垂直深度方向下大孔隙度的变化趋势进行统计学分析，得到大孔隙度随深度方向的变化趋势，见图7。从图中可以看出，两个试样在不同深度处大孔隙度变化范围介于11.9%～33.8%之间，大孔隙度随着深度先增加，到深度180 mm后开始逐渐减小的趋势，均是位于较浅的土样大孔隙度大于位于较深位置的大孔隙度，并且两个土样在同一深度处的大孔隙度也是较为接近的。

2.1.2 不同孔径级别的大孔隙数量分析

为研究不同孔径大孔隙在垂直深度方向上的分布情况，本文从两个土样中各选择10张CT切片图像，分别进行处理、计算和统计。通过Image J软件对大孔隙数量按孔隙尺寸将其分为不同的孔隙等级进行统计。计算结果显示土样A和土样B的大孔隙等效直径上限值分别为23.21 mm和25.67 mm。在已有研究基础上，本文将大孔隙按等效直径大小划分为8个级别，即：0.32 mm～2 mm(1级)，2 mm～5 mm(2级)，5 mm～8 mm(3级)，8 mm～11 mm(4级)，11 mm～14 mm(5级)，14 mm～17 mm(6级)，17 mm～20 mm(7级)和大于20 mm(8级)。

图7 大孔隙随深度方向上变化趋势图

图8为不同孔径级别大孔隙数量占总数比值及累计占比曲线图。

图8 不同孔径级别大孔隙数量分布及累计曲线占比图

从图8中可以看出大孔隙数量占比主要集中在第1—第3孔径级别，且两个土样在第1—第3孔径级别中大孔隙数量占该总数量比重分别为73.11%和75.65%，可以得出碎石土中有75%左右以上大孔隙直径在0.32 mm～8.00 mm范围内。随着深度的增加，第4—第6孔径级别的大孔隙数量比重明显高于表层，达到一定深度后，其大孔隙数量比重也越来越少。与以往研究结果相比，植物根系腐烂和干湿循环等因素引起的孔径大都是小于1 mm，动物活动引起的孔径在2 mm～4 mm范围内[19]，由此可见，表层碎石土的大孔隙主要是由于植物根系腐烂、干湿循环等因素造成的，而颗粒粒径相差悬殊是造成孔隙大于8 mm的主要原因。

2.2 大孔隙形状特征分析

对两个土样的各10张CT切片图像二值化处理并进行计数统计后，对各CT切片中各大孔隙形状参数的平均值进行了计算：圆度、成圆率、整体轮廓系数和扁平度，并分析各参数在垂直深度上的分布趋势。

从图9可以看出圆度、成圆率、整体轮廓系数在深度0 mm～90 mm范围内随深度增加，其平均值无明显变化；当深度超过90 mm左右时，圆度和成圆率均随着深度增加有轻微降低趋势，到深度大于210 mm左右后，圆度和成圆率逐渐恢复为浅层深度的变化规律。同一试验区距离相近的两个土样在同一深度上的形状参数差别不大，沿垂向深度的变化趋势也基本一致。

图9 各大孔隙形状参数平均值随深度变化趋势图

2.2.1 扁平度分析

图10为大孔隙扁平度随深度变化趋势图。从图可知，两个土样的大孔隙扁平度在深度0 mm～90 mm范围内变化不大，但超过90 mm左右后发生较为明显波动，呈现增加的趋势；在深度180 mm后，扁平度又逐渐减小趋于稳定。认真对比90 mm～180 mm范围的CT图像，如图11所示，可见该范围内切片图像的大孔隙形状呈狭长且不规则，故扁平度大于其它区域。

图10 大孔隙扁平度随深度变化趋势图

图11 A#z=120 mm和B#z=150 mm的二值化图像

2.2.2 圆度分析

为探究土样浅层与深层大孔隙形态的差异性，故选择深度为60 mm和210 mm处切片图像，分别对两个土样进行大孔隙圆度统计分析。结果表明，圆度均介于0.07～1.00范围内。为此，以0.1为一个级别，将圆度划分10个等级，对每个圆度等级的大孔隙数量进行统计分析。图12为圆度分布图，在不考虑0.9～1.0之间形态分布的情况下，其余大孔隙均服从正态分布，并且期望值均为0.6左右。对比浅、深层的切片图像可以看出，深层切片图像的大孔隙在0.2～0.4级别的频数明显大于浅层切片，这表明深层土体的大孔隙形态比浅层的更复杂，大孔隙形状呈狭长且不规则。

2.2.3 成圆率分析

在得知土样深度对大孔隙成圆率无显著影响后，为探究不同深度处碎石土大孔隙成圆率与孔隙直径的关系，在上文的基础上统计不同孔径级别(与2.1.2节大孔隙直径级别划分相同)的平均成圆率，统计结果见图13。从图中可知，大孔隙直径对成圆率有显著影响，且随着大孔隙直径的增大，其成圆率呈线性降低，即孔径越大，其成圆率越小，大孔隙形状愈不规则。为此在图13基础上，利用一次函数y=ax+b进行拟合，进一步探究大孔隙直径与成圆率的关系。根据拟合结果可知，两个土样的大孔隙直径与成圆率之间的平均相关性均大于0.9，各土样不同深度的大孔隙成圆率与孔径级别的拟合函数较为相似，表明碎石土的大孔隙成圆率与孔径级别之间存在一定的线性关系。

图12 圆度分布图

图13 不同孔径级别大孔隙成圆率变化趋势图

2.2.4 整体轮廓系数分析

图14为不同孔径级别孔隙整体轮廓系数变化趋势图。对两个土样各10张CT切片图像上的大孔隙整体轮廓系数按不同孔径级别(与2.1.2节大孔隙直径级别划分相同)进行统计计算，并求出各孔径级别大孔隙整体轮廓系数的平均值。由图可知，两个土样均整体表现出大孔隙整体轮廓系数随着孔径级别增大而逐渐减小的趋势，这表明大孔隙的等效直径越大，整体轮廓系数越小，其轮廓的凹凸程度越明显。另外，在总体变化趋势表现平滑情况下，而深度240 mm和270 mm处切片的大孔隙整体轮廓系数受孔径级别影响的变化趋势波动程度较大，主要是由于这两处的大孔隙数量相对来说较少，导致其统计分析离散性也较大。

图14 不同孔径级别孔隙整体轮廓系数变化趋势图

3 结 论

(1) 碎石土大孔隙具有分布离散性大、结构差异性明显、形状不规则等特点。

(2) 碎石土在不同垂向深度处大孔隙有明显差异，大孔隙度变化范围介于11.9%～33.8%之间，大孔隙度随着深度先增加后减小的趋势，位于浅层的土样大孔隙度大于深层位置的土样大孔隙度。

(3) 碎石土样中约有75%以上大孔隙直径在0.32 mm～8.00 mm范围内，孔径越大，其孔隙数量越少。土样中部大孔径级别的孔隙数量明显高于表层，通过对比已有研究成果，表层碎石土的大孔隙主要是由于植物根系腐烂、干湿循环等因素造成，而颗粒粒径相关悬殊形成大孔隙骨架，经渗流作用后连成贯通通道是造成孔隙大于8 mm的主要原因。

(4) 大孔隙直径对成圆率影响较为明显。大孔隙直径与成圆率之间的平均相关性均大于0.9，碎石土的大孔隙成圆率与孔径级别之间存在一定的线性关系，即大孔隙直径越大，孔隙成圆率越低。

(5) 碎石土大孔隙的等效直径越大，整体轮廓系数越小，其轮廓的凹凸程度越明显。同一试验区距离相近的两个土样在同一深度上的形状参数差别不大，沿垂向深度的变化趋势也基本一致。