RVE孔隙模型细观结构特征分析与对比

杨鹏飞，李 显，阙 云，蔡沛辰

(福州大学 土木工程学院，福州 350108)

1 研究背景

我国东南沿海地区，花岗岩残积土分布非常广泛，其不同于一般黏性土，属于典型大孔隙结构[1]。除花岗岩残积土自身孔隙结构外，在自然气候等条件影响下，斜坡表面土壤内分布着许多孔隙，这将进一步改变土体孔隙结构特性[2]。在降雨等作用下水流可通过孔隙绕过大部分基质土壤，产生孔隙流，导致边坡失稳[3]。因此，如何准确刻画土壤内孔隙分布情况是揭示土壤边坡失稳机制的基础与核心内容之一。

不同学者对于界定大孔隙的标准存在一定差异，石辉等[4]定义大孔隙为等效直径在0.2～3.0 mm范围内的孔隙；时忠杰等[5]定义大孔隙直径为0.4～2.3 mm；Marshall[6]将直径>30 μm的孔隙定义为大孔隙；Jarvis[7]将直径在0.3～0.5 mm的孔隙定义为大孔隙。本研究在前人的基础上，考虑工业CT扫描仪的精度，将等效直径大于最低分辨率0.15 mm的孔隙定义为大孔隙。

目前，土壤大孔隙的研究方法多为染色示踪法[8]、穿透曲线法[9]等，但这些方法对观测的土体试样破坏性较大，无法得到大孔隙在土壤中的三维结构参数。相比之下，CT扫描技术对原状土大孔隙研究更加便捷、精确，可在不破坏土体原始结构的条件下，提取大孔隙的三维结构特征参数。Anderson等[10]最早利用CT扫描仪对干土和湿土进行扫描；冯杰等[11]采用CT扫描得到了大孔隙在截面上的分布情况；薛华庆等[12]采用微米CT扫描表征了页岩的发育情况；戚楠[13]采用不同的孔隙结构表征技术对储层孔喉进行了定量表征和二维平面刻画。但通常，对大孔隙空间结构参数的计算往往还需将扫描图像作进一步处理，常用的图像处理软件有AVIZO、Image J等。现有的大孔隙结构研究已取得了一定进展，但由于量化技术、地域性差异等原因还存在一定问题。很多学者对大孔隙的研究多集中在二维结构参数上，且研究多针对农业土壤成岩石结构，很少有系统讨论三维结构特征参数的影响，并且在不同的研究条件、研究区域的影响下得出的结果也不尽相同。同时，现有的大孔隙结构参数获取过程中存在计算机运算速度慢、所需内存容量大等特点。

鉴于此，本文从三维层面分析原状花岗岩残积土大孔隙分布特征，基于工业CT扫描图像和AVIZO软件构建了孔隙网络模型，选取代表性体积单元(Representative Volume Element，RVE)，对花岗岩残积土大孔隙三维细观结构进行了定量表征，并对比了不同孔隙处理方法下的孔隙结构特征。研究成果可为短时强降雨下含大孔隙斜坡的水分运移机制研究奠定基础，并可进一步丰富短时暴雨型滑坡理论基础及支撑灾害防治。

2 试样与方法

2.1 原状土取样

原状土取自福州某地山坡，分别选取两个试验区A和B，沿深度范围进行取样(A区试样1和B区试样3在深度20 cm处取样，A区试样2和B区试样4在深度70 cm处取样)，最终获取4个原状土样，如图1[14]所示。 测试土样的基本物理参数，如表1[15]所示。

图1 原状土取样[14]Fig.1 Sampling of undisturbed soil

表1 基本物理参数[15]

2.2 CT扫描试验

在原状土大孔隙结构特征观测方法中，CT扫描属于无损检测方法，可在不破坏土体结构的情况下获得土体内部孔隙结构图像。其他学者研究的结果表明，在扫描土壤横断面图像质量方面，工业CT的图像质量要明显高于医学CT图像[16]。在三维重构方面，工业CT解析出的三维图像能更好地体现出土壤孔隙分布。鉴于此，对获取的原状土样进行工业CT扫描，设备名称为C450 kV高能量工业CT，将扫描试样放置在该CT系统工作台的中央，工作电压为450 kV，电流为63 mA，扫描最低分辨率为150 μm。得到典型CT扫描图像如图2所示。

图2 典型CT扫描图像Fig.2 Typical CT scanning images

2.3 三维重构预处理

本次研究采用VG studio MAX 3.0软件对CT图像进行三维重构，并结合AVIZO软件进而提取孔隙面积、孔隙等效半径、孔隙配位数、喉道面积、喉道等效半径、喉道长度等参数。

由于收集到的CT图像与真实情况有一定的差距，因此在三维重构前需要对CT图像进行归一化和环状伪影去除等操作。其中归一化主要用来纠正扫描过程中的光束强度变化，能对不同时间收集的CT投影图像的强度进行补偿，使强度的变化充分均匀。环状伪影是CT投影成像中较为常见的一种现象，使用平面场校正功能可以有效去除环状伪影。

