基于CT扫描技术预测露天煤矿区不同压实土壤水分特征曲线

王金满,张丽娜,冯 宇,荆肇睿

(1.中国地质大学(北京) 土地科学技术学院，北京 100083；2.自然资源部 土地整治重点实验室，北京 100035)

土壤压实是一种严重的环境地质现象，尤其是大型机械碾压在工程和农业中的应用。压实作用被定义为土壤致密化和土壤孔隙度的减小，会改变土壤结构，增加土壤强度，并降低土壤的水力传导性，这又反过来对粮食生产、土壤生物栖息和生态功能产生负面影响。一般来说，当土壤密度在黏土上超过1.47 g/cm，在淤泥上超过1.75 g/cm时，植物根系部分通气不良，土壤渗水吸水能力不足，植物生长可能会受到限制。对于土壤破坏严重的矿区，大型机械造成的土壤压实极为突出，成为限制当地农业发展的重要因素。

露天开采会对当地生态环境造成严重而迅速的影响，挖掘和倾倒导致陆地生态系统的快速变化，土壤剖面层紊乱并加剧土壤水土流失和退化。其中严重的压实效应导致土壤孔隙结构变化，限制了土壤的养分循环、储水、气体交换和能量运输。研究土壤水力特性对恢复矿区土壤质量有很大帮助。水分特征曲线是模拟土壤水和溶质传递的主要工具之一，大量实验结果表明，土壤变形和压实的程度对水分特征曲线有决定性作用，通过降低土壤大孔隙度和孔径分布限制土壤含水量。此外，土壤温度、质地也会对其产生影响。不同压实效应的孔隙参数和水力特性的差异需要量化，以诊断矿区因机械负荷导致的土壤退化，做好土地复垦工作。传统的实验室方法，对于孔隙的描述缺乏足够的细节，且测定水分特征曲线耗时并存在许多不确定性。因此，表征和量化土壤孔隙结构是确定压实作用对土壤结构和水力性质影响的关键。亟需新高精度技术的支持。

X射线CT扫描技术已被证明是一种非常强大且非破坏性的方法，以高分辨率来量化土壤结构。近年来，CT扫描技术已被用于确定溶质运动、孔隙度和植物根系发育等方面。该技术在不破坏样品的情况下分析研究对象的内部特征，具有快速成像的优点。CT方法在分析土壤科学参数方面具有广阔的应用前景，如土壤孔隙特征、孔隙形状、大小、分布状况；土壤水力特性、饱和导水率、水分特征曲线。从微观尺度更好地解释人为过程(如压实)引起土壤性质的变化。

由于土壤大孔隙结构与多种土壤性质和过程密切相关，例如水力性质和渗透过程，因此CT扫描技术的应用有助于这些过程的可靠预测。以往土壤大孔隙的研究多集中在农田、牧场和森林的土壤压实中，对矿区重型机械造成的严重压实研究有限。因此笔者通过CT扫描技术重建矿区不同压实土壤的三维大孔隙结构，量化土地大孔隙结构特征，并结合物理经验模型预测土壤水分特征曲线，探讨压实对矿区土壤水力特性的影响机理。

1 研究区概况

研究区位于山西省平朔露天煤矿，地理坐标为112.17°～113.50°E，39.38°～39.62°N。土壤样本采集点位于平朔矿区安太堡露天煤矿的排土场，经历了堆积和大型机械反复碾压的过程。该区域属于典型的温带干旱至半干旱大陆性季风气候和脆弱的生态系统。土壤类型主要是栗钙土，土质偏砂土质，土壤结构差，土壤密度介于1.27～1.74 g/cm，有机质质量分数低，土壤有机质质量分数一般在5.0～9.0 g/kg，对水蚀和风蚀的抵抗力弱，由于地处黄土高原干旱区，缺少良好的植被覆盖。此外，排土场内高度压实使土壤取样过程困难，且提取土壤岩心通常会对样品造成巨大的干扰。为了避免误差，在实验室进行模拟实验，根据排土场土壤的压实程度，实验室采取手工压实制备土壤模拟样品，将野外排土场表层土壤收集，并作为压实实验的材料送回实验室。

2 材料与方法

2.1 样品采集

