基于CT扫描技术的土壤孔隙定量表达优化

赵 玥 韩巧玲 赵燕东

(北京林业大学工学院，北京 100083)

基于CT扫描技术的土壤孔隙定量表达优化

赵 玥 韩巧玲 赵燕东

(北京林业大学工学院，北京 100083)

现有土壤孔隙量化方法主要通过图像处理软件实现孔隙结构的辨识与分析，此类通用的图像处理软件或医学处理软件未考虑土壤内部物质的复杂多变性以及孔隙结构的细小和不规则性，从而导致孔隙分割精度低进而量化误差大，为解决这一问题，本文针对土壤CT图像的特点提出了一种孔隙量化方法。该方法主要包括图像处理和量化分析两部分：选用自适应中值滤波算法去除噪声对孔隙边缘的影响，并采用迭代最佳阈值法与Canny边缘检测算子相结合的方法，准确识别出土壤孔隙结构及轮廓线；运用数学统计方法定量研究土壤孔隙率、孔隙数目、分形维数、成圆率等几何指标，用以揭示孔隙结构的复杂性和不规则性，实现对土壤孔隙的量化分析。最后，以冻融循环作用下的土壤为应用对象验证该方法性能。结果表明，本文方法能精确地定位孔隙轮廓，有效地分割孔隙结构，而且通过多种孔隙几何指标的量化可揭示出冻融循环作用对土壤结构的影响，为孔隙几何特征和空间特征的量化表达奠定了基础。

土壤断层扫描图像； 孔隙结构； 图像处理技术； 孔隙量化

引言

土壤孔隙是指土壤颗粒之间、团聚体之间或团聚体内部的空隙。孔隙结构包括孔隙数目、大小等几何形态，其特征决定了土壤水分保持能力与传导能力，对土壤水分和养分迁移等生态过程产生重大影响[1-4]。因此，对土壤孔隙结构进行量化分析，是从根本上认识土壤内部结构的必要前提。

CT扫描技术是孔隙量化研究的有效技术手段[5-8]，而孔隙的辨识与分割则是其重要前提。目前，研究者主要采用图像处理软件完成土壤孔隙结构的辨识，如文献[9-11]等利用与CT机配套的医学图像处理软件获得土壤孔隙结构，文献[12-14]采用Image J和Photoshop软件实现土壤CT图像的二值化，获得孔隙结构图像。其中，医学图像处理软件集成了针对人体骨骼和器官特点的处理程序，适用于识别具有较大连通域的病灶；而公共图像处理软件主要针对前景和背景具有较大差异的图像，实现目标的识别。而土壤CT图像中孔隙细小不规则(孔径多为毫米级)，与固相物质的灰度值差异也较小，受部分容积效应的影响在图像不同区域的灰度值不等，因此，采用现有图像处理软件实现土壤CT图像中孔隙的精确辨识具有较大难度，容易错误判别部分区域的孔隙结构，导致出现孔隙过分割和欠分割现象，使得实验结果无法真实反映孔隙的性质，从而对后续孔隙量化研究产生不可避免的影响。

为解决上述缺陷，本文将通过理论分析与实验比较两方面的研究，探讨现有的图像处理手段在土壤CT图像中的表现，并结合孔隙率、孔隙分形维数等量化指标提出一种针对土壤CT图像的孔隙量化方法，该方法针对土壤CT图像前景孔隙细小、连通性低与背景固相物质差异不一致等特点完成孔隙结构的高效分离，将为土壤孔隙的量化研究奠定技术基础。为验证所选方法的有效性，以冻融循环土壤为应用对象，通过孔隙量化结果的分析，揭示出土壤冻融循环对孔隙结构产生的影响。

1 孔隙量化系统

高精度土壤CT图像是获取孔隙几何特征和分布状态的数据基础，由于人眼对于灰度值的不敏感性，仅通过目视法无法得到孔隙的有效信息；而且原始土壤CT图像在获取和传输过程中易受到噪声、伪影等干扰，无法直接用于孔隙结构的量化研究[15]。因此，借助计算机图像处理技术加深对CT图像深层信息的理解是获得精确实验结果的必要前提。

