基于X-ray CT的古土壤孔裂隙识别与表征

李 鑫, 卢玉东, 张晓周, 卢阳春, 杨亚慧

(1.长安大学 环境科学与工程学院, 陕西 西安 710054; 2.旱区地下水文与生态效应教育部重点实验室, 陕西 西安 710054)

古土壤是黄土高原漫长的形成过程中在湿热环境中形成的棕褐色沉积物[1]。由于其黏粒含量高[2]，往往在黄土地区滑坡地质灾害研究中被当作隔水层[3-5]。按此观点，上覆黄土水分下渗过程中，当遇到古土壤层时会形成局部饱水带，使得黄土软化，强度降低[6-7]，形成滑动面，致使上覆黄土沿着古土壤接触面下滑[8-9]。但在对已经发生的黄土滑坡调查时，发现古土壤层往往被切穿，多数的解释是因为滑坡发生过程中，由于重力作用，使得古土壤层被切穿。然而，由于古土壤的黏性较高，在外界干湿交替、冻融循环、震动等过程中容易形成裂隙。如果是因为古土壤中存在的孔隙、裂隙优势渗流通道，使得水分穿过古土壤层，在其下形成更深层次的滑动面，也不无可能。因此，定量化研究古土壤孔隙、裂隙结构，有助于揭示古土壤中水分赋存及运移特征，对进一步开展黄土斜坡地质灾害机理研究具有重要的实践意义。目前在黄土—古土壤层序空隙微观结构研究方面，有大量针对马兰黄土、离石黄土的研究成果。高国瑞[10]利SEM和X射线衍射分析了黄土骨架颗粒形态、连接与排列方式，并提出了影响湿陷性的12种组合类型；雷祥义[11]利用压汞法将西安市附近黄土按照成因分为原生和次生孔隙，按尺寸分为特大孔隙(>250 μm)、大孔隙(16～250 μm)、中孔隙(4～16 μm)、小孔隙(1～4 μm)和微孔隙(<1 μm)；李晓军[12]在国内较早利用CT对压实黄土的微结构进行初步探讨；谷天峰[13]通过计算马兰黄土荷载前后SEM图像中孔隙微观结构参数的变化，揭示了土体变形过程中孔隙及颗粒变化特征；陈琼[14]利用氮气吸附法研究了滑带土微孔隙吸附特性；程亚楠[15]利用CT扫描的方法对黄土水力学参数进行了提取，利用编程建立了孔隙网络模型，并进行了水力学预测；潘网生[16-17]利用SEM对泾阳南塬黄土滑坡上样品进行扫描，获取黄土微观孔隙参数及非均质性和分形分维特征，并利用逾渗理论和网络分析方法对黄土微观渗流进行模拟；Li[18]利用SEM对马兰黄土进行扫描，并分析了黄土微观孔隙特征与湿陷性的关系；Li[19]利用CT对马兰黄土大孔隙特征进行研究，并对黄土大孔隙渗流进行了模拟；延恺[20]利用CT对黄土二三维微观结构进行表征，得到了团聚体中土颗粒三维微观形态特征和接触方式。另外，FIB-SEM[21]、核磁共振[22]等方法也开始应用于土体微结构的研究。上述方法基本都是针对黄土层序中马兰和离石黄土进行的研究，而针对古土壤孔隙裂隙特征和孔隙尺度水力特性的研究却少见报道。为此，本研究利用工业微焦点X射线计算机断层扫描技术(X-ray ICT)对泾阳南塬原状S5古土壤样品进行扫描，对古土壤中发育的孔隙、裂隙进行识别提取和定量化表征，为进一步进行古土壤水力学特性预测提供依据。

1 材料与方法

1.1 样品采集与处理

