基于CT的黄土孔隙尺度优先流特性

2022-09-29宁瑞浩冷艳秋何芝远李泽坤马哲

科学技术与工程 2022年23期

宁瑞浩，冷艳秋*，何芝远，李泽坤，马哲

(1. 长安大学地质工程与测绘学院，西安 710054； 2.西部矿产资源与地质环境教育部重点实验室，西安 710054)

黄土在中国主要分布在山西、陕西、甘肃、宁夏的黄土高原一带，是一种疏松多孔、竖向节理、裂隙发育，具有很强的湿陷性和水敏性的结构性土；因其特殊的湿陷性和水敏性，水便成为了黄土地区诱发地质灾害的关键因素。常规的渗水试验表明降雨或者灌溉的入渗深度只有2～3 m[1]，但是黑方台和泾阳等地的滑坡实例却显示水运移到坡体以下数十米，优先流便是这类地质现象的合理解释，因此探索黄土孔隙、裂隙中优先流的渗流特性，对防控黄土地质灾害具有重大意义。

优势流又称为优先流，是一种非平衡流，主要特征为部分水流绕过大部分土体基质沿着优势通道迅速入渗到土体深部[2]。目前对优先流的研究主要集中在两个尺度，即宏观尺度和孔隙尺度[3]，Watanabe等[4]通过构建局部非孔洞模型，研究了碳酸盐岩的孔隙度与渗透率和内部优势流的关系，为碳酸盐岩储层的开发利用提供了依据。赵宽耀等[5]通过物探探测与数值模拟相结合的方法，分析了黄土边坡在优势流与基质流共同作用下不同灌溉强度水的渗流过程。陈思婕等[6]将双重渗透模型与无限边坡稳定分析法相耦合，根据土壤含水量及孔隙水压力观测数据，模拟了强降雨条件下滑体内的水动力过程，以此分析了优先流对滑坡触发机理的影响。然而这些研究都是集中在宏观层面的。潘网生等[7]认为定量化描述黄土微细观孔隙、裂隙结构模型,建立从微细观到宏观不同尺度下研究黄土滑坡机理的系统理论方法，对认识优先流对滑坡的触发机理尤为重要，仅仅开展宏观尺度的研究并不能揭示黄土中优先流的本质规律。

研究孔隙尺度内的优先流特性首先要获取真实的黄土内部孔隙结构，在众多方法中计算机断层扫描(computed tomography,CT)技术作为一种无损检测技术正被广泛用于岩土材料的孔隙定量化表征。Li等[8]利用CT扫描技术定量化研究了马兰黄土的大孔隙结构特征，结果表明马兰黄土在垂直和水平方向上的孔隙形状、联通特征有很大差异，这与室内试验与数值模拟结果较为吻合。王超等[9]采用微米CT建立了黄土岩的孔隙网络模型，研究了黄土岩中矿物成分、形状对黄土岩的渗透率各向异性的影响。Wang等[10]通过CT扫描技术研究了黄土区露天煤矿排土场重构土壤孔隙三维分布的多重分形表征，为矿区土壤重构提供了依据。上述成果说明利用CT技术获取黄土内部真实孔隙结构并进行定量化的研究是可行的。

目前已有众多学者在真实孔隙结构的基础上开展优先流渗流特性研究，Li等[11]、Larsbo等[12]在孔隙尺度上研究了孔隙与孔隙网络特征(平均直径、配位数等)对土体内部优势流渗流规律的影响。Kim等[13]通过数值模拟的方法研究了孔隙尺度内河床附近优势流及湍流对溶质异常输运特性的影响。Zhang等[14]采用低场核磁共振技术对纳米颗粒在非均质多孔介质中的运移方式实时监测，观察到了明显的优势流现象。鲁拓等[15]研究了基于多种分形模型，计算了不同尺度范围孔隙的分形维数，并讨论了二者之间的联系。Sarah等[16]利用高分辨率的微米CT研究了土壤内部的细微观结构并利用贝叶斯统计方法探讨了土壤的微观特征与其饱和导水率等参数的相关性。Li等[17]利用AVIZO软件研究了马兰黄土的大孔隙特征并基于形状因子对其进行分类，为研究孔隙尺度内优势流特征提供了依据。