基于上述原因，本文提出一种针对土壤CT图像的孔隙量化系统，用以深入探索孔隙几何特性。该系统主要分为2层：第1层为土壤CT图像的处理，首先扫描土壤样本获取土壤CT图像，这是提取孔隙有效信息的前提；其次，针对土壤CT图像特征进行相应的处理，主要选用自适应中值滤波、迭代最佳阈值和Canny边缘检测等算法实现，从而得到用于孔隙量化分析的土壤二值图像，为后续孔隙量化分析提供精准的数据。第2层为孔隙参数的量化，主要运用数学统计的方法获得土壤孔隙率、孔隙数目、成圆率、分形维数等指标，完成对土壤孔隙的量化分析。两层是逐渐递进的，前者是后者的基础，后者是前者的目的。如图1所示。

图1 孔隙量化系统流程图Fig.1 Flow chart of pore quantitative system

1.1 土壤CT图像处理

由于CT机的技术限制，原始土壤CT图像边界处易产生畸变，因此，需对其进行剪裁、存储等预处理，以提取用于计算机处理分析的土壤有效区域。但预处理阶段并没有去除噪声和伪影对图像的干扰，仍需通过自适应中值滤波的方法来抑制无用信息，突出有效信息，以便于对土壤孔隙结构进行精确分割。本文孔隙结构分割主要涉及到选用迭代最佳阈值法将土壤固体物质与孔隙结构分离，并通过Canny边缘检测算子准确提取孔隙轮廓线，从而得到只有孔隙结构的二值图像，为后续孔隙参数的量化分析提供精确的数据基础。

1.1.1土壤CT图像预处理

原始土壤CT图像包含两部分信息：第一部分，如图2a周边黑色区域所示，为CT机器的扫描时间、电压、窗宽和窗位等基本信息。另一部分为图2a圆圈部分的土壤有效图像，直径约为100 mm。由于扫描基本信息与后续研究无关，且矩形图像更易于计算机的处理，本文基于圆的内切正方形算法对原始土壤CT图像进行剪裁处理，保留正方形目标区域(图2a正方形)，并存储为如图2b所示的易于计算机理解的图像格式，用于土壤孔隙的量化研究。

图2 土壤CT图像预处理Fig.2 Preprocessing of soil CT images

1.1.2自适应中值滤波

由于CT机本身放射束硬化、部分容积效应等的约束，导致土壤CT图像在成像、传输等过程中存在噪声，在图像上表现为与孔隙边缘有一定的相似性，会影响孔隙轮廓线的提取，因此，对图像进行滤波处理是准确获取孔隙结构的前提。

现有研究中，滤波方法主要有均值滤波、维纳滤波、中值滤波等，其中均值滤波是以邻域内均值代替原像素值，会模糊孔隙与固相物质的边界，破坏土壤图像中原始细节信息；而维纳滤波是基于最小均方误差准则的线性滤波方法，无法去除土壤CT图像中非平稳的随机噪声，因此，这两种方法均不适用于土壤CT图像的滤波处理。而中值滤波采用模板中各像素点的中值代替原始像素值，能有效消除孤立噪声点，较好保留了土壤大孔隙边缘[16]，但由于恒定模板使得滤波效果和细节信息保护相矛盾，因此，中值滤波会去除图像中土壤小孔隙，不适用于细小孔隙研究。在此基础上，本文采用自适应窗宽的滤波算法，实现在消除图像噪声的同时保留孔隙的边缘和细节信息。

自适应中值滤波采用可变化的滤波窗，通过判断滤窗中心像素是否为脉冲噪声，进行不同的输出操作，实现对图像的滤波处理。算法主要包括噪声检测、滤窗选择和滤波输出3部分。实现过程为：若一定大小的滤波窗口Sij内中值滤波值Zmed与该点像素值均介于最大滤波值Zmax与最小滤波值Zmin之间，则表明该像素点与中值点都不是噪声，直接输出该点灰度Zij，以最大程度保留原始信息；若该像素点是噪声，则选择滤窗内中值滤波值Zmed输出。在这一算法中，采用尺寸变化的滤窗是为了针对滤窗中不同类型的像素点分别进行处理，从而实现在保留孔隙细节信息的同时去除高低脉冲噪声的干扰。