试验所用古土壤原状样品采集自泾阳南塬高庄镇修石渡村取土场剖面上比较明显的红3条S5古土壤层，为减少采样过程中对样品的干扰，首先利用刀和锯切割30 cm ×30 cm×30 cm的样品，用保鲜膜包装好后，利用模板箱搬运回室内。然后利用切割机进行加工，最终试验扫描样品尺寸为5 cm×5 cm×10 cm。样品干密度ρd为1.62 g/cm3，比重Gs为2.76，孔隙比e为0.7，饱和称重法测得孔隙度为41.18%，塑限16.18，液限32.39，黏粒含量33.80%，粉粒含量63.88%，砂粒含量2.32%，粒度不均匀系数Cu=30.11,曲率系数Cs=1.93，该古土壤为级配良好的不均匀土。

1.2 X射线计算机断层扫描(X-ray CT)

利用X-ray工业CT设备(图1)进行古土壤样品扫描，原理是利用X射线穿过样品中不同物质时具有不同的衰减系数，即致密组分衰减系数大，疏松组分衰减系数小，X射线穿过样品，由探测器接收到衰减程度不同的X射线，将其转换为电信号，在计算机上呈现为灰度值(0～255)不同的CT切片图像[23]，其中孔隙的灰度值最小，在图像上表现最暗，高密度物质灰度值最大，表现为最亮，中间灰度阶则为致密程度不同的固体土颗粒，通过沿样品竖向扫描可以获得一系列连续切片(图1)。本次CT扫描，试样尺寸为5 cm×5 cm×10 cm，空间分辨率为102 μm，共获得1 728张像素为1 024×1 024 pixel的二维图像。

图1 X-ray CT实验装置及扫描切片示意图

1.3 CT图像滤波处理与二值化分割

由于古土壤样品在CT扫描过程中仪器系统运行和环境噪声的影响，所得到的图像上存在噪点。为了获取高质量的图像，图像处理的第一步就是要对连续CT切片进行滤波处理，达到既消除噪点，又保留孔隙边缘的目的。AVIZO软件中提供了20多种滤波算法，按功能主要分为平滑滤波(如Median Filter，Bilateral Filter等)、边缘检测滤波(如Sobel Filter， Gaussian Filter等)，锐化滤波、频域转换滤波和灰度转换滤波[24]5大类，其中前3类在CT图像处理中较为常用。

本文经过对原始图像(以第330张切片为例)进行众多滤波算法(图2只展示部分滤波结果，包括中值滤波、高斯滤波、盒式滤波、双边滤波)计算比对发现，双边滤波器(Bilateral Filter)在降噪和边缘保留方面比其他滤波算法效果好，主要是因为双边滤波在对周边像素灰度值进行加权平均时，所采用的权重值综合考虑了像素欧几里得距离和像素局部灰度差异，使其同时具备降噪平滑和边缘保留的功能[25]。利用该滤波得到了既可以有效降噪，又能很好地凸显孔隙边缘的图像，为下一步进行图像准确分割提供了良好的基础。

图像阈值分割的主要目的是将样品中空隙空间和黄土颗粒分割开。本文利用阈值分割法(Threshold tool)和顶帽分割(Top hat)相结合，较为精准的区分出空隙和固体颗粒。其中，Threshold tool主要是基于影像灰度进行分割，但是分割过程中会出现阈值选取过小时，某些阈值边缘灰度较浅的孔隙被分割为颗粒，阈值过大时，又出现过度分割的现象；Top hat分割方法是基于黑顶帽或白顶帽算法计算相邻像素间灰度对比值来找出更“黑”或者更“白”的聚集像素，以识别出孔隙和土壤颗粒。Tophat识别出的空隙可以加入或者减除Threshold tool分割得到的孔隙或者颗粒中，以弥补或减轻因孔隙识别不足或者过度分割造成的偏差。通过两种工具的联合使用，对双边滤波法得到的灰度图像进行分割，小于设定阈值的部分为孔隙，赋值1以黑色表示，大于该阈值的部分为颗粒，赋值0以白色表示，由此获取了较为精准的二值化图像，如图2所示。

