基于工业CT 扫描的不同类型土壤孔隙结构研究

阚晓晴 程金花 王 葆 郑 欣

（1.北京林业大学水土保持学院，北京 100083；2.国家林业局水土保持与荒漠化防治重点实验室，北京 100083；

3.北京林业大学水土保持学院建水研究站，北京 100083）

土壤结构是由生物和非生物因素（包括气候、矿物成分、有机质、根、真菌菌丝、土壤动物和耕作）相互作用，形成的固体土壤成分和孔隙在不同尺度上的三维排列[1]。其中，大孔隙是出现在土壤中大体积、连续的结构性开口，为流经土壤的物质提供快速移动路径[2]。大孔隙在土壤中普遍存在于植被生长的非饱和带[3-5]，分布相对密集[6]，包含土壤裂隙、根系通道、生物洞穴等，土壤中的大孔隙长度、直径、体积、表面积、网络结构等发育特征决定了土壤的连通性，进而影响水分运动[7-8]及溶质运移[9]，因此土壤大孔隙特性对水资源[10]、水分运移[11]、分子物质[6]和污染物[12]在土壤中的迁移及滞留等问题的研究至关重要。土壤大孔隙特性及其与水分运移的关系是土壤水分运动的难点和前沿问题之一。陈晓冰[13]和姚晶晶[14]的研究结果表明，重庆红壤区优先流现象明显，与土壤总孔隙度相关性显著；王葆[15]和郑欣[16]的研究结果表明，北京褐土区土壤大孔隙形态特征显著，水分优先运移现象与土壤总孔隙度呈正相关；Kan 等[17]、彭旭东[18]和朱晓锋[19]的研究结果表明，喀斯特红壤区水分入渗机制受可溶性岩石影响，水分沿土–石裂隙优先下渗现象显著。因此，不同土壤类型条件下的土壤孔隙量分别对当地土壤水分运移影响显著，然而土壤结构的异质性会导致土壤水力性质的变化[20]，孔隙尺度的非均质性结构也会控制着土壤的水分运移过程[21]，但水分运移及滞留严重制约着植被恢复[22]，因此在不同土壤类型条件下的土壤孔隙结构特征亟待研究。

土壤大孔隙的研究方法多为染色示踪法[23]、穿透曲线法[24]等，但因其破坏性较强，无法对土壤结构进行连续性观测。Jarvis 等[25]发现X 射线的CT 成像可表示土壤中的孔隙；Soto 等[26]发现了CT 扫描获得的孔隙路径在流体运移中的作用，CT 断层扫描原状土柱的大孔隙研究结果精准且具有普适性[27]，为污染物在大孔土壤中迁移的研究提供理论基础[28]。有研究发现，大孔隙弯曲度对水分运移有较大影响，但多为模拟弯曲度测算值，针对弯曲度的实际计算方法仍需研究[29]。本研究旨在采用非破坏性CT 断层扫描技术及结构方程模型（SEM），量化不同土壤的孔隙结构特性，定量分析不同土壤类型、植被类型条件下土壤大孔隙特征，以期为大孔隙数学建模[30]、土壤孔隙水分运移模型[10]、河岸包气带及边坡稳定性[31]等研究提供科学支撑，为不同土壤类型条件下人工林的恢复、合理配置、抚育管理等问题提供参考。

1 研究区概况

本研究的11 个原状土柱分别采集于典型褐土区（北京鹫峰试验场）、红壤区（重庆四面山试验场）、喀斯特地貌红壤区（云南建水试验场）。研究点基本情况见表1。

表1 试验区概况Table 1 Basic situation of the experimental sites

1）北京鹫峰试验场（116°08′～116°28′E，39°54′～40°10′N）位于华北平原西部，属于太行山系最北段与燕山山脉相会处的北京西山。海拔高度60～1153 m，受华北暖温带半湿润半干旱大陆性季风气候影响，年平均气温11.8 ℃，年平均日照时数2000～2800 h，年平均降水量550.3 mm。其中乔木群落为栾树（Koelreuteria paniculata）、侧柏（Platycladus orientalis）混交林，灌木为紫薇（Lagerstroemia indica），草本群落为紫花苜蓿（Medicago sativa）。