1.1.3迭代最佳阈值法

图像滤波后，孔隙边界较为清晰，与固体物质的区别更为明显，便于孔隙的提取。研究者普遍通过手动设定固定阈值的方法提取土壤孔隙，如PIERRET等[17]、 MONGA等[18]均通过该方法对孔隙的大小、数量等信息进行量化分析。手动设置阈值的方法需要根据具体的图像设定相应阈值，使得图像处理的效率较低，也不具备普遍适用性。另外，该方法中阈值的选择是通过人眼观察决定的，而人眼对于灰度的敏感性要远远低于计算机[19]，因此，针对土壤CT图像完成阈值的自动选取，是实现对土壤孔隙信息自动判读的关键。考虑到算法的实时性与精准性，本文选用基于逼近思想的迭代最佳阈值法，实现计算机自动选择最佳阈值，以完成孔隙的精确分割，为孔隙轮廓线的提取奠定技术基础。

迭代最佳阈值法主要思想：设定初始阈值，然后按照一定方式对其进行改进，直到满足设定的准则，迭代结束，所得数值即为完成孔隙分割的最佳阈值。其具体实现过程为：

(1)由图像最大灰度与最小灰度求得初始阈值T0。

(2)根据初始阈值T0将图像分割为固体区域与孔隙区域，分别求得两区域的平均灰度T1和T2。

(3)通过计算步骤(2)中T1和T2均值的方式确定迭代阈值T。

(4)比较迭代阈值T与初始阈值T0的差值d是否在允许范围内，若是，则认为该迭代阈值可以将图像中固体颗粒与孔隙目标有效地分离；相反，则将迭代阈值T赋予初始阈值T0，继续上述步骤(2)和(3)的迭代计算，直到新旧阈值差满足迭代条件。

(5)最终迭代阈值则为图像进行二值化处理的最佳阈值。

迭代最佳阈值法原理简单，易于实现，能自动找到每幅图像的精确阈值，对噪声具有一定的鲁棒性，可以准确地将土壤孔隙结构分离，适用于大批量的土壤CT图像。

1.1.4Canny边缘检测

基于土壤CT图像，提取孔隙的轮廓线，有助于精确计算孔隙结构的周长、面积等参数。边缘检测算子主要为基于微分法的一阶微分和二阶微分算子，其中Robert、Prewitt、Sobel等一阶算子采用代表梯度和方向的2个模板对像素点进行卷积运算，将最大值作为边缘点输出，而Log等二阶算子将二阶导数的过零点判为边缘。Robert算子适用于边缘陡峭的低噪声图像，Prewitt和Sobel算子的平滑操作会模糊图像边缘，适用于图像噪声小的图像，Log算子在边缘检测和抑制噪声方面存在矛盾。而在土壤CT图像的应用中，需要一种既能有效抑制图像中的噪声，又能精确检测出孔隙细小边缘的算法。而Canny边缘检测算子根据对信噪比与定位乘积进行测度，得到最优化逼近算子，能满足孔隙轮廓线提取对定位精度和边缘响应的要求，实现对孔隙边缘的精准检测，是最适用于土壤CT图像特点的算子，其实现过程如下：

(1)利用高斯滤波器对原图像进行滤波。

(2)通过在邻域内求有限差分来计算图像各点的梯度幅值和梯度方向。

(3)对梯度幅值进行非极大值抑制，剔除非极大值像素点。

(4)设置高低阈值，并结合连接分析方法，确定图像边缘。

Canny算子通过设定高低阈值，能够避免边缘信息丢失或者伪边缘的情况，从而准确检测出孔隙边缘。综上所述，采用Canny算子进行孔隙边缘的检测，能够得到精确的孔隙图像，不仅为孔隙结构的深入研究提供技术指导，也为孔隙参数的精确量化奠定数据基础。

1.2 孔隙参数量化

二值化土壤CT图像显示出孔隙的几何形态及分布，但由于人眼观察获取到的孔隙信息有限，因此，需借助计算机对孔隙信息进行深入挖掘，从而实现参数的量化分析。由于图像在计算机中是以矩阵的形式存储的，基于二值化图像对孔隙参数进行量化分析的实质就是对矩阵中不同像素点运用数学的方法进行分类与统计，而孔隙率、孔隙数目、孔隙周长、孔隙面积、成圆率以及分形维数这6个指标分别从不同方面表现了土壤中孔隙的拓扑结构及分布，对土壤的物理性质具有重要的影响。