1.4 三维重建及可视化

基于移动立方体(MC， marching cubes)算法[26]，对上述阈值分割的系列二值化图片进行三维重建，得到孔隙和古土壤颗粒三维模型，为了更加直观地观察孔隙形态和三维展布情况，利用分水岭算法对三维孔隙体进行分割。基于光的发射和吸收原理，模拟预设光源通过孔隙体素并产生不同程度的发射和吸收，从而呈现出不同明暗程度的映射色系，达到可视化的目的。通过三维渲染效果图(图3)，可清晰地看到孔隙、裂隙的几何形态、空间展布、连接方式、连通程度等信息。

图2 图像滤波处理及二值化分割

图3 古土壤孔裂隙结构及基质颗粒

1.5 古土壤空隙参数计算及分类

为了对古土壤空隙结构进行定量表征，在图像分割及三维重建的基础上，首先计算了孔隙度(n)和形状因子(SF)，并利用三维空隙形状因子对古土壤中各类空隙进行分类，在分类的基础上，计算了不同类型空隙的等效孔径(d)。

孔隙度n为空隙体积(某切片内空隙面积)与土体(某切片)总体积(总面积)之比，见公式(1)：

(1)

式中：AV，AT——二维层面空隙面积和总面积(μm2)；VV，VT——三维空隙体积和总体积(μm3)。

形状因子(SF, shape factor)是表征三维空间几何体相对于球体褒瘪长短的形状特性，见公式(2)：

(2)

式中：A3D——三维空隙表面积(μm2)，V3D——三维空隙体积(μm3)。

等效孔径d是表征古土壤中空隙尺寸的参数，由于样品中空隙类型多样，而AVIZO软件提供的等效孔径计算方法是利用球形直径计算的，对非球形空隙 (如裂隙、枝杈状孔隙、长柱状孔隙)进行计算，难以真实准确表征所有空隙孔径特性。实际上，对于裂隙来说，等效孔径代表裂隙开度，在计算时应利用垂直于最大费雷特径的宽度计算；对于枝杈状孔隙来说，实际上是由若干长柱状孔隙交叉连通形成的，而长柱状孔隙局部为圆柱体，因此枝杈状孔隙和长柱状孔隙的等效直径应利用圆柱体体积公式来推求；对于椭球和似球体孔隙则可按照球形直径进行计算。本文对AVIZO中算法进行改进，针对不同类型空隙，应用不同算法来求取等效直径〔公式(3)〕，以更加准确地表征空隙特征：

(3)

式中：B3D——裂隙开度(μm)；V3D——孔隙体积(μm3)；L3D——孔隙长度(μm)。

2 结果与分析

2.1 二维层面结构定性分析

从图2中古土壤原始二维切片图像可以看出，图像以0～255的不同灰度组成，可以明显辨识出3类组分，其中黑色为孔隙、白色为高密度物质(主要为钙质成分)、中间灰色部分代表土颗粒。从图2中二值化分割结果可以看出，古土壤中孔隙和裂隙共存，呈现随机分散分布，形状各异，有狭长延伸型、似椭圆或圆型、不规则型孔隙，局部连通性较好，但是从整体上看，大部分都是孤立存在。

2.2 三维孔裂隙结构宏观定性分析