2）重庆四面山试验场（106°17′～106°30′E，28°31′～28°43′N）位于我国地势二级阶梯东缘三峡库区内。海拔高度500～1780 m，位于中纬度亚热带，受亚热带湿润季风气候影响，雨水充沛且雾气多。全年雨量较大，但时空分布差异较大，年平均气温13.6～18.4 ℃，年平均日照时数888.5～1539.6 h，年平均降水量1221.5～1530.7 mm[17]。其中灌木群落为细枝柃（Eurga loquaiana）。

3）云南建水试验场（102°55′～102°60′E，23°37′～23°44′N）位于云南省建水镇西南岩溶区羊街集水区，为典型的喀斯特盆地地貌。海拔高度1350～1700 m，受西南印度洋季风影响，可形成明显的干湿季。该地区季节性干旱严重，说明雨季（5 月至10 月）潮湿，旱季（11 月至次年4月）降水稀少。年平均气温19.8 ℃，地表年平均气温20.8 ℃，年平均日照时数2322 h，年平均降水量805 mm。其中乔木群落为桉树（Eucalyptus maideni），灌木群落为车桑子（Dodonaea viscosa）。

2 研究方法

2.1 样地选取

为实现数据量的充足性，获得了另外两位研究者的许可和原状土柱数据，进行大孔隙特性分析。试验选取了3 个样地，包括北京褐土混交林（B−F−1，B−F−2）、灌木（B−S−1，B−S−2）、草地（B−G−1），重庆红壤灌木（C−S−1，C−S−2），云南喀斯特桉树人工林（Y−EF−1，Y−EF−2）、灌木（Y−S−1，Y−S−2），所有样地均选自坡脚处，消除山体隆起高度对土壤特性的影响[32]。为减少误差，首先将土壤表层枯落物剥离，再采用按压法取原状土柱，同时根据对角线布点法采集环刀样本。

2.2 土壤性质

采用环刀浸透法获得土壤并根据公式（1）～（3）计算总孔隙度（P）、最大持水量（P1）、自然含水量（P2），在每个样地10、30、50 cm深处使用环刀钻将环刀压入土壤5 cm 深，所得环刀土样分别代表0～20、20～40、40～60 cm 土层的土壤，并分别取3 个重复样品。环刀体积（v）100 cm3把每个土样的环刀下盖（带网孔端）打开，放入滤纸，盖好后立即称质量（m1），取下环刀上盖，将环刀带网孔端放入水盆中，盆中水层高度至环刀上沿（不淹没），吸水 12 h 后取出环刀，盖上上盖，擦净水分，立即称质量（m2），放入 105 ℃烘箱内烘烤 24 h，盖好上盖，立即称质量（m3）。取出土样，称取带滤纸环刀质量（m4），计算公式如下[18-19]：

式中：P为总孔隙度，P1为最大持水量，P2为自然含水量，v为环刀体积，m1为原质量（含环刀、滤纸质量），m2为饱和土壤环刀质量（含环刀、滤纸质量），m3为干质量（含环刀、滤纸质量），m4为环刀、滤纸质量。

2.3 CT 扫描及图像分析

本研究的原状土柱（管内直径10 cm，高30 cm）中有5 个取自北京鹫峰试验场，2 个取自重庆四面山试验场，4 个取自云南建水试验场，各重复2 个样品，均采用管壁5 mm 的PVC 管采集，因金属管壁会对扫描结果产生影响。11 个原状土柱均采用IPT4106D 扫描设备进行断层扫描，4 MeV 加速器面阵CT 扫描系统，将原状土柱以0.2 mm/pixel 扫描为2048×2048 的CT 图像矩阵，进而使CT 图像堆栈通过VG Studio MAX 2.2软件对孔隙进行缺陷自动检测计算，利用高斯滤波，在不破坏土壤整体结构的情况下获得其内部孔隙三维结构，以及大孔隙的直径、体积、表面积、分布。