前人所做研究采用CT扫描的分辨率多数大于50 μm，雷祥义[18]根据压汞实验的结果将大孔隙划分为孔隙半径为16～250 μm，中孔隙为孔隙半径为4～16 μm，采用50 μm以上的分辨率显然并不能满足精度需求，因此在提取黄土内部真实孔隙结构的基础上开展微细观尺度的优先流渗流特性研究。采用高精度的微米CT，扫描分辨率为4 μm，以此构建真实孔隙结构，研究延安黄土孔隙尺度内的优先流渗流特性。

1 实验材料与方法

1.1 实验材料

实验所用土样取自陕西延安边坡坡脚，取样深度约为75 m，所取原状土样为圆柱体，直径和高度分别为25 cm和50 cm，在密封和标记后运至实验室，运输过程中采取合理的避震措施从而避免扰动。通过对所取试样随机测量，研究区黄土试样的基本物理性质见表1。图1为试验所采用黄土试样的粒度分布曲线，其中粉粒(2～75 μm)含量81.72%，黏粒(<2 μm)含量为6.65%。结合粒度分布曲线及塑性指数，按照文献[19]对黄土的分类，研究区黄土属于黄土质砂黏土。

表1 土样基本物理性质Table 1 Basic physical properties of soil samples

图1 粒度分布曲线Fig.1 Particle size distribution curve

1.2 微米CT扫描

实验所采用CT扫描设备为Phoenix v|tome|x s工业CT(图2)，该CT系统组合有180 kV/15 W 高功率nanofocus X射线管和240 kV/320 W 的微焦点管，是一种灵活可靠的扫描分析工具。本次测试试样为圆柱体，直径7 mm，高15 mm(图3)，扫描分辨率为4 μm，共获取X、Y、Z方向切片14 200张，每张二维切片像素点的灰度值代表着该点的X射线衰减系数，密度越大衰减系数越大，借助这一原理再结合适当的分割方法，三维重建等手段，便可以把试样内部的孔隙裂隙识别重建出来以用于二维、三位参数的提取，具体技术路线见图4。

图2 CT扫描设备Fig.2 CT scanning equipment

图3 扫描土样Fig.3 Scanned soil sample

图4 技术流程图Fig.4 Technical flow chart

1.3 室内渗透实验

为研究黄土微细观与宏细观渗流之间的联系并验证计算结果可靠性，进行室内渗透实验与数值计算结果进行对比。本次室内变水头渗透实验采用土样为延安黄土原状土样与重塑土样(图5)，各自设置5组平行试验，重塑土样采用分层静压法制备，将原状黄土进行碾碎、过筛后放入涂好凡士林的压样筒内，控制重塑土样干密度与原状样相等，分五层压实，层间接触部位进行刮毛，每层压实时间不少于45 min。试样为标准尺寸，直径、高度分别为61.8、40 mm，在真空缸中饱和后进行渗透实验。渗透实验流程按照《土工试验方法标准》(GB/T 50123—2019)规定，待水头管中水头稳定后在起始时刻标定水头高度，每隔20～40 s测定水头和时间的变化，如此连续测记10次后再使水头回升至适当高度，重复试验5次，所得结果如表2所示。

图5 实验土样Fig.5 Experimental soil sample

表2 室内渗透实验所得渗透系数Table 2 Permeability coefficient obtained from indoor permeability test

2 CT图像预处理与三维重建

2.1 建模软件选取

在目前常用的多孔介质二、三维孔隙结构及颗粒定量化分析平台中，AVIZO软件可以满足图像降噪处理-图像二值分割-三维体重建-数据分析-渗流模拟-绝对渗透率计算等各方面的需求，在科学和工业领域、生命科学和生物医学、材料科学、土壤学、岩土领域有着广泛的应用，能够满足本研究的相关需求，因此选用AVIZO平台进行优先流的相关研究。

2.2 图像预处理