由于古土壤湿热气候中受环境改造明显，且黏性较高，质地较黄土脆硬，容易产生裂隙及孔道，孔裂隙相连接共同构成空隙网络。从三维数字模型中(图3)可以看出，古土壤中空隙在形态上主要由裂隙、枝杈状孔道、柱状孤立孔隙、椭球及似球状孤立孔隙组成，其中裂隙和枝杈状孔道连通性较好，为古土壤中水分运移提供优先通道，是古土壤中产生优先流的重要原因，很大程度上决定了古土壤的宏观渗透性质；长柱状孔隙虽然分散孤立存在，难以导水，但是在水分入渗时，裂隙和连通孔隙拓展发育过程中，可能会将该部分孔隙贯通，形成更大的范围的水力联系；孤立的椭球或球状孔隙通常为储气空间，一般起到阻水的作用，在极限条件下可能会以逾渗的方式发生水分交换。从发展演化的观点来看，裂隙发育部位在样品中上部，总体沿垂向延伸，开度随深度减小，表明裂隙与特大孔隙是孤立大孔隙拓展的驱动条件，后者是前者产生更大水力贡献的基础。通过以上定性描述，大致了解了古土壤中孔隙裂隙形态、赋存方式，为认识孔隙裂隙结构提供了良好的基础。但是要深入研究古土壤内部孔裂隙结构，还需要对其结构参数进行定量化表征。

2.3 孔隙度定量表征

在图像分割的基础上，由式(1)分别计算了Z，Y，X三个方向上的二维孔隙度值，得到了3个方向上连续切片的孔隙度分布曲线(图4)，以及三向孔隙度基本统计表(表1)。从图4和表1中可以看出，总体上，3个方向上孔隙度平均值为9.68%，与三维孔隙度相等；XZ和YZ平面上的孔隙度曲线形状和分布趋势相近，部分切片段内YZ孔隙度波动趋势比XZ滞后，表明有倾斜孔隙存在；XY平面上的孔隙度值与前两者差异较大，说明孔隙展布在三维空间上具有明显的垂向性或近垂向性。由表1知，XZ和YZ平面变异系数较接近，都小于10%，XY平面变异系数是前两者的2倍之多，由于样品中孔隙和裂隙并存，沿Z方向(XY平面)的孔隙多寡变化强烈，造成XY平面变异系数较大，以上现象都说明该古土壤空隙具有较强的垂向性和空间异质性。

图4 X，Y，Z三向切片孔隙度分布曲线

平面孔隙度/%平均标准差最小值最大值变异系数/%XY9.681.866.6714.1619.20XZ9.680.908.1112.059.27YZ9.680.756.8111.987.78三维9.68

2.4 空隙类型定量表征与分析

由图3直观地看出古土壤中空隙的形态，但是要将每一类提取出来，分类表征，需要借助三维形状因子加以区分，因此，依据公式(2)计算了所有空隙的形状因子，并利用自然断点法，将系列值分成了5类(图5)，即将前述五类空隙分别提取出来。当形状因子SF>40.18为裂隙，7.01

图5 基于形状因子的不同类型空隙分布频率及累积频率

2.5 孔径分布定量表征

在空隙类型划分的基础上，利用公式(3)计算了不同类型空隙的等效直径，得到古土壤孔径分布信息，并绘制了各孔径区间孔隙数量及其体积百分比柱状图(图6)。统计知,通过CT扫描识别出的最小孔径为0.109 mm，最大为5.059 mm，平均值为0.477 mm。从图6可以看出，样品中孔径主要分布在0.1～0.7 mm之间，该部分孔隙占到总孔隙数目的82.90%，0.7～1.0 mm占13.98%，1.0～5.0 mm占3.09%，大于5.0 mm的占0.02%，然而体积百分比方面，0.1～0.7 mm区间上孔隙体积占总孔隙体积的26.9%，0.7～1.0 mm占42.6%，1.0～5.0 mm占26.6%，大于5.0 mm的占3.9%。该结果与前述基于形状因子分类得到的孔隙分布和体积占比结果一致，说明基于形状因子对空隙类型进行分类表征可以达到对古土壤孔隙定量表征的目的。

图6 孔径分布与对应体积百分比直方图

3 讨 论

3.1 古土壤孔裂隙识别的水力学意义