需要计算的大孔隙特征参数还包括大孔隙角度、扭曲度、大孔隙实际长度（Lti）和弯曲度（δ）。为计算这些大孔隙特征，假设某段孔隙如图1 所示，若将其无限拆分为n层，则该孔隙可近似看做一段线段，将每一维段的数值进行积分计算即可获得特征指标。

图1 大孔隙计算过程示意图Fig.1 Diagrammatic sketch calculation process of macropores

大孔隙角度（θ）为大孔分支的倾角[33]，公式如下：

式中：θ为大孔隙角度，x、y和z分别是x、y和z方向上的投影大孔长度。

扭曲度（ξ）可描述土壤内大孔隙的扭曲程度，为实际长度（L0）的土壤内共n层大孔隙角度（θi）的累积值，计算方法见公式（5）：

式中：ξ为扭曲度，n为土壤总层数，θi为第i层大孔隙角度的累积值。

第i层大孔隙实际长度（Lti）为假设大孔隙被拆分后每一微段孔隙为圆柱体，则实际大孔隙长度（Lt）为大孔隙体积除以横截面积，计算方法见公式（6）：

式中：Lti为第i层内实际大孔隙长度（mm），V为土壤孔隙体积（mm3），S为土壤孔隙表面积（mm2），r为土壤孔隙横截面半径（mm），Vi为第i层大孔隙的体积（mm3），d i为第i层大孔隙的直径（mm）。

弯曲度（δ）为大孔隙路径的卷积程度，可反应大孔隙通道内流动物体的运移路径长度，计算方法见公式（7）：

式中：δ为弯曲度，Lti为第i层内实际大孔隙长度（mm），Lli为第i层内垂直距离长度（mm）。

2.4 结构方程模型统计分析

为定量描述土壤大孔隙特征参数与土壤自然含水率、土壤最大持水量间的关系，采用结构方程模型（SEM），在考虑变量间的相关性的同时也会考虑误差，因此可以用来分析优先流的试验数据[34]。与传统的Pearson 相关系数统计方法相比，结构方程模型SEM 在分析相关性过程中将可测变量及潜在变量分开考虑，允许参数估计中的测量误差，减少由测量误差引起的潜在变量之间相关性估计的干扰[35]。同时，结构方程模型不仅可以得到因子间的相关系数，还可显示相关关系中的方向性，因此本研究选择结构方程模型分析各因子间的相关关系。SEM 由结构方程和测量模型组成，通常由3 个矩阵方程表示，见式（8）～（10）[36]。公式（8）反映了实际中无法测量的潜在变量之间的结构关系，公式（9）、（10）为测量模型，可以反映观测变量和潜在变量之间的关系。

式中：η为内在潜变量，B为内在潜变量之间的相关系数矩阵，Г为内在潜变量对外在潜变量的结构系数矩阵，ξ为外在潜变量，ζ为测量误差，Y为内在观测变量矢量，X为外在观测变量矢量，Λx为外在观测变量矢量的因子载荷，Λy为内在观测变量矢量的因子载荷观测变量向量，分别是外在观测变量和内在观测变量的度量。潜在变量可以通过测量模型由可测量变量反映出来。通过求解方程（8）～（10），可以得到每个潜在变量和可测变量之间的关系。在建立初始模型之后，路径系数表示变量之间关系的程度，通常用极大似然法计算。

3 结果与分析

3.1 不同土壤类型条件的大孔隙特征对比

3 个研究区植被群落土壤的三维可视化图像见图2。由图2 可知，这些大孔隙既连续又独立，而大孔隙的连续性在土壤的渗透过程中极为重要[37]。

图2 11 个样本的三维可视化图像Fig.2 3D visualization images of 11 samples