在CT图像的采集和获取过程中，会因为采集设备和仪器的少量干扰和周围环境因素的影响，使得采集的图像或多或少的存在部分噪声，为了进一步提高扫描图像的质量，突出图像中孔隙的边缘轮廓特征，以便能够准确分割目标孔隙和背景，提取出感兴趣的孔隙部分，通常采用图像噪声滤波方法和边缘增强进行预处理。

目前常用的滤波方法有高斯滤波、中值滤波、非局部均值滤波法等，图6分别为采用上述滤波方法的预处理结果。由图6可知，非局部均值滤波法几乎能完全去除图像中的高斯噪声、脉冲噪声等，同时还能够很好地保留图像细节，很大程度地减轻边缘模糊[图6(d)]；其他方法虽然能够消除噪声但是不能很好的保留图片信息，致使边缘模糊，因此选用非局部均值滤波法对所获取的CT图像进行降噪处理。非局部均值滤波法的基本思想为根据图像的自相似性来计算邻域像素权重，通过加权的形式将最相近的几个像素块中的中心点结合起来估计真实值，缺点是对计算机的计算能力要求较高；虽然采用的滤波方法在有效去除噪声的同时能够很好地保留图像细节，但仍然存在一定的边缘模糊现象，即黄土固相和孔隙之间交界处的过渡带[图6(e)]。采用AVIZO软件中的unsharp mask模块对图像进行边缘增强处理[图6(f)]，消减固相和孔隙之间交界处的过渡带，大大增加分割精度。

2.3 图像分割

在黄土孔隙结构的三维重建及要素提取研究中，针对黄土CT扫描图像的分割技术是整个研究的核心与关键点之一，其分割结果直接影响到后续孔隙结构三维模型重建的精确性。到目前为止还没有任何一种完备的分割方法能够适用各种应用环境，常用的图像分割方法可以划分为:基于阈值的分割方法、基于区域的分割方法、基于边缘的分割方法以及基于某些特定理论的分割方法等。本研究采用基于阈值的分割方法。阈值分割法的基本原理为给定一个恰当的灰度阈值，图像某点处的灰度值大于该阈值认为是土颗粒，小于该阈值则认为是孔隙，因此阈值的选取尤为关键，过大会导致过分割，过小则会导致欠分割；常用的阈值分割方法有最大熵法、Ostu法、矩量法等，对同一张灰度图像采用不同的分割方法的局部放大图如图7所示，最大熵法与矩量法确定的阈值偏大，存在明显的过分割现象，Ostu分割法的分割效果与实际情况最为接近，然而还存在轻微的过分割现象，需手动调试阈值。基于此，分割阈值采用Ostu法确定参考阈值，在AVIZO软件中通过手动调整确定最佳阈值[图7(d)]，对部分难以识别小孔隙、微小裂隙通过AVIZO软件中的Top-hat算法进行识别[图7(e)]，最终取二者并集[图7(f)]作为最终孔隙分割结果。经计算分割后的微观孔隙度约为19%，小于实际孔隙度，这是由于小于4 μm的孔隙无法被识别出来。

图6 图像预处理结果Fig.6 Image preprocessing results

2.4 孔隙结构三维重建

2.4.1 代表性体元确定

目前，基于CT扫描图像三维重构的基本方法主要包括体绘制和面绘制两种方法，采用可以避免图像信息丢失，能够保留数据完整性的体绘制方法对土样中的孔隙网络进行三维重构。考虑到计算机的计算能力有限并且必须兼顾样品具有够的代表性，需确定代表性体积单元(representative elementary volume，REM)[20]。本文研究确定代表性体积单元的方法为在土样任意一点处取一系列不同边长的立方体(图8)，计算其孔隙度，当孔隙度趋于稳定时即认为该边长的立方体为最小代表性体积单元。如图9所示，在立方体边长达到1 250 μm时，代表性立方体孔隙度已达到稳定值，即立方体边长达到1 250 μm以上时即可认为该子体积具有代表性。

图7 图像分割结果Fig.7 Image segmentation results

图8 REV选取示意图Fig.8 Schematic diagram of REV selection