1.2.1孔隙率

土壤孔隙率是指土颗粒间的孔隙体积占土壤总体积的比率，表征了土壤中自由水和空气的体积，是评价土壤质量和农作物生长环境的重要指标。由于在图像中，二值化图像的灰度值只有2个状态：代表土壤孔隙的黑色区域(像素为0)和代表土壤固态物质的白色区域(像素为1)，因此，黑色区域(像素值为0)的像素点占总像素点的比重，即为该土壤图像的孔隙率。

1.2.2孔隙数量

孔隙数量即图像中所包括的孔隙数目，在二值图像上表现为像素值为0的集合个数。由于孔隙的大小、形态各异，每个孔隙所包含的像素点数目不同，因此，本文基于连通域的方法来计算孔隙的数量。其基本思想为若中心点像素值与周围区域的8个像素值相同，则称这些像素点是连通的，将其判断为一个连通域，图像中连通域的数量即为孔隙数量。

1.2.3孔隙成圆率

孔隙成圆率是指孔隙形态与标准圆形的接近程度，表示的是不规则孔隙的形态特征，其值越接近1，孔隙结构越接近标准圆形。孔隙的形态特征对土壤的水分运输能力和透气性能具有较大影响。由于成圆率无法直接计算，因此，本实验采用间接法来求得孔隙成圆率

C=4πA/P2×100%

(1)

式中A——孔隙面积，像素×像素

P——孔隙周长，像素

1.2.4孔隙分形维数

土壤孔隙是细小的不规则结构，要对其复杂程度进行准确的描述，找到定量的指标是必要的。孔隙分形维数描述的是土壤孔隙的自相似特性，是土壤孔隙不规则性的综合体现[20]，可用于表征土壤孔隙结构的复杂性。

由于土壤CT图像是由像素点组成的灰度图像，每个像素点都有相应的灰度值，表示为灰度函数Z=f(i,j)，在三维空间中，灰度函数可等价为一个曲面，因此，本文利用双毯法计算孔隙分形维数。双毯法通过求取以灰度曲面为中心的具有一定高度的立方体体积，来估计灰度曲面的分形维数，并且只适用于求取图像表面分形维数。其基本思想如下：

将土壤CT图像的灰度函数f(i,j)视为1个曲面，其中(i,j)为图像上的坐标，选定灰度曲面上的高度ε处，分别建立覆盖该曲面的毯子，上毯子表示为u0(i,j)，下毯子表示为b0(i,j)，则构成1个厚度为2ε的“毯子”，毯子的表面积为体积除以2ε。对于不同的高度ε，可分别计算出分形表面积。

初始情况下，令

u0(i,j)=b0(i,j)=f(i,j)

(2)

上下两张“毯子”分别按如下原则生长

(ε=1,2,…)

(3)

(ε=1,2,…)

(4)

则“毯子”的体积为

(5)

式中w——图像宽度，像素

l——图像长度，像素

对于二维曲面，当毯子高度由ε-1增至ε时，曲面表面积为

(6)

根据分形定义，分形表面积满足关系

A(ε)=Fε2-D

(7)

lgA(ε)=c1lgε+lgF

(8)

式中F——常数D——分形维数

c1——直线斜率

对于不同的ε，可以计算出不同的lgA(ε)，通过最小二乘法拟合直线斜率，可以计算出c1，由

c1=2-D

(9)

即可计算出分形维数D，作为表征土壤图像表面纹理不规则性和复杂度的参数。

2 实例应用与验证

为测试本文选用的基于CT扫描技术的孔隙量化方法在孔隙结构提取与孔隙参数量化上的应用效果，以对孔隙结构影响较为强烈的冻融循环土壤为研究对象进行结果分析，并通过不同方法的对比实验，证明本文方法的优越性。