通常认为，古土壤相对于黄土具有较小的渗透性能[4-5]，对垂向入渗的水分产生明显的阻滞效应，并引起水分的侧向渗透，形成局部饱和带，改变渗流场，降低坡体强度[16]。上述认识是基于野外和室内渗透试验建立的，渗流特性也是基于宏观达西定律或Richard方程描述的，如何在孔裂隙尺度上对古土壤水分入渗规律进行精细刻画，对于揭示古土壤在斜坡地质灾害中水力学贡献至关重要。本文无论从二维切片分割(图2)还是三维孔裂隙模型结果(图3)，都可以看出古土壤中发育规模不同、产状各异、形态多样的复杂空隙，具有较强的非均质性和各向异性，说明利用适用于均质假定的达西定律描述其渗流规律是有偏差的，因此有必要在孔裂隙尺度(微米尺度)上开展研究。本文采用新技术新方法对空隙类型进行了有效划分，并提出各类空隙的度量算法，得到准确的孔径分布结果。有助于进一步开展孔隙网络建模、孔隙裂隙尺度渗流模型的建立与渗流模拟，从而从本质上揭示古土壤在斜坡地质灾害中的重要作用。

3.2 分辨率对CT法测定孔隙度的影响

本文从CT扫描连续切片三维重建得到的孔裂隙结构，计算得到古土壤的孔隙度为9.68%，而利用饱和称重法计算的孔隙度为41.18%，压汞法计算得到的孔隙度为30.2%，3种方法测得的孔隙度值差异较大，主要原因是:饱和称重法测定时，由于土颗粒的亲水性，水分子直径0.4 nm，可以进入所有孔隙，因此得到的孔隙度最大，压汞法最佳测量范围为0.003～400 μm的孔隙[27]，对于特大孔隙测量不准确，因此比称重法偏小。而本文CT法测得的空隙受本次扫描分辨率影响，测得的最小孔隙为109 μm，最大孔隙达到5 mm，得到的孔隙仅仅为大孔隙和特大孔隙，因此可以看出，古土壤中等效直径小于100 μm的空隙占总空隙的比例较大。

在空隙识别与分类中，由于本次扫描分辨率仅为102 μm，对于微小孔隙难以识别，但是可以识别出绝大多数有水力学意义的大孔隙(d≥100 μm)[28]，为了系统性的对古土壤空隙特征精确识别，在今后应开展多尺度、多分辨率CT扫描，从微—介—宏观系统精确表征古土壤的内部空隙特性，为进一步研究古土壤孔裂隙双重介质渗流及其水力学特性研究提供科学基础。

4 结 论

本文为了进行古土壤孔隙裂隙结构识别与表征，优选双边滤波法对古土壤CT扫描图像进行滤波处理，进一步结合阈值分割和顶帽分割较为精准地区分了空隙和固体颗粒，利用移动立方体(MC)算法重建了三维孔裂隙模型，经过光学渲染得到了古土壤内部三维可视化结构，通过定性分析与定量表征，对于空隙二三维形态特征、连接方式、空隙类型与孔径分布规律等进行了分析。主要研究结论为：

(1) 依据二三维定性分析可知，古土壤孔隙裂隙结构并存，整体上呈分散分布，局部孔裂隙连通性较好；在形态上识别为裂隙、枝杈状孔隙、长柱状孔隙、椭球及球状孔隙；空间展布上具有明显的垂向性和空间异质性。

(2) 根据定量表征结果可知，古土壤大孔隙度(>100 μm的孔隙度)为9.86%，仅占总孔隙度(饱和称重法孔隙度41.18%)的23.4%，说明古土壤以等效直径在100 μm以下的孔裂隙为主；形状因子是有效提取与识别孔裂隙结构的关键参数，依其可对空隙进行精确分类表征；可识别孔径分布在0.1～5.0 mm之间，以0.1～0.7 mm孔隙居多，而0.7～1 mm孔隙体积贡献率最大。

(3) CT无损扫描及三维可视化重建技术可以有效识别与表征古土壤内部孔裂隙结构，为进一步开展孔隙尺度渗流模拟提供可靠基础。