北京褐土、重庆红壤、云南喀斯特红壤3 种土壤条件下的灌木土壤大孔隙形态特征清晰。北京褐土和重庆红壤的灌木群落土壤大孔隙主要为纵向发育，通道形态清晰且分布较均匀，而云南的灌木群落土壤大孔隙主要为横向发育且多为小体积裂隙（孔隙体积<1 mm3），分布较均匀，这与Kung[38]的研究结果相似。

北京褐土、云南喀斯特红壤的乔木土壤大孔隙连通度高。北京褐土区的乔木群落土壤大孔隙明显下沉且纵向连续性强，而云南喀斯特地貌的人工林群落土壤大孔隙集中区域明显上移且破碎化程度极高，横向连续性强。

北京褐土的草本群落土壤大孔隙分层明显。草本土壤的大孔隙明显分为上、中、下3 层，上层近地表孔隙可能与草本植物浅根相关，大孔隙通道直径较小且不聚集成团，中层存在连接上、下层大孔隙的通道，但分布数量较少，而下层土壤中接近地下30 cm 深处的大孔隙均为垂直形态且分布较均匀，该孔隙通道的形成可能与地下裂隙或水分入渗侵蚀相关。

选取11 个样品中较具代表性的横截面见图3。北京褐土地区的大孔隙呈网状结构且密集，重庆红壤地区灌木的大孔隙呈通道状，较分散且孔隙直径明显大于北京褐土地区的灌木样地，云南喀斯特地区的人工恢复林大孔隙呈聚集状且覆盖面积较大，该地区土壤破碎程度较大。

图3 11 个样本的代表性横截面Fig.3 Representative cross section of 11 samples

不同土壤类型条件下灌木土壤大孔隙特征值对比见表2。灌木群落土壤大孔隙的平均直径为Y−S>C−S>B−S，标准差较小；孔隙总体积为C−S>B−S>Y−S，标准差较大；孔隙总表面积与累积扭曲度趋势相似，为B−S>C−S>Y−S，标准差较大；累积弯曲度与平均角度的趋势相反，标准差较大。

表2 灌木大孔隙特征值Table 2 Characteristic value of macropore of shrub

不同土壤类型条件下乔木土壤大孔隙特征值对比见表3。乔木群落土壤大孔隙的平均直径、累积扭曲度与平均角度的趋势相似，为B−F>Y−EF，标准差为Y−EF>B−F；孔隙总体积、孔隙总表面积与累积弯曲度趋势相似为Y−EF>B−F，标准差为Y−EF>B−F。

表3 林地大孔隙特征值Table 3 Characteristic value of macropore of forest

3.2 不同植被类型条件的大孔隙特征对比

不同植被类型条件下北京褐土区的大孔隙特征值对比见表4。北京褐土区的林地、灌木、草地土壤大孔隙直径平均值及平均角度均为乔木最大，灌木和草本相似且较小；土壤孔隙总体积为B−S>B−G>B−F，灌木土壤的标准差较大；孔隙总表面积、累积弯曲度、累积扭曲度均为B−G>B−S>B−F。

表4 北京褐土不同植被类型条件下大孔隙特征值Table 4 Macropore characteristic value of cinnamon soil in Beijing

不同植被类型条件下云南喀斯特红壤区的大孔隙特征值对比见表5。云南喀斯特红壤区的林地、灌木土壤大孔隙直径平均值、土壤孔隙总体积、孔隙总表面积及平均角度均为Y−EF>Y−S；累积弯曲度、累积扭曲度均为Y−S>Y−EF。

表5 云南喀斯特不同植被类型条件下大孔隙特征值Table 5 Macropore characteristic value of karst soil in Yunnan

3.3 基于结构方程的土壤孔隙特征参数相关性计算