目前对优先流的常见分类为：裂隙流、管流、指流、侧向流，指流和侧向流分别发育在上细下粗和含有较大阻水块体的土体中，而本次所取土样土质较均一且不含岩石、树根等阻水块体。基于此，在扫描土样中分别提取含裂隙(C1、C2，图10)、含管状通道(C3、C4，图10)与不含优势通道的(C5、C6，图10)2 000 μm×2 000 μm×2 000 μm的立方体，排除孤立孔隙，在各立方体连通孔隙基础上对比分析裂隙流、管流的渗流特性。

图9 孔隙度与立方体边长关系曲线Fig.9 Relation curve between porosity and cube size s

图10 各代表性立方体Fig.10 Representative elementary volume

2.4.2 隙网络模型的构建

孔隙网络模型(pore network model,PNM)的基本假定为一个相互连通的孔隙网络有多个孔隙组成，不同孔隙之间通过孔道连接，孔喉是孔道中横截面积最小的地方。在AVIZO中将孔隙网络模型用球棍模型来表示，分别代表孔隙与孔隙之间的联通关系，通过构建孔隙网络模型可以获取孔隙之间的孔道长度、孔喉面积以及配位数等关键信息，孔隙网络模型如图11所示。

图11 孔隙网络模型Fig.11 Pore network model

2.5 重建结果分析

2.5.1 孔隙与孔喉分布特征分析

根据所建立孔隙网络模型所得数据，孔喉和孔隙等效直径的分布范围基本一致，对孔隙和孔喉的等效直径以10 μm为间隔进行直方图统计，发现孔隙与孔喉近似服从对数正态分布。根据图12，孔隙主要分布在孔径为40～140 μm的区间内，占总数目的91.83%，孔喉主要分布在孔径为10～50 μm的区间内，占总数目的85.07%。

图12 孔隙与孔喉分布特征Fig.12 Distribution characteristics of pores and pore throats

2.5.2 孔隙连通特征分析

孔隙的连通性可以用配位数表征，一个孔隙的配位数是指与其他孔隙连接的孔喉的数量，配位数越多则代表着与周围的孔隙连通性越好，根据构建的孔隙网络模型，对获取的孔隙配位数进行统计，发现配位数近似服从对数正态分布(图13)，其值介于1～28之间；含优势通道的立方体中平均值为7.17，众数为6，不含优势通道的立方体中平均值为6.15，众数为5，这说明在含优势通道的立方体中，孔隙的整体连通性要好于不含优势通道的立方体。

图13 孔隙配位数分布特征Fig.13 Distribution characteristics of pore coordination number

3 渗流模拟

3.1 基本原理

Navier-Stokes方程是用于描述不可压缩黏性流体的动量守恒方程，其矢量形式为

(1)

式(1)中:ρ流体密度，kg/m3；V为速度矢量；P为流体压力,kPa；f为单位体积流体所受外力，kN/m3；μ为流体动力黏度,Pa·s。对V有

(2)

根据式(2)将式(1)展开，得

(3)

一般条件下很难得到其精确解，考虑以下条件将Stokes方程简化：流体为不可压缩流体；流体为牛顿流体，动力黏度μ为常数；流体为稳定流，流速不随时间改变；流体流动方式为层流；对黄土孔隙尺度内流动的水，上述条件均可满足。

简化后的Stokes方程为

(4)

为保证方程在整个体积的有效性，将其转化为体积平均形式，使方程在体积空间内更加平滑，即

(5)

式(5)中:D为速度扰动场张量；d压力扰动场向量；I单位张量。

通过求解系统体积v上速度扰动场张量的平均值来求得渗透率张量k，即

(6)

在AVIZO中用有限体积法求解方程组，方程在交错网格上离散，假定体素为立方体，待求解压力位于体素中心，而速度在体素面上分解，待求解量的时间导数趋于零时认为收敛，误差计算方法为

(7)

式(7)中:n为当前迭代次数；∂t为时间增量；c2为虚拟压缩系数[21]。

对于本次计算，当误差小于10-4时即认为收敛，当求得绝对渗透率后，渗透系数计算公式为

(8)

