砾石土特征分析及数字岩心重构

Analysis of Gravelly Soil Characteristics and Digital Core Reconstruction

摘要: 砾石土中的孔隙结构对土体的渗透性和稳定性有很大影响，砾石土微观结构特征直接影响其宏观性能，为了更深刻认识砾石土内部细观结构，为砾石土中的管涌、流土等渗透侵蚀破坏防治提供参考。对于砾石土制取样品后进行CT扫描获取图像，对图像进行预处理，基于图像重建出三维数字岩心。通过其孔隙度、孔隙连通性、孔径和粒径等对砾石土的结构特征进行分析，从而得到砾石土的细观结构特征。选取砾石土体孔隙充分发育的一张二维图像，采用模拟退火法进行重构，以两点相关函数、线性路径函数作为约束函数。成功重构出与真实土体孔隙结构特征相似的数字岩心。模拟退火法建模所需建模资料要求较低，是一种经济、高效的建模方法。采用两点相关函数和线性路径函数对于重构图像进行评价，可以得到模拟退火算法重构的砾石土孔隙结构与真实土体细观结构较为相似，重构出的砾石土数字岩心效果较好。

Abstract: The pore structure of gravelly soil has a great influence on the permeability and stability of soil mass, and the microstructural characteristics of gravelly soil directly affect its macroscopic properties, which provides a reference for preventing and controlling the seepage erosion damage of piping and flowing soil in gravelly soil. For the gravelly soil samples after the acquisition of a CT scan image, image pre-processing is based on the image reconstruction of a three-dimensional digital core. The microstructure characteristics of gravelly soil are analyzed by porosity, pore connectivity, pore size and particle size. A two-dimensional image of gravelly soil with well-developed pores was reconstructed by simulated annealing, and the two-point correlation function and linear path function were taken as constraint functions. The digital core is reconstructed successfully, which is similar to the pore structure of real soil. Simulated annealing (SA) is an economical and efficient modeling method, which requires less modeling data. By using the two-point correlation function and linear path function to evaluate the reconstructed image, the pore structure of gravel soil reconstructed by simulated annealing algorithm is similar to that of real soil, and the reconstructed gravel soil digital core has a better effect.

1. 引言

砾石土是多孔介质结构，主要有孔隙和土骨架两部分组成。内部孔隙结构相互交错、错综复杂形成孔隙结构网。土体的孔隙分布和颗粒分布对于其渗透性和稳定有很大影响 [1] 。孔隙是指土壤中颗粒之间的空隙，它们可以影响土体的透水性、通气性和保水性等特性。砾石土中的孔径分布直接影响着土体的渗透能力和排水能力。当孔隙较大时，土体渗透性较好；而当孔隙较小时，土体的保水性较强。粒径分布描述的是砾石土中颗粒大小的分布情况，不同粒径的颗粒对土体的物理性质产生重要影响。通过粒径分布，可以评估砾石土的颗粒组成及比例，进而确定土体的工程性质。颗粒和孔隙几何信息、尺寸分布、形状、体积和表面积都会对渗透性产生影响 [2] 。砾石间接触次数的增加，砾石土的抗液化能力随砾石含量的增加而显著增加 [3] 。对于砾石土的研究，Yusuke Fujikura等提出从细砂到粗砂砾的孔径分布模型和两种基于孔径分布估算渗透率的方法 [4] 。Qian Zhai等通过孔径分布函数来估算砂土的渗透性 [5] 。金爱兵等对不同粒径分布的尾砂进行了研究，探讨了其充填体的强度和损伤特性。细尾砂的充填体中，粗颗粒与细颗粒相互连接紧密，并被水化产物完全包裹，使得孔隙度较低，从而形成了稳定的骨架结构。此外，细尾砂的微观结构也表现出较好的致密性 [6] 。

