基于最大球算法定量表征土体空间孔隙网络*

解朝阳 张 巍 姚东方 王 坚 朱雨辰

(南京大学地球科学与工程学院， 南京 210046， 中国)

0 引 言

孔隙网络是流体在土体内部赋存的物理空间，它决定着土体的渗透性，控制着土体内部流体及物质的迁移速率，并使土体在变形阶段产生超孔压及排水固结行为； 此外，砂土液化、非饱和土的毛细管力及基质吸力等重要的工程性质本质上也都由土体孔隙网络所控制(Mitchell et al.，2005)。

分析扫描电镜(SEM)数字图像可获得土样的孔隙结构定量特征。施斌等采用莫斯科大学研制的Videolab图像分析系统，获得了郧县膨胀黏土、太平洋洋底软泥等样品SEM图像的孔隙数量、孔隙相对面积、平均孔径、孔径分布以及孔隙平均形状系数等形貌特征与定向性指标(施斌等， 1995； 施斌， 1996)。Liu et al. (2011)基于数字图像理论，自行编制了孔隙与裂隙分析系统(PCAS)，并用PCAS分析了8张SEM图像，研究黏土样在直剪试验过程中孔隙结构变化的定量特征，所获得信息包括孔隙面积、孔隙数、平均形状系数、孔隙周长、定向性以及概率熵、概率分布指数与分形维数等。但需指出，分析SEM图像仅能对土样在平面内的孔隙结构加以表征。压汞(MIP)与氮气吸附法(NA)能获得土样的孔径分布与比表面积(李志清等， 2017)，却难以进一步表征土样孔隙及其网络的定量几何特征。与以上方法相比，工业显微CT(μCT)具有无损、非侵入式透视成像的优点，可提供微米尺度的高分辨率二维层析图像(刘治清等， 2017)。张巍等(2017)使用南京粉砂土样，基于三维二值矩阵运算，对CT成像表征单元体内沿空间任意方向切面的孔隙定向性进行了表征，同时还获得了沿土样高度方向孔隙率的空间分布。尽管使用显微CT能够获得土样的三维重构模型，但要从中有效提取出孔径大小、分布、连通性、形状因子等表征孔隙网络的空间几何特性参数，仍需进一步建模。

近年来，在页岩油、气开发与二氧化碳地质封存等能源、环境工业背景驱动下，对岩石孔隙尺度成像与建模(Pore-scale imaging and modelling)的技术与方法得以不断深入(Blunt et al.，2013； 赵斌等， 2018)。最大球算法(Principle maximal ball algorithm)是一种多孔介质内部孔隙网络的建模方法，该方法建立在工业显微CT扫描成像的基础之上，Fatt(1956)首先提出了孔隙网络模型的概念，用于表征多孔介质的毛细管压力。在此基础上，Silin et al. (2006)提出了最大球算法，采用球-杆模型表征孔隙网络，其核心思想是，将多孔介质的孔隙空间视作由孔隙体与孔喉所组成的几何体，孔隙体为流体所赋存较大开放孔隙空间，孔隙体决定着多孔介质的孔隙率； 而孔喉则是连接孔隙体之间相对狭长的孔隙通道，可视作供流体迁移的毛细管，孔喉决定着多孔介质的渗透率； 配位数为单个孔隙体与其他孔隙体相联接的孔喉数量。Al-Kharusi et al. (2007)发展了最大球算法，加入了最大球簇的定义，对孔隙网络的分析结果更加精细化，但该算法限制最大孔隙数量不能超过1000个，实用性不强。

Dong et al. (2009)改进了最大球算法，打破了之前的主球、从属球和最大球簇的系统分层结构(Fatt， 1956)，将最大球仅区分为孔隙和喉道，孔隙与喉道相连组成孔隙网络，开发了相应的算法，可统计分析多孔介质的孔隙网络的拓扑结构； 并基于最大球算法，提取出标准砂岩的孔隙网络参数，与实验数据比较结果表明，计算与实验结果一致性程度较高，证明了最大球算法的可靠性。此外，王合明(2013)利用最大球算法提取了4种不同的填砂模型的孔隙网络模型，计算所得配位数与理论配位数值十分接近，进一步验证了最大球算法提取砂土孔隙网络模型的可靠性。目前，最大球算法已被广泛地应用于岩石孔隙网络提取及岩石内部单相流与多相流模拟(Bijeljic et al.，2013； Raeini et al.，2014)，但应用于土体孔隙网络提取的研究刚刚起步。肖瑞等(2017)采用自主研发的土体细微观三轴压缩在线观测控制系统，得到了南京粉细砂饱和试样在不同围压下加载过程中孔隙结构的实时扫描图像并提取出孔隙网络模型，采用孔隙率、孔隙形状等指标分析了孔隙结构的演化过程，研究了孔隙结构的演化规律与其力学性质和强度特征之间的关系。刘语等(2019)从南京粉细砂三维重构模型中提取出5组立方体孔隙表征单元体(REV)，采用最大球算法对每个REV建立孔隙网络模型，从中提取孔隙率、单位体积孔隙数、孔隙平均体积等8个孔隙结构参数，利用假设检验确定了表征样品孔隙结构参数的统一REV边长尺寸。曹剑秋等(2018)获取了南京粉砂试样在三轴加载试验过程中不同应力状态下的微米级孔隙结构图像，基于最大球算法，构建出空间孔隙网络模型并统计了土体的细观孔隙参数，分析了三轴加载过程中细观孔隙结构演化特征与其宏观力学性质之间的相关性。许林等(2018)对砂雨法成样与振捣干法成样的南京粉砂、鄂尔多斯油砂、自贡砂岩和砂雨法成样玻璃珠等试样依次进行了显微CT扫描成像，基于最大球算法，提取出16种空间孔隙REV细观孔隙结构参数，并分析了各参数对多孔介质渗透率的影响。