结合表1～5 数据可知，不同土壤类型、不同植被类型条件下的土壤大孔隙平均直径、孔隙总体积、孔隙总表面积、累积弯曲度、累积扭曲度、平均角度等特征参数与土壤孔隙度、土壤自然含水量、土壤最大持水量之间存在相同或相反的变化趋势，为定量研究大孔隙特征参数与土壤含水量间的关系，本研究利用AMOS 7.0 软件中的结构方程模型（SEM）建立各指标间的关系，分析土壤总孔隙度、最大持水量、自然含水量、平均直径、孔隙总体积、孔隙总表面积等可测变量与累积弯曲度、累积扭曲度等潜变量间的关系[34]，模型结构路径经过多次调整后得到图4。各特征参数之间的路径系数、标准误差、临界比和模型影响因素的显著性见表6～7[39]。本研究采用绝对拟合指数和相对拟合指数对模型进行检验，表7 中SEM 各指标均符合标准，且临界比绝对值大于2，说明可以利用该模型研究各特征参数的相关性。在本研究中，采用绝对拟合指数和相对拟合指数对模型进行检验（表7）。每一个都符合标准，说明图4 的结构方程模型可用。经SEM 去误差化分析土壤总孔隙度、最大持水量、自然含水量、平均直径、孔隙总体、孔隙总表面积、累积弯曲度、累积扭曲度共8 个指标之间的关系。由表6 可知，孔隙总体积与平均直径呈显著负相关（标准化估计为−0.56），说明在相对固定的土体内大孔隙总体积的增大可能抑制大孔隙平均直径的大小。累积弯曲度与孔隙总表面积呈显著正相关（标准化估计为0.69），说明弯曲度高的土体内，大孔隙的累积表面积更大，水分流经该部分土体时可以与土壤有更充分的接触。

表6 SEM 协方差结果Table 6 Regression results of SEM variables

表7 SEM 整体适配度评价及拟合结果Table 7 Evaluation and fitting results of SEM overall fitness

图4 土壤孔隙参数典型特征值的结构方程Fig.4 SEM of preferential flow characteristic parameters

4 结论与讨论

同一立地植被条件下不同土壤类型的土壤孔隙度及弯曲度均与土壤持水量成正比，即乔木和灌木土壤在不同土壤类型条件下，土壤总孔隙体积与累积弯曲度越大，土壤最大持水量和自然含水量越小。云南喀斯特红壤区乔木群落土壤大孔隙的总体积、总表面积、累积弯曲度、累积扭曲度均大于北京褐土区乔木群落土壤，但两地区的灌木群落土壤大孔隙特征值相反，云南喀斯特红壤区的乔木土壤大孔隙总体积和表面积分别是北京褐土区乔木土壤的2.25 倍和5.64 倍，说明云南喀斯特红壤区人工恢复地区乔木群落土壤的大孔隙分布较北京褐土地区更大，喀斯特红壤区人工林土壤破碎程度高且水平连续性强，同时其内部孔隙表面积更大，因此可有效增加土壤的蓄水能力，减少水分滞留在地表的时间[40]，进而减少水分挥发，提高地表水分的利用。云南喀斯特红壤区乔木土壤的累积弯曲度和累积扭曲度分别是北京褐土区乔木土壤的2.78 倍和0.17 倍，说明北京褐土区和云南喀斯特红壤区人工恢复乔木群落土壤的大孔隙扭曲度近似，但水分在云南喀斯特红壤区中下渗的路径更长，因此在土壤中的滞留时间可能更长。

