基于CT成像技术的盐渍土壤孔隙结构识别与分析

洪明海，黄介生，曾文治，伍靖伟，王雪扬(.武汉大学水资源与水电工程科学国家重点实验室，武汉 43007；.江苏省盐城市第三人民医院，江苏 盐城 400)

0 引 言

土壤是植物根系从中获取营养物质和水分的基础，是植物赖以生存的根基，而土壤孔隙结构和孔隙分布是影响土壤质地、肥力、含水量等土壤性质的重要因素。土壤孔隙是土体内部的空隙，是土壤结构的重要组成部分，土壤孔隙的数量、大小、分布及其连通性对土壤持水能力有着十分重要的作用[1]，影响着土壤有机污染物的吸附/解吸过程[2]。已有的研究多关注降雨[3-6]、土壤改良措施[7]、耕作方式、种植模式等对土壤孔隙的影响，很少考虑盐分对土壤结构的影响。然而，土壤盐渍化已经成为土地退化和影响作物生长的主要威胁。位于华北的内蒙古河套灌区，是我国最大的自流灌区，由于强烈的蒸发、较少的降水量和较浅的地下水位，约70%的耕地受到土壤盐渍化的影响。因此，有必要就盐分对土壤结构的影响机理展开研究，这对进一步认识土壤和植物的作用过程、改良和利用盐渍土有着重要的意义。

要定量研究盐分对土壤孔隙结构的影响机理，必须获取高精度的原状土壤孔隙结构并对其进行预处理和提取分析。获取土壤孔隙度最为经典的方法是根据土壤容重和密度进行计算[8]，但是这种方法的工作量大、耗时长。另一种较为常用的获取原状土壤孔隙结构的方法为序列数字图像技术[9]，但该方法操作复杂并且获取的切片图像精度较低，对原状土壤带有一定的破坏性。近年来，CT成像技术[10]以其快速、准确、便捷等优点广泛应用于医学、农业、生物等多个领域，由于CT成像技术可以在不破坏土壤样品的情况下对研究对象的内部结构进行扫描和测定，且成像及分析速度快，操作过程相对简便，该技术已经被逐渐应用于土壤孔隙结构的可视化和定量化研究[11]。本文通过CT扫描获得内蒙古河套灌区6个不同盐分水平下原状土壤的连续切面图像，利用数字图像处理技术并结合Arcgis 10.1软件以及python程序批量提取连续图像的土壤孔隙结构参数，并对获取的原状土样图像和孔隙参数进行对比分析，探求土壤盐分对土壤孔隙的影响机制，为盐渍土的改良和合理利用提供借鉴。

1 材料和方法

1.1 供试土壤及CT扫描

供试土样取至内蒙古自治区五原县义长试验站永联基地(108°00′16.14″ E，41°04′11.57″ N)，海拔约为1 030 m，地处黄河北岸，为河套灌区(处106°20′ E～109°19′ E， 40°19′ N～41°18′ N)盐渍土地带。

首先利用土壤三参数传感器(W.E.T Sensor, DELTA-T Co. Ltd)选择6块不同盐分水平的田块，并对每个田块取样带回实验室用雷磁DDSJ-318电导率仪准确测定土壤水溶液的EC1∶5值，最后通过公式[12,13]换算，以饱和土壤水浸提液的电导率ECe表示，各田块盐分值如下：田块1(S1:ECe=1.51 dS/m)；田块2(S2:ECe=4.33 dS/m)；田块3(S3:ECe=15.33 dS/m)；田块4(S4：ECe=21.31 dS/m)；田块5(S5：ECe=28.89 dS/m)；田块6(S6：ECe=54.92 dS/m)，其余田块各层ECe值如表 1所示。

然后利用Φ50 mm×50 mm的PVC管在上述6个田块分别采集原状土样，将取出土样密封后用冷藏箱运至江苏省盐城市第三人民医院进行CT扫描，所用CT机器型号为Discovery CT750 HD，CT扫描参数为：电压140 kV，电流220 mA，层厚0.63 mm，设置窗宽(Window width)、窗位(Window level)为1×580，对每个土样连续扫描80张CT图片，图像输出DICOM格式医学图片[图1(a)]。