本文图像为由飞利浦Brilliance 64排螺旋CT机扫描土壤样本所得的CT图像。扫描样本分别为经历0、1、3、6、9次冻融循环的土壤样本，每次冻融循环选取3个重复样本，单个样本可得220张扫描图像，单次冻融循环共有660幅扫描图像，因此，本实验的土壤图像数据库共包括3 300幅土壤CT图像。

根据图像中土壤有效面积的位置, 采用最大内切正方形的方法，将每张土壤CT图像剪切为211像素×211像素后进行实验。

2.1 土壤CT图像滤波

为避免土壤CT图像中的噪声干扰，精确后续的边缘化检测与定量分析，本实验采取常用的中值滤波、均值滤波和维纳滤波3种滤波方法对土壤CT图像进行滤波处理，与本文自适应中值滤波效果进行对比，并将实验结果直观地显示在图3中。由于篇幅有限，图3为随机选取土壤图像数据库中一幅典型图像。

图3 土壤CT图像的滤波处理效果Fig.3 Filtering effects of soil CT image

由图3可以看出，相较于图3a所示的原图，图3b和图3d所示的中值滤波算法和维纳滤波算法较好地保留了土壤孔隙的边缘，但丢失了图像的细节信息；图3c所示的均值滤波算法平滑了土壤孔隙的边缘，使图像整体信息变得模糊。上述3种滤波算法虽然能够有效去除噪声，但是对于孔隙边缘的检测和细节信息的保护不尽人意，不适用于土壤CT图像孔隙结构的精确提取。而本文的自适应中值滤波算法，在去除噪声污染的基础上，完整保留了土壤孔隙的弱边缘和细节信息，如图3e所示。

经实验证明，本文采用的自适应中值滤波算法，能最大程度地保证孔隙结构，有利于进行孔隙结构的分割处理。

2.2 土壤CT图像二值化

目前常用于土壤孔隙研究的方法为全局固定阈值法，该方法主要通过阈值的调节，选取使得孔隙直径与真值最相近的数值作为二值化的阈值，本文采用此方法提取出的土壤孔隙结构图如图4b所示。由图4b明显观察到提取的孔隙数目偏多，孔隙连通域明显偏大，这是因为图像不同区域的灰度差不同，适合特定区域的阈值并不适用于整幅图像，以此为依据的孔隙结构提取易出现过分割或欠分割的现象，从而使得后续孔隙参数量化分析存在偏差。

由上文可知，阈值选取不理想，易导致孔隙连通域失真，而本文选取的迭代最佳阈值法则克服了这一难题。该方法选取的阈值综合了图像整体信息，是像素点共同作用的结果，适用于描述土壤CT图像各方面信息的变化，采用该方法提取的土壤孔隙结构如图4c所示。通过与图4a对比可知，采用迭代最佳阈值法提取孔隙结构的形态、数目更加符合实际情况，分割效果更为理想。

图4 土壤CT图像二值化结果Fig.4 Soil CT image binarization results

图5 土壤CT图像边缘检测结果Fig.5 Soil CT image edge detection results

2.3 土壤CT图像边缘检测

对孔隙参数进行量化分析，不仅需要提取孔隙结构，还需明确孔隙轮廓线。本实验采用经典的Roberts、Prewitt、Sobel、Log、Canny 5种边缘检测算子对孔隙结构进行边缘检测，并随机选取一幅图像用于检测效果的比较分析，检测结果如图5所示。

由图5可以看出，相较于图5a，图5b和图5c所示的Roberts算子和Prewitt算子对孔隙边缘定位较准，但提取的孔隙轮廓线不具备完好的连通性。图5d所示的Sobel算子对孔隙边缘定位不准，得到的孔隙结构有一定程度的失真。图5e所示的Log算子平滑过程中使得边缘具有延展性，导致孔隙结构被明显放大，孔隙形态失真。上述4种边缘检测算子都能够提取孔隙轮廓线，但是对于孔隙边缘的定位和形态的判断存在一定缺陷，不符合精确提取土壤CT图像孔隙轮廓线的要求。而图5f所示的Canny边缘检测算子具有定位精准和单边缘响应强的优点，提取的孔隙边缘连续性和清晰度比较理想，还可以检测出图像中的弱边缘，更适合土壤CT图像的细节刻画。