CT扫描得到的原始切片需后续进一步处理才能进行细观结构识别。VG studio是目前最常用的图像三维重构软件之一，主要功能为显示、处理和分析工业CT扫描后得到的图像数据。由于本次试验着重于研究花岗岩残积土大孔隙的三维形态，在较高的分辨率情况下，细微的扰动对试验结果可能造成较大的影响，因此需要选择适当的分析区域，去除各类干扰，尽可能还原花岗岩残积土大孔隙的三维形态。样品尺寸为15 cm×15 cm×40 cm，选取的分析区域为13 cm×13 cm×36 cm。

2.4 三维模型构建

CT图像三维重构是采用连续的CT图像通过插值运算生成三维立体结构的一个过程[17]。本次试验根据文献[14]、文献[18]、文献[19]中的方法，选取了300×300×300体素的代表性体积单元(RVE)，其中体素是数字数据于三维空间分割上的最小单位。相邻切片之间的间距与CT图像的分辨率相同，即每隔0.15 mm扫描一张CT图像。

AVIZO可视化软件是一套全方位的三维数据处理软件，可对三维图像进行阈值分割、曲面重建、网格划分、孔隙网络模型建模等操作[14]。对图像进行三维分析时，采用拟合法[20]，在每个试样中随机挑选10张CT图像导入Image pro-plus软件，然后利用Image Histogram功能提取出图像灰度值，绘制灰度值分布图，并通过软件进行拟合，得到每张图的阈值，对阈值求平均值，即可得到4个试样的阈值分别为15 470、16 505、19 920和16 528。

图3为模型中间位置处4个试样三维图像采用中值滤波器降噪后的阈值分割图，其中蓝色部分代表孔隙域，灰色部分代表基质域。将分割后的三维图像进行对象分离后，对孔隙结构进行分析，即可得到三维孔隙的特征参数。图4为4个原状试样的代表性体积单元三维重构结果，图5为删除孤立孔隙后的简化孔隙网络模型。

图3 原状土阈值分割Fig.3 Threshold segmentation of undisturbed soil samples

图4 原状土三维重构模型Fig.4 Three-dimensional reconstruction model of undisturbed soil

图5 简化孔隙网络模型Fig.5 Simplified pore network model

3 结果与讨论

3.1 孔喉特征参数分析

采用孔隙模型分析大孔隙结构特征时，主要将孔隙空间划分为孔隙结构和连接孔隙之间的喉道结构。因此孔喉的尺寸特征分析主要包括孔隙数量、喉道数量、孔隙结构参数和喉道结构参数。孔隙特征参数见表2，喉道特征参数见表3。

从表2可以看出：不同试样的连通孔隙数量差值较大，试样3的孔隙数量是试样1的2倍左右，而不同试样的孔隙半径差值较小。图6为4个试样的孔隙平均半径分布，可以看出4个试样连通孔隙的半径多介于1～4 mm之间。

从表3可以看出：4个试样的喉道数量差距较大，其中试样3的喉道总数为3 765个，而试样1的喉道总数仅有819个，说明不同试样之间的孔隙连通情况差异性较大。4个试样的喉道等效半径、喉道长度和喉道表面积平均值差异性较小，表明在花岗岩残积土中虽然孔隙数量差异性较大，但由于孔隙形成原因类似，孔隙结构之间差异性较小。

由2.4节可知，本文REV模型尺寸为300×300×300体素，而文献[14]中为大尺寸三维重构模型，其相比于小尺寸模型，能够更真实地反映土体的原始结构，但同时计算量更大，运算也更加耗时。再对比二者的孔隙特征参数结果发现，4个试样的孔隙平均半径、喉道长度和喉道表面积相差较小(平均误差分别为4%、0.8%和2.9%)，孔隙表面积和喉道等效半径平均误差在5%和17%左右，分析原因可能是原状土孔隙三维结构较为复杂，选取区域与整体区域结构特征存在不同。此外，由于REV选取的自身特点使得二者在孔隙配位数方面差异性较大[14]，但总体REV选取的方法仍可在一定程度上对孔吼结构参数进行分析。

图6 孔隙半径分布Fig.6 Distribution of pore radius

3.2 孔隙拓扑结构量化表征

土壤孔隙空间的拓扑结构是指土壤内部孔隙与喉道连接关系，即连通率[21]。孔隙拓扑结构量化可由孔隙配位数来表征。由表2可知，4个试样的孔隙配位数最大值分别为24、31、47、38，经过统计发现花岗岩残积土中孔隙配位数主要集中在20以内，超过20的孔隙配位数分布呈离散点状，且数量皆为个位数，其与大尺寸模型统计结果较为相近[14]。因此本文在研究孔隙配位数分布情况时，只统计20以内的数值，并将其绘于图7。

表2 孔隙特征参数

表3 喉道特征参数