从表1 可看出，除田块5(S5)的20～30 mm的17.61 dS/m和田块6(S6)的0～10 mm的68.97 dS/m两个值比较异常外，其余各田块各层ECe值都比较相近。而按照美国农业部[14]根据饱和土壤水浸提液电导率ECe值对土壤盐分水平的划分，所取土样盐分水平也列于表 1，2个非盐田块(S1、S2)，1个中盐田块(S3)，3个高盐田块(S4、S5、S6)，无低盐田块。

表1 各田块不同深度电导率值

1.2 图像处理

由于考虑扰动因素，去除土样表层的5张和土样底部的5张CT图片，选择余下的70张CT图像作为此次分析：①需要利用专业的医学图像处理软件Sante DICOM Viewer FREE软件打开CT扫描得到的DICOM医学图片，并进行模式选择和阈值调整处理，最终转化为tiff格式保存[如图1(b)所示]；②为了能使同一批图片都能用Arcgis10.1地理软件进行处理，先用ArcCatalog 10.1对tiff格式的图片构建金字塔显示(Buid Paramids)；③再利用Arcgis 10.1对图片定义投影(Define Projection)；④使用裁切工具将tiff格式的栅格图片剪切(Clip)为我们所要的半径为20 mm的研究区域(以图像中心为原点)；⑤用栅格转多边形(Raster to Polygon)工具将校正后的图片矢量化，得到如图1(c)所示的图片；⑥通过多边形的属性表提取出每一层孔隙的数量、单个孔隙的周长、单个孔隙的面积等土壤孔隙结构的基础计算参数作为后续的分析使用。

图1 土层切片中的孔隙分布情况

然而，在实际操作过程中，这一列的操作繁琐、易出错、处理周期长，且中间可能出现的问题较多，这对于大量的田间试验的图像处理是难以接受的；另外在矢量化后的图片中往往会产生大量的小多边形，这主要是因为CT得到的土样图片中会有很多微小的小孔隙还有非二值化的一些图像域，都会导致这种错误，比如，某个土样剖面的实际孔隙数只有52个，通过矢量化后得到的多达613个，这对于孔隙结构参数的提取有很大的影响，精度会显著下降。针对这些问题，在图像处理初期，用python编程统一将420张tiff图片编号排序(如：“S1_1.tif”表示的是盐分水平S1土样的第一张CT扫描图片的tiff格式)，并除去那些因为黑白阈值而导致的多边形使计算结果更加准确，由于Arcgis 10.1版本可外置python窗口处理，因此，前面所介绍的整个过程均通过python编程调用Arcgis 10.1中相应的模块实现批量处理，且后面的土壤孔隙结构的基础参数提取(见1.3)亦通过python编程实现并保存为“*.txt”文档。

1.3 图像分析

经过图像处理后提取到的孔隙面积、孔隙周长，按照以下公式[9]计算土样横断面孔隙的成圆率：

C=4 πA/P2

(1)

式中：C为成圆率，其值介于0～1之间，当C=1时表示孔隙为圆形；A为孔隙面积，mm2；P为孔隙的周长，mm。

统计土样剖面孔隙面积的成圆率的最大值、最小值、平均值、标准差，分析它们在土壤中深度的变化以及在不同盐分水平之间的差异性。

孔隙数是土壤某层CT切片中孔隙的数目，记为n。

孔隙度的定义为单位容积土壤中孔隙容积所占的百分数如下：

f=Vp/V

(2)

式中：f为土壤的孔隙度，无单位；Vp为土壤孔隙的体积，mm3；V为整个土壤的体积，mm3。

Vp=Sph

(3)

V=Sh

(4)

式中：Sp、S分别表示某一层孔隙的面积和某一层整个截面的面积，mm2。

将公式(3)、(4)代入式(2)中，消去相同的层厚h=0.63 mm，则有：