经实验证明，Canny算子能够准确定位和检测土壤孔隙的轮廓线，其与迭代最佳阈值法的结合有利于提取孔隙的真实结构，为孔隙的量化研究提供精确的数据基础。

2.4 土壤孔隙量化分析

有研究表明，季节性冻融循环对土壤团聚体结构与物理性质的影响较为强烈，而团聚体结构与物理性质的改变直观地表现为土壤孔隙数量与形态的变化，因此，为探明冻融循环对土壤结构的影响，本文基于土壤CT图像，实现了对孔隙结构的量化研究。本文主要通过对孔隙率、孔隙数量、孔隙成圆率、分形维数等参数进行具体的分析，实现对土壤孔隙结构细致深入的研究，从而加强对土壤孔隙几何特性和空间特性的了解，进而判断冻融循环对农业系统、生态系统等方面的影响。

基于上述土壤CT图像数据库，本实验所得孔隙参数如表1所示。由于篇幅有限，部分参数以多幅图像均值形式表示。

表1 不同冻融次数孔隙参数的比较Tab.1 Comparison of pore parameters of different freezing-thawing cycles

由原始土壤CT图像可知，土壤孔隙表现为CT图像中不同范围的黑色区域，土壤孔隙率则体现为黑色区域所占图像比例，由表1比较得出，随着冻融循环次数的增加，孔隙率逐渐减小。孔隙率对冻融循环次数的响应主要受两方面影响：孔隙数量，如冻融循环0次所示，孔隙数量最大，说明黑色像素点覆盖的区域较多，像素点占图像的比例也就大，孔隙率也较高；孔隙面积，如冻融循环1次所示，虽然孔隙数量较冻融循环3次少，但由于孔隙平均面积大，单个孔隙所占像素数多，其孔隙率仍然较大。土壤孔隙的分形维数则反映了孔隙的复杂程度，其变化与冻融循环次数没有直接的正负相关关系，但其受孔隙率与孔隙成圆率综合因素的影响。如比较表1中冻融循环3次和6次的参数，虽然前者的孔隙率与孔隙数量都较大，但是由于其孔隙成圆率较小，两参数综合作用下，使得前者的孔隙分形维数略小于后者。冻融循环对孔隙数量、孔隙面积和孔隙成圆率虽然产生显著影响，但并没有明显的变化规律。综上所述，各参数对冻融循环次数不是单一响应的，而是通过互相影响与制约决定了土壤结构与物理性质的变化。

由表1数据可知，同一图像的孔隙参数在不同冻融次数间的差异较小，因此，研究冻融循环对土壤孔隙结构的影响需要保证孔隙量化过程中所产生的误差。而本文方法针对孔隙量化表达的所有环节，采取不同处理算法确保高获取数据的精度，从而得到较为真实的孔隙参数，完成土壤孔隙的定量表达。

3 讨论

本文基于CT扫描技术的孔隙量化系统综合了3种处理算法，以冻融循环土壤CT图像为研究对象进行不同算法实验结果的比较，从而选择适用于土壤孔隙量化表达的算法，增强了孔隙形状的还原度与孔隙轮廓的精确度。但由于算法的比较主要通过目视法定性分析实现，由此获得的孔隙参数未与真值进行对比，缺乏一定的说服力。

由于土壤孔隙结构的不规则性，目前尚未有获取孔隙真值的标准方法，难以准确获取土壤孔隙参数的真值。但已有研究人员采用不同方法进行孔隙结构的量化分析，因此，研究者可以通过与现有量化方法的对比实验证明所用方法在孔隙还原与轮廓提取方面的优势。

4 结束语

土壤孔隙状况是土壤的物理特性之一，它决定了土壤中空气和水分的含量，进而影响土壤各种变化进程及作物的生长，是目前评价土壤质量优劣的指标和研究热点之一。本研究针对土壤CT图像特点，提出了一种针对土壤CT图像处理的孔隙量化系统，该系统通过自适应中值滤波、迭代最佳阈值法与Canny边缘检测算子的结合，提取清晰的土壤孔隙轮廓线，增强了孔隙形状的还原度与孔隙轮廓的精确度。在此基础上，通过对像素值的统计分析，获得土壤孔隙率、孔隙数目、成圆率、孔隙分形维数等几何指标，为后续土壤孔隙连通性、土壤水分运移状况分析等研究奠定了数据基础。并经实验证明，本文提出的将数字图像处理技术与高精度土壤CT图像相结合的孔隙量化方法，能够有效解决无法准确区分具有相似性背景与目标的缺陷，较好地从土壤固相中分离出孔隙结构，准确识别和定位孔隙轮廓线，为土壤孔隙研究提供了一种有效技术手段，对土壤孔隙几何特征和空间特征的研究具有重要影响。