模拟实验的土样采集自安太堡露天煤矿排土场0～30 cm土壤(图1)。取样后，在室内风干、压碎，混合并过2 mm孔筛。通过手工压实将土壤装进直径60 mm、高70 mm的玻璃管中，每组密度设置3组重复。为得到均匀压缩的土芯，每个岩心的压实分为4层进行，从300 mm的固定高度锤压。样品分为2组，每组18根土柱，密度为1.3～1.8 g/cm，密度间隔为0.1 g/cm，其中G1.4略微压实，其压实水平与现场未受干扰的原状土壤相似，G1.3和G1.4被认为是未压实的土壤，G1.5，G1.6，G1.7和G1.8被认为是压实土壤。2组实验样品中一组用于实验室离心机实测土壤水分特征曲线，另一组送至中国航天特征材料与加工技术研究院实验室进行CT扫描。

2.2 水分特征曲线的测定

将按照设定密度准备的土样在水中浸泡饱和24 h，称重饱和后土样，置入离心机中，设定各基质势对应的转速(表1)，离心机设置旋转时间为60 min，取出离心桶称重，依此进行下一转速，离心旋转结束后，将所有样品在120 ℃烘箱内干燥至恒重，计算个基质势下土壤体积含水量。最后得到土壤水分特征曲线(SWCC)。

图1 土壤样品Fig.1 Soil sample

表1 离心机转速与基质势

2.3 CT扫描

在螺旋扫描模式下扫描岩心，为提高图像精度，扫描土柱中部和5～55 mm区域，扫描图像的分辨率为30 μm ×30 μm×30 μm，扫描间隔为连续0.05 mm。CT扫描配置参数值：管电压100 kV, 扫描分辨率25 μm,像素数2 284×2 284。CT扫描后，有1 200张灰色图像。然后，通过计算机处理，获得土壤切片图形，基于Avizo对图像切片进行三维重建。

在土壤孔隙研究中，一般将土壤孔径大于30 μm的孔隙定义为土壤大孔隙，相对于土壤中的微孔隙，土壤大孔隙更容易受到土壤压实的影响。

在ImageJ软件中处理切割采集的图像堆，尺寸为10 mm×10 mm。首先排除岩壁的空隙，以尽量减少射束硬化带来的影响。采用最大熵阈值算法分割图像，通过选择最大类间熵来确定合理的阈值，结合目视检查确保正确性。分割后，使用Avizo软件的孔隙网络模型重建并可视化大孔隙结构，获取土壤大孔隙结构，不同的颜色代表不同的孔径。图像处理流程如图2所示。

图2 图像处理技术流程Fig.2 Workflow of image processing

模型中土壤孔隙的当量直径()的计算公式为

(1)

其中，为土壤中大孔隙体积，μm。由于图像分辨率的限制，本研究中当量直径()为大于30 μm的大孔隙。在Avizo软件中，采用Volume Fraction模块提取土壤大孔隙度结构信息。

土壤大孔隙度()计算公式为

=

(2)

式中，为模型中统计的大孔隙体积；为土样样品的总体积。

2.4 水分特征曲线预测模型

在预测土壤水力特性模型中常用的有物理经验模型。物理经验Arya-Paris模型最早于1981年提出，该模型根据土壤密度和土壤大孔隙的分布特征来模拟土壤水分特征曲线。其基本假设是土壤粒径组成与大孔隙半径相关。基本原理是计算每个孔隙级别范围内大孔隙的含水量及其对应的基质势，因此，获得土壤含水量和基质势，预测整个范围内的水分特征曲线。

该研究中，土壤大孔隙半径通过CT扫描土壤重建土壤三维结构后获得。土壤孔隙分布的直方图可以在Avizo 软件中直接获得。据此，可以估算出土壤大孔径的分布。具体公式为

=4cos(1-)

(3)

式中，为每组直方图内土壤大孔隙平均直径；为水的表面张力，N/m；为水分接触角；为水密度，kg/m；为重力加速度，m/s；,为表征土壤孔隙形状的经验参数。

土壤大孔隙分布曲线将土壤孔径分为份：

=(=1，2，3，…，)

(4)

式中，为级土壤中的大孔隙体积总和；为级土壤的质量；为土壤颗粒密度。

土壤大孔隙度()可以通过以下公式计算得到：

=(-)

(5)

其中，为土壤密度；将土壤的累积粒度分布曲线分为份，相应的土壤大孔隙半径也分为份；假设在水体流动过程中，土壤水先填充满土壤小孔隙，然后填充较大孔隙，因此，土壤含水量小于第个土壤孔隙的积累，其计算公式为