f=Sp/S

(5)

因此，如式(5)所示，计算土壤的孔隙度则是计算孔隙的面积占整个土样剖面面积百分数；分析的剖面取半径为20 mm的圆形区域，因此整个剖面的面积为1 256.64 mm2。

2 结果与分析

2.1 土壤孔隙成圆率的变化

土壤孔隙的成圆率(C)是表征土壤孔隙的形态特征，成圆率值介于0～1之间，其值越接近于1说明土壤的剖面的孔隙越接近圆形。一般认为，成圆率受到土壤孔隙形状和孔隙大小的影响，较小的孔隙其成圆率一般更接近1，因此根据成圆率可以判断土壤孔隙的大小。因考虑表层和底部的扰动影响，从表层和底部各去掉了5张切片，故实际分析范围为3.15～47.25 mm。本次试验的6个土柱不同深度的成圆率分布如图2所示。总体而言，6个土柱的成圆率值均介于0.5～0.75之间，并且没有表现出随土壤深度明显的规律性变化，都在0.5～0.75之间波动。

图2 土壤孔隙成圆率随深化度的变化

此外，盐分影响土壤的成圆率。具体而言，S6土样的盐分含量最高(ECe=54.92 dS/m)其各层的成圆率也最大，并且成圆率的最大值0.742 5也在S6土样中出现；S1土样的盐分含量最低(ECe=1.51 dS/m)，其平均成圆率也最小。此外，尽管中间盐分水平的成圆率没有呈现出明显的规律性的变化，但总体而言，除S5土样，土壤成圆率随盐分含量的变化(S2～S4)不断增大，更接近S1。

2.2 土壤孔隙度的变化

图3给出的是6个不同盐分土柱的孔隙度f随土层深度的变化，其对应的描述性统计分析参数和各层的平均孔隙度值如表2所示。具体而言，土壤孔隙度在土层深度上的变化是中间(约15～35 mm)比较大，两头值较小，而35 mm以后各处理的孔隙度也明显小于上层土壤的，这在表2中也可以体现。此外，与成圆率规律类似，盐分最高的S6处理(ECe=54.92 dS/m)孔隙度值最小(0.007)并且与其他盐分处理的孔隙度差异明显。同时，根据图3和表2中给出的标准差0.005 9的值，发现，S6的波动是最小的，S1处理(ECe=1.51 dS/m)的平均孔隙度最大且孔隙度值在整个深度范围内最不稳定，波动最大。除S5外，孔隙度随盐分的升高而表现出减小的趋势，孔隙度的最大值(0.258 9)和最小值(0.000 1)也分别出现在S1处理和S6处理中。

图3 土壤孔隙度随深度的变化

表2 土壤孔隙度描述性统计

2.3 土壤孔隙数的变化

图4给出的是土壤孔隙数随深度的变化图，从图4中可以看出，土壤孔隙数在土层中的变化存在一定的差异，约在35 mm后都表现为下降的趋势，深层(约40 mm以后)的土壤孔隙数明显少于土壤表层和中间土层的土壤孔隙数，因此，从整体趋势来看，随着土层深度的增加孔隙数是有先增大后减小的趋势，位于土层表面的孔隙数较土壤深层的孔隙数多，而中间(约15～35 mm)深度最多，各土层深度范围内的孔隙数变化都比较集中出现在20～60区间，并且变化差异较大，除S2(ECe=4.33 dS/m)处理出现81的高孔隙数值外，其他均未超过80的孔隙数值。考虑不同盐分的影响时，图像最左侧高盐的S6(ECe=54.92 dS/m)处理的孔隙数最少，这与孔隙度和成圆率的变化规律都是一致的，因此可以认为高盐分的土壤孔隙较少，尤其不利于大孔隙的形成。

图4 土壤孔隙数随深度的变化