相较于MIP与NA等实验方法，基于最大球算法建立土体空间孔隙网络模型，不仅能获得土样的孔径分布，还能对孔隙连通性、形状因子等进行定量表征； 不同于数字图像处理技术只能获取平面内的孔隙分布特征，最大球算法提取的三维孔隙网络模型能够更加全面地刻画出砂土空间孔隙结构特征。本文首先介绍了最大球算法的基本原理，以振捣干法生成南京粉砂试样，通过三维二值矩阵运算，获得了显微CT扫描试样的空间孔隙REV； 并基于最大球算法，建立了空间孔隙REV相应的孔隙网络模型，对空间孔隙的配位数、孔隙半径，喉道内切圆半径、喉道长度等一系列孔隙网络拓扑参数进行了统计分析。

1 最大球算法

1.1 基本原理

以下简述Dong et al.(2009)改进后最大球算法的基本原理。

最大球算法即利用若干个球体占据多孔介质的孔隙空间，以达到区分孔隙和喉道的目的。若某个球体完全处于某个孔隙空间中且不被其他的球完全包含，则该球可以被称作一个最大球，即一个局部的最大球被定义为1个孔隙体，搜索最大球之间的链接即得到喉道。

最大球与骨架概念如下所示：将二值化切面图片重构为三维数据矩阵后，对于三维空间中的每个点的坐标，均可用向量p=(x，y，z)的形式来表示，则任意两点p1、p2间距离Dist(p1，p2)为：

(1)

=(xg-x0)2+(yg-y0)2+(zg-z0)2

(2)

p0

(3)

为表征不同大小的孔隙空间的结构特征，采用最大球簇来对不同大小的最大球进行分类。可分为最大球单簇和最大球多簇。其中最大球单簇为一个最大球与所有半径小于它并且与它重合的最大球的集合(图2)。

图1 离散体素所构成的不同半径的最大球(Al-Kharusi et al.，2007)Fig. 1 Different radius of the maximal ball with discrete voxels

图2 最大球单簇算法示意图(Al-Kharusi et al.，2007)Fig. 2 The illustration of the maximal cluster algorithm

生成最大球单簇的方法是：首先将最大球的两倍半径设定为搜索范围，然后在上述范围内搜寻半径小于该球且与之相切或相交的所有最大球，即：

Dist(p1，p2)

(4)

再将搜寻到的最大球归入最大球单簇。

每个最大球都可以重叠与它相邻的半径小于它的最大球，而这些被重叠的最大球又可以再这样去重叠其他的最大球。这样的过程最终生成最大球多簇(图3)。而在该多簇中，半径最大的球被称作为该簇的祖先，对应地，每个祖先都代表了1个孔隙空间，而当一个最大球有两个祖先时，该球就对应了一个喉道，喉道最大球到孔隙最大球之间的多簇被定义为孔隙喉道链。

最大球算法是利用若干个球体去占据孔隙空间，而孔隙空间源自于通过计算机断层扫描成像(CT)等技术获得的数字图像，多张连续数字图像的像素点构成了孔隙空间的体素点，体素点是孔隙空间的基本单元。因为数字图像的质量受图像分辨率大小的影响，所以在基于体素点进行最大球搜索以及孔隙与喉道的识别划分时，更高分辨率成像可以识别出更精细的孔隙与喉道特征。

1.2 孔隙网络模型

孔隙网络模型即将多孔介质孔隙空间等效为一系列通过喉道相连的孔隙体，孔隙体被简化为圆球，喉道被简化为细杆，孔隙空间被简化为一系列球-杆相连的网络状结构(Silin et al.， 2006)，所构建的土体孔隙空间中，输入的最大球的数据是一个三维二值矩阵，其中1表示骨架， 0表示孔隙。孔隙网络模型构建具体方法如下：