3.2 模拟结果

模拟渗流方向为自上而下，上下边界为定水头边界，四周为不透水边界，上下边界的水头差为0.002 m，即代表性立方体的高度，各立方体内部流速分布如图14～图16所示。根据图14(a)(e)～图16(a)(e)，在含优势通道的代表性立方体中(C1、C2、C3、C4)高流速的渗流路径沿着裂隙展布方向集中，裂隙、管状通道内的流线更加顺直平滑，渗流压力大，而其周围的大孔隙内渗流路径则较为分散，渗流压力也更小，压力分布更加不均匀；而在不含优势通道立方体中(C5、C6)，流线及压力的分布较为均匀，流速也更慢，表现为均匀入渗,各代表性立方体渗透系数及绝对渗透率如表3所示。

表3 渗流模拟所得渗透系数Table 3 Permeability coefficients obtained from seepage simulations

图14 含裂隙代表性体元模拟结果Fig.14 Representative voxel simulation results with fractures

综合室内渗透实验及数值计算结果可知重塑黄土由于内部结构被完全破坏导致孔隙之间连通性极差，渗透系数最小，要比原状黄土小近一个数量级。不含优势通道的代表性立方体渗透系数与原状黄土室内渗透实验所得渗透系数在同一个数量级上，含优势通道的代表性立方体渗透系数要比原状黄土大近一个数量级，这与前人所得结果一致[22]，说明优势通道对黄土的渗透性起着主导作用，也间接证明了黄土微细观尺度与宏细观尺度渗流的统一性。

为详细研究黄土内部优先流的分布及流速特征，借助前述对黄土内部孔隙进行概化的孔隙网络模型将孔喉按流量大小进行展示[图14(b)(e)～图16(b)(e)]，通过孔喉截面面积及流量换算流速并分组统计，依据数值计算及室内渗透实验所得渗透系数将孔喉按流速划分为三部分，不流动-低流速孔喉(T1)：流速小于1 μm/s，多位于孔隙连通性较差部位，流速与重塑黄土渗透系数相当，处于同一数量级；均匀流速孔喉(T2)：流速介于1～100 μm/s，距离优势通道较远，但孔隙连通性一般部位，与原状黄土及含优势通道的代表性立方体渗透系数处于同一数量级；优先流发育孔喉(T3)：流速大于100 μm/s，多位于优势通道及孔隙连通性较好部位。

图15 含管状通道代表性体元模拟结果Fig.15 Representative voxel simulation results with pipes

对三个流速段的孔喉数量进行统计，统计结果见图17。三种立方体中均以均匀流速孔喉为主，其中无优势通道立方体中的T2占比最高，含裂隙立方体与含管状通道立方体中的T2通道数量几乎相同。将含裂隙立方体与含管状通道立方体进行对比分析发现前者T3占比为30%，其含量超过T1，同时也比后者T3占比高，因此含裂隙立方体的绝对渗透率比含管状通道立方体大。将含裂隙与不含优势通道的代表性立方体孔喉流速分布进行对比，发现裂隙发育的立方体中T1所占百分比明显降低，T2数量约减少10%，而T3所占比例显著增加，约增加15%，说明裂隙对整体的渗透特性影响为均匀促进型，即裂隙的存在会整体地提升立方体的渗透性；而管状通道发育的立方体与不含优势通道的立方体相比T1占比约增加5%，T2数量约减少10%，T3所占比例增加幅度较裂隙发育的立方体小，但管状通道内的最大流速显著高于裂隙内部，这说明管状优势通道对整体的渗透特性影响为渗流集中型，即渗流路径沿着管状通道集中，对非优势通道在渗流过程中所发挥的作用有所削弱。在不含优势通道的代表性立方体中仍有相当一部分孔喉流速大于100 μm/s，这说明局部联通性较好的大孔隙也可形成优势通道。对三种流速段孔喉直径进行统计发现T1部分孔喉平均直径为11.2 μm, T2部分孔喉平均直径为19.4 μm, T3部分孔喉平均直径为29.1 μm,由此可见孔喉流速与孔喉直径关系为正相关。