(6)

其中，为第级土壤孔隙充满水时的体积含水量；为第级土壤孔隙体积；为单位土壤的体积，可通过以下公式计算得到：

(7)

其中，为级土壤颗粒的质量，因此式(6)可改为

(8)

取2个相邻等级的土壤体积含水量的平均数作为对应等级的较大孔隙体积水分含量：

(9)

根据毛细管理论，可以计算相应的基质势：

=2cos()

(10)

式中，为第级土壤孔隙对应的基质势；为第级土壤孔隙半径。

本研究中CT扫描所能识别的最小孔径是30 μm，因此模型所能预测的最大基质势是102 cm。由于实测值和预测值的基质势范围不同，无法比较整个范围内CT扫描技术结合物理经验模型所提出预测方法的预测能力。因此运用Matlab软件基于Van Genuchten模型对预测值进行整个范围内的拟合，将得到的水分特征预测曲线与实测值的拟合曲线进行误差分析。

用实验结果的实测值和基于CT扫描技术得到预测值的相关系数和均方根误差 (RMSE)来描述土壤水分特征曲线预测性能。

(11)

研究中的统计分析通过Spass软件进行。

3 结果与分析

3.1 实测土壤的水分特征曲线

图3显示了实验室离心法测定的土壤水分特征曲线对压实作用的响应程度。观察发现，在同一压实状态下，随着土壤基质势的增加，土壤含水量先迅速降低而后缓慢降低。这种在低基质势显著降低，高基质势降低趋势缓慢的现象，可能与在土壤压实过程中，大孔隙对压实作用的敏感性有关，随着压实程度的增加，土壤大孔隙中的水分优先排出。土壤水分特征曲线在低基质势时不同密度土壤的含水量差距较大，在中基质势时，除密度为1.8 g/cm时，其他压实状态下，土壤含水量的差距较小，在高基质势下，不同密度土壤含水量又显出差异。比较同一基质势下，随着密度的增加，土壤含水量减少，且含水量的变化范围减小。压实程度增加，土壤结构发生重组，进一步统计土壤大孔隙发现，当密度较小时，大孔隙数量占比要比密度较高时提高30%～50%。这表明，土壤大孔隙特征与土壤含水量有直接关系，可进一步基于土壤孔隙结构参数预测土壤水分特征曲线。

图3 水分特征曲线实验室测定Fig.3 Soil water characteristic curves by experiment

3.2 不同密度土壤的大孔隙特征

..土壤大孔隙分布特征

矿区土壤不同压实程度对土壤大孔隙的分布存在显著差异，从图4可以看出，不同密度下土壤大孔径分布及其孔隙总体积占比特征。在0～500 μm孔径范围内，各密度下大孔隙体积占比分别达到了11.66%，24.86%，49.11%，52.43%，95.84%，100%。

尽管不同密度时土壤大孔隙数量存在差异，但孔径在50～150 μm内的大孔隙数量最多，且随着压实程度增加，土壤孔径范围降低，且同一孔径范围内的土壤大孔隙数量减少。当密度增加到1.8 g/cm时，土壤中大孔隙的孔径均小于500 μm；进一步比较各孔径范围对总体积的贡献率，在各压实状态下，土壤大孔对体积贡献率均处于较高水平，当密度在1.3～1.4 g/cm内，随着孔径增加，对孔隙总体积贡献率增加，1.5～1.6 g/cm内，随孔径增加，体积贡献率波动增加，并当孔径较大时，贡献率达到峰值，1.7～1.8 g/cm内，随孔径增加，体积贡献率先增加后减少，且100～200 μm孔径范围的大孔隙对其贡献率最大。随着密度增加，大孔隙数量在可测量的范围内有效减少，这表明，压实作用优先破坏土壤中的大孔隙结构，在土地复垦管理工作中，可通过犁耕等疏松表土，减少压实对土壤带来的负面影响，改善土壤环境，提高土壤质量。

图4 压实土壤大孔隙数量与体积占比Fig.4 Number and volume proportion of macropores in compacted soils

..压实对土壤大孔隙结构的影响