通过统计相关函数可以来定量表征砾石土微观结构。对多孔介质和复合材料等随机非均质材料的有效输运、电磁和力学性能的基本理解，依赖于定量表征介质微观结构的能力 [7] 。事实上，微观结构的完整表征需要无限一组统计相关函数的知识。在实践中，无论从实验上和理论上都只能获得低阶信息。确定是否存在同时包含重要微观结构信息的低阶函数仍然是一个重要的研究方向。B Lu等提出了线性路径函数，并给出了函数表达式。两相异质材料的一个基本形态测量是所称的线性路径函数L(z)。这个量给出了一个长度为z的线段完全在其中相中的概率 [8] 。A. Derossi等人认为食品微观结构与其质量之间的关系的重要性被认为是至关重要的。在面包质量方面，面包屑结构的精确表征被认为是评价其感官特性的必要条件，线性路径函数用于统计面包的内部结构特征 [9] 。孔隙率、两点相关函数线性路径函数可以反映多孔介质的基本特征。通过统计相关函数来定量表征多孔介质的细观结构特征，可为土体的重构提供控制函数。

传统的岩心分析需要进行物理切割、处理和测量，而数字岩心可以通过CT扫描技术快速获取高精度的样本数据，进行可视化和重构。目前重构数字岩心方法有截断高斯模拟法 [10] 、蒙特卡洛法 [11] 、多点地质统计法 [12] 和模拟退火法等 [13] 。而模拟退火法所需建模资料要求较低，数字岩心重构较为经济、高效。借助于CT扫描出土体的二维平面图像，对图像进行去噪、二值化等处理。然后就可以提取重构所需要的孔隙率、两点相关函数和线性路径函数，利用模拟退火法可重构出土体的数字岩心，就可以得到与原始CT图像较为相似的细观孔隙结构。使用孔隙率、两点相关函数和线性路径函数对于重构数字岩心的优劣进行评价。

2. 图像处理

将砾石土样放入CT进行扫描，获得1014张1004 × 1024像素的CT二维灰度切片图。在Image J软件中调整其对比度，在孔隙发育丰富的地方批量截取700 × 700像素大小图像共700张。在图像成像和传输过程中，会受到设备与周围环境，总会或多或少的带来一些不能避免的噪声。对图像进行去噪处理，采用中值滤波法对图像进行去噪 [14] ，用MATLAB软件中的medfilt2函数，滤波器大小为3 × 3。最后在Avizo软件中将700张图像进行叠加，重建成为三维图像，如图1所示。接下来就可以对三维图像进行细观结构特征提取。

Figure 1. Reconstruct a three-dimensional image

图1. 重建三维图像

3. 细观结构分析

3.1. 孔隙度

孔隙度是一个非常重要的土体物理参数，它决定了土体中水分、空气的存储和运移能力。孔隙度计算表达式见下式(5)：

$n=\frac{V_v}{V_{总}}$ (1)

式中： $V_v$ 表示孔隙的体积； $V_{总}$ 表示总的体积。

在Avizo软件中利用Interactive Threshold模块将二维图像进行阈值分割，分为孔隙和土体颗粒两部分。二维图像重建成三维图像的孔隙分布图如图2所示，用Volume fraction模块可计出体孔隙度为0.214，逐层面孔隙度大小分布图如图3所示。可以看出该砾石土土样是非均质体，最大孔隙度为0.391，最小孔隙度为0.085，平均孔隙度为0.214。

Figure 2. Distribution of pores

图2. 孔隙分布图

Figure 3. Porosity layer by layer

图3. 逐层面孔隙度

3.2. 孔隙连通性

砾石土孔隙连通率是指砾石土中孔隙之间的连通程度。如果孔隙之间相互连接，并且形成连续的通道，水分能够更容易地渗透和流动。相反，如果孔隙之间存在较多的隔离或者孤立，水分渗透性将受到限制。通过孔隙连通性分析，了解到孔隙的连通性情况，进而评估土体的渗透性。在Avizo软件中先利用Interactive Threshold模块选择孔隙部分，然后用Axis Connectivity模块做砾石土的连通性分析。连通孔隙分布如图4所示。