进一步将土壤以层为单位分析每层土壤孔隙的数量以及整个土柱的孔隙数量，其结果如图5所示。整体而言，S2(ECe=4.33 dS/m)处理的孔隙数量最多，为3 458，S4(ECe=21.31 dS/m)处理、S5(ECe=28.89 dS/m)处理的总孔隙数紧随其后，S6(ECe=54.92 dS/m)的孔隙数最少[图5(a)]。观察图5～(f)也可发现，表层土壤(3.15～10 mm)和深层土壤(40～47.25 mm)的孔隙数较中间的少。而不同土层的孔隙数量分析也表明，除了表层3.15～10 mm的S1土样，其余各层的孔隙数量随盐分的变化规律均与全剖面的规律一致。而S1土样3.15～10 mm土层中孔隙数较大可能是由于取样受到扰动所致。

3 结 语

(1)CT扫描技术能够非破坏性地获取土壤结构特征，与其他方法相比，CT扫描方便快速并且能够准确描述孔隙的数目、大小等几何参数。

(2)土壤孔隙结构的形成受到盐分的影响。具体而言，当土壤盐分较高S6(ECe=54.92 dS/m)时，土壤孔隙的成圆率明显增加，土壤孔隙度明显降低并且孔隙度的波动随盐分增加逐渐减小，而盐分最低的S1(ECe=1.51 dS/m)处理，其成圆率明显较小、孔隙度较大。

(3)所取土柱的总孔隙数以及层间孔隙数随着盐分水平的提高表现为先增加后减小的趋势，较高盐分条件下(S4、S5)的土壤孔隙数较大，高盐条件下(S6)的土壤孔隙数最少。

□

[1] 孙 梅, 黄运湘, 孙 楠, 等. 农田土壤孔隙及其影响因素研究进展[J]. 土壤通报, 2015,(1):233-238.

[2] 王震宇, 于晓冬, 许 颖, 等. 土壤微孔对有机物吸附/解吸的影响及其表征[J]. 生态学报, 2009,(4):2 087-2 096.

[3] Marshall T J. A relation between permeability and size distribution of pores[J]. Journal of Soil Science, 1958,9(1):1-8.

[4] Belnap Jayne, Büdel B, Lange Otto L. Biological soil crusts: characteristics and distribution[M]. Springer, 2001.

[5] Fox D M, Bryan R B, Fox C A. Changes in pore characteristics with depth for structural crusts[J]. Geoderma, 2004,120(1):109-120.

[6] 郭 磊. 降雨过程中典型红壤孔隙时空变化规律研究[D]. 武汉：华中农业大学, 2007.

[7] 杨永辉, 武继承, 毛永萍, 等. 利用计算机断层扫描技术研究土壤改良措施下土壤孔隙[J]. 农业工程学报, 2013,(23):99-108.

[8] 程亚南, 刘建立, 张佳宝. 土壤孔隙结构定量化研究进展[J]. 土壤通报, 2012,(4):988-994.

[9] 李德成, Velde B, 张桃林. 利用土壤的序列数字图像技术研究孔隙小尺度特征[J]. 土壤学报, 2003,(4):524-528.

[10] 周 虎, 李文昭, 张中彬, 等. 利用X射线CT研究多尺度土壤结构[J]. 土壤学报, 2013,(6):1 226-1 230.

[11] Munkholm Lars J, Heck Richard J, Deen Bill. Soil pore characteristics assessed from X-ray micro-CT derived images and correlations to soil friability[J]. Geoderma, 2012,181:22-29.

[12] Hachicha M, Mansour M, Rejeb S, et al. Applied research for the utilization of Brackish/Saline water in center of Tunisia: water use[C]∥ Salinity evolution and crop response: Proceedings of International Salinity Forum, Riverside, 2005.

[13] Kanzari Sabri, Hachicha Mohamed, Bouhlila Rachida, et al. Simulation of water and salts dynamics in bouhajla (Central Tunisia): Exceptional rainfall effect[J]. Soil and Water Research, 2012,7(1):36-44.

[14] Department Natural Resources Management And. CHAPTER 7-SALTY SOILS[EB/OL]. http:∥www.fao.org/docrep/R4082E/r4082e08.htm.,2016-04-17.