土壤样品的孔隙度见表2。其中土样的平均孔隙度是各土样检测到的大孔隙度和微孔隙度之和，统计发现平均估计孔隙度随密度的增加而减小，当密度为1.8 g/cm。时，孔隙度降低了37.02%；对比土壤大孔隙度与微孔隙度对压实作用的响应，大孔隙度随密度增加迅速下降(图5)，微孔隙度随土壤密度增加下降不明显。土壤中的连接孔隙只有在密度为1.3 g/cm才存在，进一步说明压实作用对土壤大孔隙有显著的影响。

表2 估算和检测土壤样品的孔隙度

图5 不同密度土壤大孔隙度Fig.5 Macroporosity of soils with different bulk densities

图6显示了土壤大孔隙数量和密度间的关系，2者间的相关性高达0.945，随着密度的增加，大孔隙数量减少，当密度达到1.8 g/cm时，土壤大孔隙数量降低至1 254，降低了78.72%。

土壤中的欧拉数是描述压实土壤连通性的有效指标，其值越大，表示孔隙间的连通性越低。由图7可知，随着密度的增加，土壤孔隙间的连通性逐渐减小，且当密度高于1.5 g/cm，时，平均欧拉数缓慢增加，这表明当压实程度的增加，土壤中的连接孔隙迅速减少，当压实达到一定程度时，连接孔隙极少，压实对土壤孔隙结构参数影响极为敏感。研究发现。土壤大孔隙结构特征与土壤压实存在着强相关，CT扫描技术在微观尺度更好解释了压实对土壤的负面影响，矿区土壤面临的土壤压实威胁更多，机械的多次碾压，矿山废弃物的堆积等，为避免压实带来的各种威胁，有必要根据实际情况采取恢复措施。

图6 不同密度土壤大孔隙数量Fig.6 Macropores number of soils with different bulk densities

图7 不同密度土壤平均欧拉数Fig.7 Average euler number of soils with different bulk densities

3.3 CT扫描技术预测土壤水分特征曲线

基于CT扫描技术结合物理经验模型对土壤水分特征曲线进行了预测，图8显示了不同压实程度下土壤水分特征曲线的实测值和预测值的拟合曲线。结果显示，不同压实状态下土壤水分特征曲线预测值和实测值的形状变化一致。在可预测的基质势范围内，土壤体积含水量随着基质势的增加先缓慢减少后迅速减少，当基质势达到一定高度时，土壤体积含水量减少趋势又趋于平缓。对比同一基质势下，土壤保水能力随压实增加而下降。通过对于预测值和实测值可验证本研究所提出预测方法的可行性，在模拟结果中，预测值和实测值接近，2者的变化趋势一致。密度越小，拟合效果越好，说明预测精度越高。当压实处于1.3～1.6 g/cm时，在整个基质势范围内，预测值先高于实测值，后低于实测值，整体上预测能力较高，当压实达到1.7～1.8 g/cm时，在低基质势时，土壤含水量预测值与实测值的偏差较大，随着基质势的增加，2者间的差值逐渐减小。预测的土壤水分特征曲线的相关性均大于0.8，见表3，统计了实验测定值的拟合曲线与模型预测值拟合曲线的均方根误差(RMSE)，其值越小，表明CT扫描技术的预测能力越高。本研究结果显示，基于CT扫描技术获取土壤大孔隙特征参数，结合物理经验模型预测土壤水分特征曲线达到较高的精度，可应用于实际问题的研究。

图8 土壤水分特性曲线预测值与实测值的拟合Fig.8 Fitting between predicted and measured values of soil water characteristic curve

表3 不同密度土壤水分特征曲线的预测精度

4 讨 论

4.1 压实对土壤水力特性的影响

对于传统实验室测定的土壤水分特征曲线，压实程度的增加有效减少了土壤中的含水量，且随着离心力的增加，土壤体积含水量先快速减少，当达到较高基质势时，含水量变化不再明显。通过密度与大孔隙特征相关性的分析，有效说明土壤中的大孔隙分布显著影响土壤水分特征曲线变化。压实会导致土壤紧实度增加，一旦压实作用展开，其内部颗粒会不断挤压，重新排列组合，首先大孔隙中的水分优先排出，这与土壤体积含水量迅速减少一致。土壤孔隙作为微观尺度物质交换的场所，可通过翻耕、添加生物碳等维持土壤大孔隙结构与孔隙间的连通性，这对于压实作用对孔隙结构、植被生长发育、农业生产带来的负面影响有缓解作用。其次，土壤体积含水量还会受到土壤质地、温度等因素的影响，这可作为后续的研究方向。