砾石土中的孤立孔隙是指孔隙与其他孔隙之间没有有效的连通通道，无法与周围的孔隙或介质进行交换或传递物质。砾石土中的孤立孔隙对于水分和气体的流动具有一定的限制作用。由于孔隙之间缺乏连通通道，水分在孤立孔隙中可能滞留较长时间，难以通过砾石土层快速渗透和排水。孤立孔隙影响着砾石土的渗透性。在Avizo软件中利用Arithmetic模块分离出孤立孔隙分布如图5所示。

有效孔隙就等于全部孔隙去除孤立孔隙。通过Avizo软件中Volume fraction模块计算，可得有效孔隙度为0.210，连通孔隙占比97.9%，孤立孔隙占比2.1%。有效逐层面孔隙度大小分布图如图6所示，最大有效孔隙度为0.391，最小有效孔隙度为0.075，平均有效孔隙度为0.210。

Figure 4. Connected pore distribution diagram

图4. 连通孔隙分布图

Figure 5. Distribution of isolated pores

图5. 孤立孔隙分布图

Figure 6. Effective layer by layer porosity distribution

图6. 有效逐层面孔隙度分布图

3.3. 孔隙定量分析

砾石土中的孔隙结构对土体的渗透性有很大影响，孔隙定量分析可以了解土体中不同尺寸孔隙的存在情况 [15] 。较大的孔隙能够容纳更多的水分，从而增加土体的渗透能力；而较小的孔隙则会限制水分渗透，降低渗透性。

Figure 7. Schematic diagram of pore segmentation

图7. 孔隙分割示意图

Figure 8. Aperture size distribution

图8. 孔径大小分布图

首先，在Avizo软件中使用Interactive Threshold模块筛选出孔隙部分，然后利用Separate Objects模块将孔隙分离。经过分离后，各孔隙分布如图7所示。最后，通过使用Label Analysis模块，对各个孔隙的等效孔径、体积和表面积进行了统计。可以得到孔隙结构特征参数，用来定量对孔隙结构特征进行描述。孔隙的最大孔径为4924.10 μm，平均孔径为93.27 μm，筛选后孔隙的孔径大小分布图如图8所示。孔隙中最大体积为2.01 × 10¹⁰ μm³，平均体积为2.62 × 10⁷ μm³，筛选后孔隙体积大小分布图如图9所示。孔隙中最大表面积为4.67 × 10⁸ μm²，平均表面积为4.12 × 10⁵ μm²，筛选后孔隙的表面积大小分布图如图10所示。

Figure 9. Pore volume distribution diagram

图9. 孔隙体积分布图

Figure 10. Pore surface area distribution diagram

图10. 孔隙表面积分布图

3.4. 颗粒定量分析

砾石土的颗粒大小分布对于其稳定性和渗透性有很大影响。若细颗粒分量较多，粗颗粒“浮于”细颗粒中易导致砾石土产生管涌甚至流土破坏；若粗颗粒分量较多，细颗粒会在粗颗粒孔隙中被过滤移动带走，但粗颗粒形成相互支撑的骨架结构而使砾石土整体保持稳定。粗颗粒含量的增加，砾石土渗透系数也随之增大 [16] 。定量分析砾石土颗粒可以了解不同尺寸颗粒的存在情况。

使用Interactive Threshold模块选择颗粒部分，然后使用Separate Objects模块将颗粒进行分离，分离后各颗粒分布如图11所示，最后使用Label Analysis模块对各个颗粒的等效粒径大小、体积、表面积分别进行统计，用来定量的对颗粒特征进行描述。颗粒最大粒径为5114.36 μm，平均粒径为1039.81 μm，筛选后的粒径分布图如图12所示。颗粒最大体积为1.06 × 10¹¹ μm³，平均体积大小为2.01 × 10¹¹ μm³，筛选后的颗粒体积分布图如图13所示。颗粒最大表面积为2.51 × 10⁸ μm³，平均体积大小为1.16 × 10⁷ μm³，筛选后的颗粒体积分布图如图14所示。