首先对土体样品的灰度图像进行二值化，随后将二值图像序列表征为三维二值矩阵，据此可对土样进行三维重构并提取孔隙REV用于空间孔隙网络模型建模。孔隙REV中，每个体素点都对应一个最大球，在孔隙空间中任选一体素点，向其周围的26个方向搜寻与其接触的颗粒体，即骨架体素。随后，依据该体素的形成范围，采用收缩算法，寻找最大球以及其半径的上下限，确定该最大球的大小和位置，再根据上述方法继续寻找下一个最大球。如此循环往复，得到所有最大球。

在以上操作过程中，可能有冗余球出现，冗余球是指完全被包含在另一个最大球中的最大球，这类球对于孔隙体的表征无效的。假定搜索到的球B被完全包含在球A中，则满足：

Dist(pA，pB)<|RLA+RLB|

(5)

式中，RLA和RLB分别为两球的半径下限，将上式中的RLB换成RUB，RUB为B球的半径上限，则可删除冗余球B。

当所有冗余球清除完毕后，孔隙REV中所有体素都已经转化为最大球，下步需对所有最大球进行孔隙或喉道的分类。首先按半径上限的递减顺序，将最大球分为相同半径的若干组。假设第1组最大球个数为N，最大球初始分组数量默认为无穷大。首先选取半径最大的一组最大球C，定义C为最大球多簇的祖先，设其级数为1并且包含所有与其相邻或是重叠的最大球，这些被包含的最大球级数设为2。按递减顺序，对第1组中余下的N-1个最大球排序，取出级别最高的最大球D，球D包含所有半径小于它的最大球，这些最大球未被赋予级数，此时赋予它们的级数比D多1。若排序前D为无穷大，则D作为一个祖先定义了一个新的孔隙，反之则D将自身祖先的编号往下一级传达。若D中包含有另一簇中的最大球，那么D就是一个公共最大球，定义其为喉道。喉道被识别之后，从喉道到每一簇的祖先会被分别定义为两个孔喉链。对最大球使用该算法进行成簇与排序，直到到达事先设定的最小半径，该情况下若按原先设定，所有最大球均会被识别为喉道。因此，将(最小可搜寻的)只含有7个体素点的最大球设定为最小孔隙，将(最小可识别)只含一个体素点的最大球设为最小喉道(图4)。

图4 最大球算法最小孔隙及喉道Fig. 4 Minimum pore and throat of maximal-balls algorithma. 最小孔隙； b. 最小喉道

区分孔隙与喉道最大球后，采用球体表示孔隙(图5a)； 采用与最大球等直径的圆柱体表示喉道(图5b)； 采用可视化软件即可输出所得到的孔隙网络模型(图5c)。

图5 孔隙网络球-棍模型Fig. 5 Pore network ball-stick modela. 孔隙最大球； b. 喉道圆柱体； c. 孔隙网络模型(球棍模型)

下文以南京粉砂为例，详述土样空间孔隙网络建模及其特征参数的提取步骤，并对所生成的孔隙网络参数进行统计分析。

2 材料与方法

2.1 试样制备与扫描成像

本试验所用砂土试样为南京河西地区的南京粉砂(图6a)，物理、力学性质及粒径分布特征参见(张巍等， 2017； 刘语等， 2019)。为保证土颗粒粒径的均匀程度，消除黏土颗粒和破碎颗粒的影响，在配置土样时只保留0.075～2imm的粉细砂，采用振捣干法配制含水量为0的重塑干样。

采用显微CT技术进行砂样的孔隙尺度成像，扫描工作在中国科学院南京土壤研究所完成。扫描设备及技术规格参见(张巍等， 2017)，使用14iμm分辨率，采用内径8imm、高度40imm的透明密封离心管装样。CT扫描工作电压为110ikV，工作电流为110iμA。样品扫描时间约为60imin，共获得2000张断层序列扫描的栅格数据(TIFF)格式图像(8位灰度)。

2.2 图像处理与三维重构

采用可视化图像处理软件ImageJ完成图像增强及三维重建，得到土样三维重构模型(图6b)。

图6 三维重构模型及表征单元体Fig. 6 3D reconstruction model and the representative element volumesa. 砂土试样； b. 三维重构模型； c. 表征单元体

受计算机硬件处理能力的限制，孔隙尺度成像与建模在REV尺度上进行。刘语等(2019)研究确定了适用于南京粉砂样品孔隙结构参数的统一REV边长尺寸为400像素(2.6imm)。由于本文采用与以上文献相同的砂土与显微CT成像分辨率，因此直接沿用以上结论，分别在样品三维重构模型的上、中、下3个位置中心各选取1个400×400×400iREV，得到图6c所示的3个颗粒REV。