然而也有研究指出，适当的压实会增加土壤含水量，如在沙地中，压实对于其有一定的保水作用，可以防止水分蒸发，以及土壤膨胀引起的土壤孔隙结构和孔径分布潜在的可逆变化。黏土在压实过程中，土壤含水量减少的程度要小于其他土质，相关文献记载，当密度从1.2 g/cm增加到1.6 g/cm时，黏土质量分数为17%的土壤体积含水量会下降50%～60%，而黏土质量分数为83%的土壤含水量仅下降20%～29%。

矿区重型机械压实不仅对矿区水分平衡和流域生态水文安全造成威胁，也影响复垦区植被恢复效果。土体表面压实，导致水分入渗率降低，地表径流增加，加剧水土流失风险，同时，土体与大气物质和能量交换受阻，植被生长受限，降低作物产量，对玉米种植的压实研究表明，土壤压实显著降低玉米产量，经过中等和严重压实后的玉米产量分别减少25%和50%。一些植物在经过土壤压实之后能够迅速恢复，但也有植被则无法恢复到正常生长状态。平朔矿区植被生长状况较好时土壤密度均值维持在1.47 g/cm，但大多数区域压实程度高于该水平。有必要预防土壤过度压实，为植被生长发育提供良好的土壤环境：首先尽量避免机械作业，其次，对土层进行适当的翻耕，合理施肥并添加土壤改良剂，如生物碳的添加，有利于恢复土壤结构，改善土壤质量。

4.2 CT扫描技术预测土壤水分特征曲线的精度

本研究基于CT扫描技术结合物理经验Arya-Paris模型预测的土壤水分特征曲线结果具有较高的精度。相比传统的实验室方法获取土壤水分特征曲线，CT扫描技术结合土壤三维重构孔隙结构预测方法更加便捷准确，是可行的。同时发现在同一基质势下，预测精度随着土壤压实程度增加而有所减小，这主要与压实后土壤大孔隙减少、CT扫描仪精度有限有关，限制了模型的预测精度。后续考虑用高分辨率CT扫描仪进行土壤结构中微孔的测定，来研究各种土壤压实条件下土壤水分特征曲线的预测。同时整个土层的分布并不是均匀的，存在一些致密层影响水分运动，微观尺度的土壤水分运动是十分复杂的，在这种情形下，只考虑土壤大孔隙的分布特征无法完全诠释土壤水分的运动规律，土壤的质地、连通性、异质性均需全面考虑在内。

CT扫描仪获取的大孔隙结构参数可有效解释压实状态下土壤水分运动的规律，但其流动的变异性有待进一步探索，实现模型预测能力的优化。

5 结论与展望

(1)通过CT扫描技术三维重建土壤大孔隙结构，发现土壤孔隙中数量最多的是50～150 μm的大孔隙，压实作用降低了土壤大孔隙度，减少大孔隙数量，使土壤孔隙间连通性降低，同时大孔隙对压实作用更加敏感，对于矿区压实要以预防和避免大规模机械压实为主，同时可进行翻耕和添加改良剂等措施。

(2)实验室测定了土壤水分特征曲线，发现压实作用通过土壤孔隙结构影响土壤水分的运输及储存能力，相比高基质势时，在低基质势下，土壤含水量下降速率快，且随土壤密度增加而降低，表明压实会降低土壤的持水能力。

(3)研究提出的土壤水分特征曲线预测方法取得了较高的精度(>0.89)，相比传统的实验测定方法，省时省力。同时CT扫描仪精度的提升，可观测到更小孔径土壤，更能真实反映土壤水分状况，提高模型的预测能力。

当前，通用CT扫描仪的分辨率为30～100μm，在获取土壤孔隙结构时有一定的限制，图像分辨率在量化孔隙结构中起着决定性作用，有必要应用合适分辨率的CT扫描仪识别土壤中的微孔隙结构，探索更小尺度孔隙结构对土壤水分运动的影响。同时，现有的经验模型预测土壤水力特性仍存在一定的局限性，一方面各模型中的土壤参数各不相同，部分获取难度大，另一方面相关模型对土壤孔隙结构做了很大简化，相比真实的土壤水分运移，吸排水过程会有所差异，进一步完善合适的预测模型是今后的研究方向。