Figure 11. Particle segmentation schematic diagram

图11. 颗粒分割示意图

Figure 12. Particle size distribution diagram

图12. 粒径大小分布图

Figure 13. Particle volume distribution diagram

图13. 颗粒体积分布图

Figure 14. Surface area distribution of particles

图14. 颗粒表面积分布图

4. 统计相关函数

统计相关函数可以准确表示微观孔隙的结构特征，可以用来重构砾石土的数字岩心。选取砾石土体孔隙充分发育的一张二维图像，截取150 × 150像素大小为原始土样图像，用于重构数字岩心，如图15所示。提取其两点相关函数和线性路径函数，为重构数字岩心过程中的控制函数。

Figure 15. Primitive soil samples

图15. 原始土样

4.1. 两点相关函数

两点相关函数定义为在介质中任意两点，这两点分布于其中一相位(孔隙或土体颗粒)中的概率 [17] 。砾石土任意一相的两点相关函数定义为：

$S^{\left(i\right)}\left(r_1,r_2\right)=\langle Z^i\left(r_1\right)Z^i\left(r_2\right)\rangle$ (2)

式中：S表示两点相关函数；r₁、r₂表示在某一项任取两个点；< >表示取平均。

根据公式可以知道，当r = 0时， $S\left(r\right)$ 就等于孔隙度，r表示像素数。原始土样图15的孔隙度为0.1509。从孔隙相两点相关函数分布图16可以得到，当任意两个点的距离超过15像素后，两点相关函数的值会趋于稳定。对目标图像重构时，选取最大的r为61像素。

Figure 16. Two-point correlation function

图16. 两点相关函数

4.2. 线性路径函数

线路路径函数 $L\left(r\right)$ 表示在某一相位找到一个长度为r的线段的概率。线性路径函数可以表征相与相之间的连通性，因此它可以区分不同的相。在砾石土中，可用于描述其连通性 [18] 。其表达式为：

$L^{\left(i\right)}\left(r\right)=\langle P\left(r_1,r_2\right)\rangle$ (3)

式中：L表示线性路径函数；r表示任意长的线段；< >表示取平均。

从孔隙相线性路径函数图17分析可得，当任意两点的长度超过8个像素时，线性路径函数的值就会趋近于0。

Figure 17. Line-path function

图17. 线性路径函数

5. 模拟退火法重构数字岩心

模拟退火算法(SA)来源于固体退火原理，是一种基于概率的算法。将固体加温至充分高的温度，再让其慢慢冷却，加温时，固体内部粒子随温升变为无序状，内能增大，分子和原子越不稳定。而冷却时粒子渐趋有序，能量减少，原子越稳定。在降温过程中，固体在每个温度都达到平衡态，最后在常温时达到基态，内能减为最小。

5.1. 算法流程

模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，即在局部最优解能概率性地跳出并最终趋于全局最优。模拟退火算法是通过赋予搜索过程一种时变且最终趋于零的概率突跳性，从而可有效避免陷入局部极小并最终趋于全局最优的结构优化算法。

是否接受新解，是根据Metropolis准则按一定概率接受。在某一温度T下，当前状态能量为 $E\left(k\right)$ ，新状态能量为 $E\left(k+1\right)$ ， $\Delta E=E\left(k+1\right)-E\left(k\right)$ ，接受k + 1状态的概率p为：