从图7可以看出：4个试样的孔隙配位数分布差异性较大，其中试样1中配位数为1～4的孔隙构成了土壤孔隙的主体，占总孔隙数量的65.52%；试样2中配位数为1～4的孔隙构成了土壤孔隙的主体，占总数量的69.93%；试样3中配位数为1～7的孔隙构成了土壤孔隙的主体，占总数量的91.73%，试样4中配位数为1～6的孔隙构成了土壤孔隙的主体，占总孔隙数量的89.22%。上述现象表明：在花岗岩残积土中，不同位置的土壤内部孔隙拓扑空间结构差异性较大，而深度对其影响较小，如试验区B的不同深度试样3和4的孔隙连通性都较好，而试验区A中试样1和试样2的孔隙连通性都较差。

图7 孔隙配位数分布Fig.7 Distribution of coordination number of pores

同时对比文献[14]中大尺寸模型，小尺寸试样1、2、4中构成土壤孔隙主体的配位数差异不大，而试样3中获得的配位数1～7的孔隙占孔隙总数量比例偏大，分析原因可能是由于REV模型试样3相比大尺寸模型孔隙数量较少，整体分布不均匀所致。综上，REV模型可一定程度上表征孔隙的空间拓扑结构。

3.3 不同处理方法孔隙参数的差异性分析

为对比不同处理方法下大孔隙特征参数的差异性，本文将二维图像处理[22]和三维重构处理分别提取出的孔隙参数平均值汇总于表4。为减少工作量，二维图像每隔1.5 cm选取一张进行孔隙参数分析，三维重构处理可分为孔隙模型和孔隙网络模型两种。

表4 2D和3D孔隙参数平均值

从表4可以看出：采用二维图像处理得到的孔隙平均直径最小，而且不同试样之间差异性较大，而通过删除孤立孔隙而建立的孔隙网络模型得到的孔隙平均直径最大。同时，对比三维孔隙模型和孔隙网络模型处理结果可知，在孔隙网络模型中通过等效半径计算出4个试样的平均孔隙体积分别为33.36 、42.68 、47.27 、49.18 mm3，通过孔隙模型计算出4个试样的平均孔隙体积分别为6.28、6.05、5.57、8.34 mm3，说明在土壤内部有大量的孤立孔隙存在，由于孔隙体积越小，和其他孔隙产生连通的几率也越小，因此孤立孔隙的体积一般远小于连通孔隙。当采用孔隙网络模型分析孔隙参数时，孤立孔隙被删除，因此造成2种统计方式平均孔隙体积相差5倍以上。孔隙的平均表面积和平均孔隙体积类似，在三维孔隙模型中4个试样的平均孔隙表面积分别为16.47、15.87 、14.80 、19.11 mm2，而通过孔隙网络模型计算出4个试样的平均孔隙表面积分别为136.38 、178.46 、249.98 、267.31 mm2，可以看出通过2种方式统计出的平均孔隙表面积差距比平均孔隙体积更大，原因在于随着孔隙的增大，孔隙表面复杂程度也随之增加。

本研究中对于大孔隙的分析采用的是三维孔隙网络模型处理方法，原因在于：二维处理方法是对众多CT扫描切片进行参数分析，只能从二维层面体现孔隙特征，而事实上，土体中的大孔隙是三维结构，二维表征的精确度有待商榷；三维孔隙模型处理方法虽是从三维层面对土体大孔隙结构进行表征，但统计时将各种孤立孔隙都计算在内，而这些孤立孔隙一般不参与渗流过程，对于后续渗流研究不能起到很好的指导作用；三维孔隙网络模型处理方法是在三维孔隙模型的基础之上，检测连通孔隙和孤立孔隙，并将孤立孔隙删除之后形成的连通孔隙网络模型，这样可以有针对性地统计有效孔隙结构特征参数，还可以提高后续渗流模拟研究中的效率。

4 结 论

(1)花岗岩残积土RVE模型中不同试样的连通孔隙数量差异性较大，而孔隙结构参数差值较小，且孔隙半径多介于1～4 mm之间，表明在花岗岩残积土中虽然孔隙数量差异性较大，但由于孔隙形成原因类似，孔隙结构之间差异性较小。

(2)不同区域的土壤内部孔隙拓扑空间结构差异性较大，部分孔隙的配位数可达40以上，但深度对其影响较小。其中孔隙配位数主要集中在20以内，超过20的呈离散点状分布。

(3)对比不同处理方法知，2D图像处理得到的孔隙平均直径最小，孔隙差异性较大，而通过删除孤立孔隙而建立的孔隙网络模型，得到的孔隙平均直径最大，孔隙分布较为均匀。对比3D处理结果表明，花岗岩残积土内有大量的孤立孔隙存在，而孤立孔隙的体积一般远小于连通孔隙。

(4)明确土壤内部孔隙特征分布是揭示土壤边坡失稳机制的基础，本文对三维层面花岗岩残积土内大孔隙分布特征进行了定量表征，考虑到土体原始三维结构计算复杂，采用REV孔隙模型代替大尺寸模型研究，其结果能较为精确地表征土壤孔隙结构，但由于土壤内部结构复杂，整体结构会与选取区域存在部分差异。REV模型可为小尺寸原状土定量化研究孔隙内部水分运移机理提供一定借鉴。