1 HILL R L, HORTON R, CRUSE R M.Tillage effects on soil water retention and pore size distribution of two mollisols [J].Soil Science Society of America Journal, 1984, 49(5):1264-1270.

2 张慧娟，孙宇瑞,林剑辉,等.不同粗糙度尺度下预测表层土壤孔隙率量化指数比较研究[J].应用基础与工程科学学报,2009,19(1):69-76.

ZHANG Huijuan, SUN Yurui, LIN Jianhui, et al.Comparative study on predicting porosity by different roughness indices [J].Journal of Basic Science and Engineering, 2009, 17(1):69-76.(in Chinese)

3 TOKUMOTO I, NOBORIO K, KOGA K.Coupled water and heat flow in a grass field with aggregated Andisol during soil-freezing periods[J].Cold Regions Science and Technology, 2010, 62(2-3):98-106.

4 STEPHENS P R, HEWITT A E, SPARLING G P, et al.Assessing sustainability of land management using a risk identification model[J].Pedosphere, 2003, 13:41-48.

5 PETROVIC A M, SIEBERT J E, RIEKE P E.Soil bulk density analysis in three dimensions by computed tomographic scanning [J].Soil Science Society, 1982, 46:445-450.

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

CHENG Ya’nan, LIU Jianli, ZHANG Jiabao.Advance in the study on quantification of soil pore structure [J].Chinese Journal of Soil Science, 2012(4):988-994.(in Chinese)

7 ELLIOT T R, HECK R J.A comparison of optical and X- ray CT technique for void analysis in soil thin section [J].Geoderma, 2007, 141(1-2): 60-70.

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

YANG Yonghui, WU Jicheng, MAO Yongping, et al.Using computed tomography scanning to study soil pores under different soil structure improvement measures [J].Transactions of the CSAE, 2013,29(23):99-108.(in Chinese)

9 WAMER G S, NIEBER J L, MOORE I D, et al.Characterizing macropores in soil by computed tomography [J].Soil Science Society of America Journal, 1989, 53(3): 653-660.

10 PEYLON R L, GANTZER C J, ANDERSON S H.Fractal dimension to describe soil macropores structure using X-ray computed tomography [J].Water Resources Research, 1994,30(3): 691-700.

11 CORTINA-JANUCHS M G, QUINTANILLA-DAMINGUEZ J, VEGA-CORONA A, et al.Detection of pore space in CT soil images using artificial neural networks [J].Biogeosciences, 2011, 8(2):279-288.

12 RANJITH P U, STEPHEN H A, CLARK J G, et al.Influence of prairie restoration on CT-measured soil pore characteristics [J].Journal of Environmental Quality, 2008, 37(1):219-228.

13 赵世伟，赵勇钢，吴金水.黄土高原植被演替下土壤孔隙的定量分析[J].中国科学：地球科学,2010，40(2)：223-231.

ZHAO Shiwei, ZHAO Yonggang, WU Jinshui.Quantitative analysis of soil pores under natural vegetation successions on the Loess Plateau [J].Science of China: Earth Science, 2010, 40(2):223-231.(in Chinese)

14 王恩姮,卢倩倩,陈祥伟.模拟冻融循环对黑土剖面大孔隙特征的影响[J].土壤学报,2014,51(3):490-496.

WANG Enheng, LU Qianqian, CHEN Xiangwei.Characterization of macro-pores in mollisol profile subjected to simulated freezing-thawing alternation [J].Journal of Soil Science, 2014, 51(3):490-496.(in Chinese)

15 王艳丽,程展林.CT扫描技术在我国土工试验中的应用研究进展[J].地震工程学报,2015,37(增刊1):35-39.