3 孔隙网络

以图7a中颗粒表征单元体REV1为例，通过三维二值反运算得到其对应的孔隙表征单元体REV1′(图7b)，基于最大球算法提取出对应的孔隙网络模型N-REV1′(图7c)。孔隙网络模型中，球体表示孔隙，圆柱体表示喉道，球体或者圆柱体半径越大，表示孔隙或者喉道半径越大。

图7 孔隙网络模型生成Fig. 7 Pore-network model generationa. REV1； b. REV1′； c. N-REV1′

从N-REV1′中可提取出多个表征参数，包括表征孔隙网络连通性的配位数、表征孔隙拓扑结构的孔喉内切圆半径、体积、长度等以及表征孔隙形状的形状因子等。

4 参数统计分析

对所取3个REV的孔隙网络模型参数数据进行统计，并统计分布特征值，绘制出各参数频率分布图(图8～图12)。

图8 REV孔隙半径、喉道半径分布频率图Fig. 8 Distribution frequency diagram of pore and throat radiusa. REV孔隙半径分布频率图； b. REV喉道半径分布频率图

图9 REV配位数分布频率图Fig. 9 Distribution frequency diagram of coordination number

在孔隙空间中，与单个独立的孔隙体相连的所有喉道的数目称为该孔隙体的配位数。配位数直接表征孔隙网络的连通性好坏，而连通性会直接对介质的输运特性造成影响。一般情况下，孔隙的配位数在1～15之间，配位数越大，孔隙的连通程度就越高。图9为所取3个REV平均配位数分布频率图，由图可知，试样配位数大部分分布在0～25区间内，在7～9区间附近出现频率分布峰值，约为8%。对分布频率进行高斯曲线拟合，相关系数为0.94，表明配位数分布近似服从正态分布。由拟合的正态分布曲线可知：配位数数学期望10.0，标准差5.1，变异系数0.51，表明配位数分布较集中，孔隙连通性较强。

图10 REV孔隙截面积形状因子、喉道截面积形状因子分布频率图Fig. 10 Distribution frequency diagram of area shape factor of pore and throata. REV孔隙截面积形状因子分布频率图； b. REV喉道截面积形状因子分布频率图

细粒土体的孔隙或喉道形状相对复杂，为了采用规则的几何形状表示其形状，引入形状因子定量表征，其取值越小说明形状越不规则。图10为孔隙截面积形状因子、喉道截面积形状因子分布频率图。对孔隙和喉道截面积形状因子分布频率分别进行高斯曲线拟合，得到相关系数为0.961和0.990，表明孔隙和喉道截面积形状因子均近似服从正态分布。孔隙与喉道的截面积形状因子分别分布在0.01～0.04与0.01～0.05区间内，其正态曲线特征为，孔隙截面积形状因子的数学期望值0.019，喉道截面积形状因子的期望值0.033，均小于正三角形的形状因子标准值0.048，表明试样绝大部分孔隙和喉道为不规则的几何形态，且孔隙几何形态较喉道更为复杂。孔隙截面积形状因子的标准差0.005i54，变异系数Cv=0.291； 喉道截面积形状因子的标准差0.006i79，变异系数Cv=0.206，表明孔隙和喉道的截面积形状因子分布较集中，孔隙截面积因子较喉道分布更为分散。

图11 REV喉道长度分布频率图Fig. 11 Distribution frequency diagram of throat length

图12 REV孔隙体积分布频率图Fig. 12 Distribution frequency diagram of the REV pore volume

5 结 论

(1)最大球算法可用于建立多孔介质空间孔隙网络球棍模型，该算法使用局部的最大球定义孔隙体，通过搜索最大球之间的链接获得喉道。

(2)对显微CT二值图像序列进行三维重构，从三维重构模型中可提取出任意表征单元体REV，通过三维二值反运算可获得其所对应的孔隙表征单元体REV′，基于最大球算法可获得REV′所对应的孔隙网络模型N-REV′。

(3)计算获得了样品REV孔隙网络模型参数，统计分布特征后发现，孔隙半径、喉道半径、孔隙配位数、孔隙截面形状因子、喉道截面形状因子与喉道长度等孔隙参数均近似服从正态分布，孔隙体积近似服从衰减型指数分布。

(4)孔隙与喉道半径分别分布在100iμm与65iμm以内，两者数学期望分别为40.0iμm与18.0iμm； 配位数分布在25以内，数学期望为5.1； 孔隙与喉道截面形状因子分别分布在0.01～0.04与0.01～0.05的区间内，两者数学期望分别为0.019与0.033； 喉道长度分布在100～800iμm以内，数学期望为292.22iμm。

(5)样品中体积小于1.5×107iμm3的小孔隙数量超过90%。