北京褐土和重庆红壤的灌木群落大孔隙分布均匀且纵向发育显著，大孔隙网络结构密集，云南喀斯特红壤的灌木群落土壤大孔隙横向发育显著。重庆红壤区灌木群落大孔隙总体积最大，分别是北京褐土区灌木和云南喀斯特红壤区的1.22倍和6.28 倍，说明重庆红壤区灌木群落的土壤大孔隙空间最大，存在大量蓄水的可能性。但重庆红壤区灌木的大孔隙表面积不是最大的，北京褐土区灌木大孔隙总表面积分别是重庆红壤区灌木和云南喀斯特红壤区灌木的1.44 倍和1.46 倍，说明重庆红壤区灌木的大孔隙表面积与云南喀斯特红壤区灌木接近，水分在北京褐土地区灌木群落的土壤中运移过程中可能会接触到更多的土壤，若浅层土壤中存在可溶解的有毒物质，则其溶解在水中并渗入土壤的可能性较大。Saravanathiiban 等[41]和徐宗恒等[42]的研究表明，水分在土壤中的入渗速度会在孔隙通道弯曲度高处降低，因此弯曲度高的立地类型更易存储水分并供植物根系生长使用。北京褐土区灌木的累积弯曲度是重庆红壤区灌木和云南喀斯特红壤区灌木的5.84 倍和0.24 倍，说明北京褐土地区和云南喀斯特人工恢复地区灌木群落土壤的大孔隙弯曲度相近，但重庆红壤地区灌木群落土壤的大孔隙弯曲度很小，而北京褐土区灌木的累积扭曲度是重庆红壤区灌木和云南喀斯特红壤区灌木的1.62 倍和4.52 倍，说明北京褐土区灌木的累积扭曲度和累积弯曲度最高，大孔隙通道的形状可近似于螺旋体，重庆红壤区灌木的累积弯曲度低且累积扭曲度高，大孔隙通道的形状可近似于螺丝钉体，云南喀斯特红壤区灌木的累积弯曲度高但累积扭曲度低，大孔隙通道的形状可近似于滑梯状。

同一土壤类型条件下不同立地植被的土壤孔隙特征差异显著。乔木土壤大孔隙垂直连续性强，乔木群落大孔隙平均直径约为灌木群落大孔隙平均直径的1.33 倍，接近草本群落大孔隙平均直径的1.87 倍，说明草本群落的土壤孔隙细小，地表水分的入渗可能会受其土壤孔隙细小的影响。乔木群落的土壤大孔隙分布最广泛且垂直连通性极强，存在明显的聚集现象，水分从地表可以优先沿大孔隙进入地下[43]，并快速运移至深层土壤，因此乔木群落相对于草地更不易形成地表径流，且蓄水性更好；灌木群落土壤大孔隙特征值与乔木相反，灌木群落的土壤大孔隙分布较均匀，且孔隙直径较小，通道形态清晰；草本群落土壤大孔隙总体积小，但累积扭曲度最高，且草地的土壤大孔隙分布极不均匀且出现分层隔离现象，连通性更差，因此草地无法将地表水分优先运移至土壤内部，因此极易形成地表积水和地表径流[44]。

根据不同土壤类型、植被群落条件下土壤大孔隙的各特征参数，褐土区的累积弯曲度和累积扭曲度均较高，故大孔隙通道的形状可能近似于螺旋体；红壤区累积弯曲度低但累积扭曲度高，故大孔隙通道的形状可能近似于螺丝钉体；喀斯特人工恢复地区土壤累积弯曲度和扭曲度都低，故大孔隙通道的形状可能近似于滑梯状。草本群落的土壤大孔隙旋转程度远高于其他立地条件。

经结构方程模型验证，土壤大孔隙弯曲度的增加会促进土壤大孔隙总面积的增大，该结果与Saravanathiiban 等[41]的结论一致。大孔隙的体积与平均直径呈负相关，与Zhang 等[43]的结论一致，表明大孔隙只占土壤总孔隙度的一部分。弯曲度会抑制土壤自然含水量和最大持水量，该结果与Meng 等[33]的结论相似，表明土壤大孔隙结构的复杂性可能增加水分流失情况发生。

土壤大孔隙结构的三维研究对于土壤水分运移具有决定性影响，但目前对于土壤大孔隙结构的研究仍较少，主要受到CT 扫描技术以及样本运移的限制，但本研究主要对比3 种土壤类型条件下乔灌草的土壤大孔隙形态特征对比，不同条件下的土壤大孔隙结构具有较高辨识度，有较明显的区别和特征，但仍需进一步更广泛的实验验证和理论完善。