WANG Yanli, CHENG Zhanlin.Progress in the application of CT scanning technology in chinese soil tests [J].China Earthquake Engineering Journal, 2015,37(Supp.1):35-39.(in Chinese)

16 华珊,陈研,梁露焘,等.利用基于偏微分方程的图像滤波技术研究土壤孔隙结构[J].农业工程学报,2014,30(3):78-85.

HUA Shan, CHEN Yan, LIANG Lutao, et al.Studying soil pore structure by using image filtering technology based on partial differential equation model [J].Transactions of the CSAE, 2014, 30(3): 78-85.(in Chinese)

17 PIERRET A, CAPOWIEZ Y, BELZUNCES L, et al.3D reconstruction and quantification of macropores using computed tomography and image analysis [J].Geoderma, 2002, 106(3-4):247-271.

18 MONGA O, BOUSSOA M, GARNIERB P, et al.3D geometric structures and biological activity:application to microbial soil organic matter decomposition in pore space [J].Ecological Modelling, 2008, 216(3-4):291-302.

19 姬伟,陶云,赵德安,等.基于CLAHE的苹果树树枝迭代阈值分割方法研究[J/OL].农业机械学报，2014,45(4):69-75.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?file_no=20140411&flag=1.DOI：10.6041/j.issn.1000-1298.2014.04.011.

JI Wei, TAO Yun, ZHAO Dean, et al.Iterative threshold segmentation of apple branch images based on CLAHE [J/OL].Transactions of the Chinese Society for Agricultural Machinery, 2014,45(4):69-75.(in Chinese)

20 王聪颖,张慧娟,孙宇瑞,等.基于分形理论的土壤粗糙指数与孔隙率映射规律研究[J].农业机械学报,2011,42(11):32-38.

WANG Congying, ZHANG Huijuan, SUN Yurui, et al.Relationship between soil surface porosity and roughness indices based on fractal theory [J].Transactions of the Chinese Society for Agricultural Machinery, 2011,42(11):32-38.(in Chinese)

OptimizationofSoilPoreQuantitativeExpressionBasedonComputedTomographyScanningTechnology

ZHAO Yue HAN Qiaoling ZHAO Yandong

(SchoolofTechnology,BeijingForestryUniversity,Beijing100083,China)

In recent years, image processing software was wisely applied to identify and analyze pore structure.However, these softwares, such as Photoshop and Image J, did not take into account the complexity of the internal material in the soil and the irregularity of pore structure, and they caused low pore segmentation precision.In order to solve the problem, a pore quantitative method based on the characteristics of soil computed tomography (CT) image was proposed.This method mainly included image processing and quantification analysis.Firstly, the adaptive median filtering algorithm was adopted to remove the effect of image noise on the edge of pore.Then, the method of iterative optimal threshold and canny edge detection was used to identify the pore structure in the soil and the contour line of the pore.Secondly, soil pore structure had evident spatial characteristics, which included soil porosity, pore number, pore radius, spore size distribution, circularity, fractal dimension, and so on.These pore geometry indicators were calculated by using the mathematical statistics method, and they could reveal the complexity and irregularity of pore structure.These geometry indicators were useful for realizing the quantitative analysis of the soil porosity.Finally, the method was applied to the soil under freeze-thaw cycle.The results showed that the method can accurately locate the pore profile, and segment the pore structure effectively.Furthermore, the effect of freezing and thawing cycles on the soil structure was revealed by quantifying the geometrical parameters of various soil pores, thus it proved the effectiveness of the method and laied foundation for quantification of soil pore geometry and spatial characteristics.

soil computed tomography image; pore structure; image processing technology; quantitation of pore

10.6041/j.issn.1000-1298.2017.10.031

S152

A

1000-1298(2017)10-0252-08

2017-01-12

2017-02-14

中央高校基本科研业务费专项(BLX2015-36)和国家自然科学基金青年基金项目(41501283)

赵玥(1986—)，女，讲师，博士，主要从事图像处理与模式识别等研究，E-mail: zhaoyue0609@126.com

赵燕东(1965—)，女，教授，博士生导师，主要从事生态信息智能检测与控制研究，E-mail: yandongzh@bjfu.edu.cn