$p=\{\begin{array}{cc}1,& \Delta E<0\\ \exp \left[-\frac{\Delta E}{T}\right],& \Delta E>0\end{array}$ (4)

当∆E < 0时，接受k + 1为当前状态；当∆E > 0时，则在0~1中取一个随机数，若接受概率 $P=\exp \left[-\Delta E/T\right]$ 于这个随机数，则接受新状态k + 1为当前状态，保留状态k为当前状态。在每个温度下进行多次系统状态的改变，使系统在这个温度下达到较为稳定的状态，然后再降低温度进入下一个循环，继续进行系统状态的优化。

5.2. 数重构数字岩心

首先随机的生成150*150像素且与原始土样图像孔隙度相同的数字岩心，1代表白色孔隙，0代表黑色土体颗粒。然后开始进行模拟退火算法。

引入更多的反应砾石土细观孔隙结构的基本特征，重构的砾石土的细观孔隙结构才会更加接近真实结构。采用两点相关函数和线性路径函数为约束函数，分别计算重建系统中的、 $L\left(r\right)$ ，在根据下式(5)计算当前时刻的系统能量E。在随机的数字岩心中任意选取一孔隙点1和一土体颗粒点0交换其位置得到新的系统。在生成图像的孔隙率始终与目标图像保持一致。然后计算新系统的 $S\left(r\right)$ 、 $L\left(r\right)$ 和系统能量E。根据Metropolis准则是否接受新系统。

$E={\displaystyle \sum {\alpha }_i\cdot {\left[S_r\left(r_i\right)-S\left(r_i\right)\right]}^2}+{\beta }_i\cdot {\left[L_r\left(r_i\right)-L\left(r_i\right)\right]}^2$ (5)

$S\left(r_i\right)$ 是目标图像的两点相关函数， $S_r\left(r_i\right)$ 是对应于 $S\left(r_i\right)$ 两点相关函数。同理 $L_r\left(r_i\right)$ 是 $L\left(r_i\right)$ 对应的线性路径函数，建模的实质就是对系统不断进行迭代，直到两点相关函数和线性路径函数足够逼近时，当前的图像就是重构出的最优图像。

模拟退火的参数设置为：初始温度T0 = 1，终止温度Te = 10⁻⁸，降温系数alpha = 0.99，同一温度下演化次数N = 100，计算中土体孔隙两点相关函数和线性路径函数的最大距离r_max分别为61、17，不同距离r处孔隙两点自相关函数对应的权重值当r ≤ 25时，α_i = 1.5，β_i = 1.5，否则α_i = 1，β_i = 1。不同距离r处孔隙线性路径函数对应的权重值当r ≤ 11时，α_i = 1.5，β_i = 1.5，否则α_i = 1，β_i = 1。

5.3. 重构结果分析

模拟退火法重构的数字岩心如图18所示，经过与原始土样图15的对比可以观察到重构数字岩心图像的细观孔隙结构与原始土样图像较为相似。图19为孔隙相重构数字岩心和原始土体两点相关函数和线性路径函数的结果分布图，可以观察到两点相关函数计算结果接近，表示原始土样图像和重构数字岩心的孔隙空间分布状态相似；线性路径函数计算结果接近，表示原始土样图像和重构数字岩心的孔隙大小和连通状态相似。模拟退火法重构效果较好。

Figure 18. Reconstruct the soil structure

图18. 重构土体结构

Figure 19. Reconstruction of digital core and original soil function contrast map (a) Two-point correlation function (b) Line-path function

图19. 重构数字岩心和原始土体函数对比图 (a) 两点相关函数 (b) 线性路径函数

6. 结论

基于砾石土的CT扫描图像，重建出三维特征图像，并对三维图像进行特征分析，并选取一张二维图像采用模拟退火法对土体进行数字岩心重构，得出以下结论：

1) 从孔隙度、孔隙连通性、孔径分布和粒径分布对砾石土的结构进行分析，定量的表示砾石土的细观结构特征，可以得到砾石土的孔隙分布和颗粒分布的特征。孔隙的最大孔径为4924.10 μm，平均孔径为93.27 μm，颗粒最大粒径为5114.36 μm，平均粒径为1039.81 μm。

2) 计算二维图像的两点相关函数和线性路径函数作为约束函数，采用模拟退火法成功重构出数字岩心，对比其两点相关函数和线性路径函数，重构出的数字岩心与原始土样结构特征相似，可为砾石土的进一步研究提供基础资料支持。