土壤孔隙是评价土壤质量的重要指标之一。结构良好的土壤通常包含两种孔隙类型：一是由气候、田间管理和生物活动等因素形成的团聚体间的结构孔隙，二是由土壤颗粒形成的微团聚体内的质地孔隙[1-2] (图 1)。土壤孔隙结构是孔隙形态大小、数量搭配和空间分布状况的综合反映，由孔隙的几何特征(孔隙体积、表面积、分形维数等)和拓扑特征(连通性、弯曲度和欧拉数等)组成[3-4]，两者共同决定了土壤水分迁移、气体扩散、生物活动和植物有效水的获取等[5-6]。因此，开展孔隙结构的量化研究对于预测土壤水力性质、控制地下水污染、增加土壤肥力和作物产量等方面具有重要意义。

当前，量化孔隙结构的方法包括间接法和直接法，间接法以水气信息的变化间接推测孔隙的分布特征，如压汞法[7]、土壤水分特征曲线法[8]和气体吸附法[9]等，能粗略获取孔隙度和孔隙大小分布信息，但无法衡量孔隙的形态参数和空间构型[10-11]。直接法以数字图像为基础，经图像处理实现了孔隙结构的直接可视化和定量化[12]。如土壤切片法[13]、扫描电镜法[14]和染色示踪法[15]等，可满足孔隙二维特征参数的提取，但却具有制样麻烦，操作费时费力，对孔隙结构破坏性大的缺点。计算机断层扫描(computed tomography，CT)技术的发展，实现了原状孔隙结构的三维构建及其拓扑结构(连通度、弯曲度、欧拉数等)的研究，具有高效、快速、无损的特点，但易受样本大小和空间分辨率之间权衡的影响，主要用于小体积样本研究[6，16-17]。

因此，开展孔隙研究方法的评述，明确已有研究方法的不足，对于孔隙研究方法的创新及其主导的生态过程的研究具有重要意义。本文基于现有孔隙研究方法的发展历程，综述了常用孔隙定量研究方法(间接法和直接法)的基本原理、优缺点及不足，并对未来研究方法的发展方向进行了展望，以期为孔隙结构量化、重建和相关土壤过程的研究提供方法和理论支撑。

土壤取样方法因研究目的和研究方法而异，主要分为原状土样和混合土样采集两大类。原状土样可通过环刀、铝盒或PVC管等工具采集，采集后的土样用保鲜膜和泡沫封装好，防止水分蒸发和运输过程中的剧烈扰动，以备后续孔隙结构研究[4，13]。混合土样用土钻法、铁铲法采集[18]。土钻法先将地表植物小心去除，然后以重锤少击的方式将土钻按压进土壤(注意锤击次数，避免过击造成土壤压实)，多点采集后，混合样本，并轻置于塑料密封袋或铝制容器中[19]。铁铲取样同土钻法类似，采集的块状土壤需轻轻剥离与工具接触变形的部分，再将其密封保存。采集运回的土样根据土壤的塑性大小，于室内风干过程中将其沿自然结构断面轻轻剥离成直径小于10 mm的小土块，并剔除里面的根系和石子，风干后，干筛成不同粒径的团聚体，以便后续团聚体孔隙研究[20]。若需挖掘土壤剖面，需划分土壤层次取样，获取剖面图像则尽可能选择高分辨率相机。

土壤环境的复杂性和不透明性，使早期的孔隙结构量化大多建立在假设的基础上，通过水、气信息的变化间接对孔隙结构进行判断[10]，如表 1所示。容重法是应用最早的孔隙研究方法，简单易行，但只能获取总孔隙度信息。压汞法因其孔隙范围测定广，常与其他研究方法结合使用，如Dal Ferro等[22]结合压汞法和微CT技术量化了不同施肥条件下的团聚体孔隙结构特征。汞为非浸润液体，需施加一定压力才能进入孔隙空间，根据Young-Laplace方程将压力转换为等效孔径，即可获得孔隙大小分布，但汞无法进入闭合孔隙，易导致孔隙数量被低估，加上“ink-bottle”效应，其测定的孔径往往小于实际孔径[7]。孔隙大小分布也可通过土壤水分特征曲线法测定，但其测定时间往往比压汞法长很多[21]。气体吸附法利用N2、CO2和水蒸气等气体分子的吸附量推导孔隙的分布特征，主要用于毫米级土壤样品分析，测定的孔径范围很小，一般在1 ~ 200 nm[23]，不能用于孔隙结构的多尺度研究。总之，间接法只能粗略地获取孔隙分布信息，无法衡量孔隙的形态参数，加上其破坏性强，测定的结果往往与实际孔隙具有较大差异。

数字图像技术的发展，使孔隙结构的研究由基于水气信息的间接推测过渡到基于图像分析的直接量化[22]。如表 2所示，直接量化一般通过专业设备获取土壤结构图像，再经过图像处理和量化分析，即可实现孔隙结构的2D、3D可视化以及孔隙结构特征值的提取，此过程明显降低了人为操作对孔隙的直接破坏，更逼近土壤结构的真实空间分布状况。挖掘拍照法是应用最早的直接量化方法，通过挖掘土壤剖面获取土壤图像，虽操作简单，但挖掘过程对孔隙结构扰动大，加上土壤基质与孔隙的对比差异小，易造成孔隙数量被低估[24]。扫描电镜法可用于土体表面精细结构的研究，但样品在扫描前需经过干燥、喷金处理，以增加其导电性，此过程会不同程度地破坏土体结构，加上设备自身限制，扫描电镜只能用于小样本或大样本的局部观测，导致测定结果的代表性受限[14，26]。用光学显微镜观察土壤切片之前，样品需经过干燥、树脂硬化、切片、抛光等处理，制备过程繁琐，且树脂不能充分填充孔隙，易造成孔隙被低估[13]。染色示踪法与拍照挖掘法类似，但需用染料提前着色样地，浸染一定时间后再挖掘着色剖面并拍照，经图像处理获得孔隙结构特征及量化信息，主要用于大孔隙研究[15]。上述方法主要用于孔隙2D特征的量化，难以满足孔隙立体空间构型信息的获取。CT技术和核磁共振成像(magnetic resonance imaging，MRI)技术的发展，直接实现了原状孔隙的3D可视化和定量化，使结果具有较高的鲁棒性[4，27]。相比于MRI技术，CT技术具有较高的空间分辨率，能获取更多的孔隙细节信息，但用于大样本检测时易受样本大小和空间分辨率权衡的影响，形成的图像对比度低[4]。相反，MRI技术在检测大样本时具有较高的对比度，能获取更多的孔隙，但易受土壤磁性粒子影响，不具备普遍适用性，可作为CT技术的补充研究方法以扩大孔隙的研究范围[25]。在孔隙研究中，CT技术是当前最广泛使用的方法，主要通过异位研究获取原状孔隙结构特征，但田间条件下的孔隙受诸多因素影响，存在时间和空间上的变异，具破坏性取样特点的CT技术仍然难以实现原位孔隙的长期观测，因此，如何可视化和定量化原位孔隙结构的动态演变过程将是未来研究的重点和难点。

除图像本身的分辨率、对比度影响之外，分割方法的选择也直接关系到量化数据的准确性和有效性。由于土壤图像背景复杂，包含的噪声较多，量化孔隙结构之前需要先进行目标对象(土壤孔隙)的提取，即图像分割。常用的图像分割算法有全局阈值分割算法、局部阈值分割算法、边缘分割算法、聚类分割算法和深度学习等[28-30]。

全局阈值分割算法是一种简单有效、应用广泛的分割方法。它依据图像的灰度特征来确定一个或多个阈值，将图像中每个像素的灰度值与阈值进行比较，根据比较结果将像素的灰度值分配到合适的类别中，适用于目标与背景灰度值具有明显差别的图像[31-33]。但该方法只考虑像素点灰度值本身的特征，对噪声敏感，容易高估或低估孔隙结构，鲁棒性不高[34-35]。常用的全局阈值分割算法包括迭代最佳阈值法和Otsu法等。

迭代最佳阈值法从设定初始阈值开始，按照一定方式对其进行改进，直到满足设定的标准，迭代结束，所得数值即为完成孔隙分割的最佳阈值[36]。该方法实现了自动选取图像阈值，原理简单，易于实现，对噪声具有一定的鲁棒性，避免了人为确定的主观性误差，能实现图像的批量处理。

Otsu法(大津法)，为迭代法的一种，通过图像灰度特性自适应地选取阈值，进而分割图像中的目标和背景[37-38]。其可在一定程度上避免过分割现象，但只适用于单一目标分割，或感兴趣目标属于同一灰度范围，若目标灰度范围过大，可能会造成部分孔隙信息的丢失。

局部阈值分割算法是以直接寻找区域为基础的分割技术，具有两种基本形式：一是区域生长，从单个像素出发，逐步合成并形成需要的分割区域；二是从全局出发，逐步切割至需要的区域。但该方法对噪声敏感，具有过分割的缺点。主要包括区域生长法、分水岭法等。

区域生长法[39]是从预先选定的一组代表不同区域的种子像素开始，根据制定的标准将具有相似灰度、强度、纹理等特征的相邻像素合并到种子像素所代表的区域，将新添加的像素又作为新的种子像素，继续对各像素遍历直到找不到符合标准的像素为止，即可完成孔隙结构的分割判断[40-41]。一些研究人员利用该方法实现了土壤大孔隙的分割，在分割过程中，初始种子点(即像素点)的位置和数量对大孔隙分割精度具有较大影响[42]。对于具有多个分割目标的图像的应用具有局限性，易造成孔隙过分割。

分水岭法是一种基于拓扑理论的数学形态学方法，通过寻找像素灰度值分布的局部极小值来确定分类阈值，以此完成图像的分割[1，43]。该方法对尺寸较小、颜色较浅的大孔隙具有良好的分割效果[44]。但在CT图像中，由于土壤基质中含有的砾石等杂质也会产生细微的灰度边界变化，从而可能导致孔隙结构的过分割现象[25，45]。

边缘分割算法通过获取不同区域的轮廓线以达到分割的目的，可通过微分算子实现。常用的一阶微分算子有Roberts、Prewitt、Sobel等[46]。一阶微分算子用梯度和方向2个模板对像素点进行卷积运算，以最大值作为输出实现边缘点提取。但因为其基于滤波器进行检测，会使检测能力受到图像特性和噪声的影响[36]。LoG算子是二阶微分算子，其采用Laplace算子计算高斯函数的二阶导数，将二阶导数的过零点作为边缘，但在边缘检测和噪声抑制方面存在矛盾[36]。Canny算子为高斯函数的一阶导数，在滤波的基础上增加了非极大值抑制和双阈值处理，通过对信噪比与定位乘积得到最优化逼近算子，实现孔隙边缘精准检测[47-48]，是一种具有滤波、增强、检测的多级边缘检测算法。

聚类分割算法是将图像空间中的像素用对应的特征空间点表示，根据它们在特征空间的聚类进行分割，然后将其映射回原图像空间，得到分割结果[30]。该算法提高了孔隙分割精度，但其运算速度和稳定性易受初始条件影响，且无法满足孔隙细节信息的获取。主要包括K-均值聚类算法、模糊C均值聚类算法(fuzzy C-means clustering algorithm，FCM)和其他改进算法(如快速模糊C均值聚类算法(FFCM))等。

K-均值聚类算法是根据数据到中心点的距离来确定数据的类别。将有N个数据的集合A分为K类，即K个簇，将簇中对象的均值表示为每个簇的中心点，通过距离计算公式计算出数据到中心点的距离，根据此距离将其归类到合适的聚类中[30，49]。聚类结束时检查聚类中心是否收敛，收敛则终止，否则就继续迭代。该方法是一种全局方法，适用于灰度区域不明显的图像，但算法具有不稳定性，易受噪声和聚类个数的影响。

FCM算法是最常用的模糊聚类算法[50-51]。它通过对目标函数的迭代优化计算图像上的像素点在各个分类中所占比率，实现像素点的初始模糊分类，对于不同的图像具有较强的鲁棒性[36，52]。FCM算法整体分割效果较好，常用于形状不规则、边界不明确的大孔隙分割，具有无监督、软分割等优点[53-54]。但每次迭代计算都需要考虑所有的像素点，耗时较长，在孔隙密集区域会出现孔隙结构相连的情况，且其分割精度、运算效率和稳定性易受初始条件(聚类数目、聚类中心等)的影响[55-56]。FFCM算法为FCM算法的改进，使用灰度级参与迭代运算，显著提高了运算速度[57]。但在描述复杂的孔隙结构时容易忽略细节信息，对于孔隙边界模糊的图像仍无法准确描述。

基于人工神经网络的深度学习是机器学习的一部分，它可以是有监督的、半监督的或无监督的。深度神经网络已应用于计算机视觉、医学图像处理、根系和孔隙的分割[28，58-59]。目前，应用于土壤领域主要有卷积神经网络(convolutional neural networks，CNN)和全卷积网络(fully convolutional networks，FCN)等。

CNN采用误差梯度反馈算法进行模型优化，通过模拟人类视觉进行逐层抽象以获取高级的图像语义信息，从而提高了分割性能。其在纳米尺度上已实现对土壤离子束扫描电镜(FIB-SEM)图像的分割[60]。但CNN模型参数复杂，易受全连接层区域尺寸和输出特征的限制，导致运行内存浪费[29，34]。FCN为CNN的改进算法，可以解决CNN的不足，其本质是将传统CNN的全连接层替换为全卷积层，以保证卷积特征能恢复为原始尺寸的二维矩阵，通过图像端对端的输出，实现任意尺寸的图像的输入都能输出相同尺寸的孔隙二值图像[61]。其包含卷积层[62]、池化层[63]和上采样层3个不同的采样层。卷积层通过对输入图像的卷积运算，可以提取孔隙结构不同层次的特征。为避免过多的网络参数造成过拟合现象，在保证特征的基础上，需要经过池化层将图像划分成小于原始尺寸的固定矩形区域(m×n)，通过平均值池化或最大值池化，完成对卷积特征的采集。上采样层是池化过程的逆运算，目的是重构原始尺寸图像，通过对不同层次的二维孔隙特征图中的像素级进行分类，完成孔隙结构分割。FCN有效利用了图像的多层次特征(如线条、形状和纹理等)，既能精确定位不规则的孔隙空间位置，又能准确分割细小的孔隙结构，实现孔隙细节信息的捕获，避免了噪声的影响，保证了图像分割的有效性。虽然深度学习在土壤孔隙、植物根系的分割提取中具有较好的输出结果，但也存在大样本依赖、易过拟合、缺乏可解释性等缺点。

数字图像的分割方法有很多，目前没有针对于土壤图像的特定分割方法。各类分割方法在应用中不断改进，常依据图像本身的特点(分辨率、类型、数量等)进行选择。阈值法仍然是当前最广泛使用的方法，部分已集成到专业的分析软件中，如Image J、VGStudio MAX等，主要应用于CT图像处理。在各类方法中，阈值法、边缘法和聚类法易受噪声影响，难以获取孔隙的细节信息，存在过分割和欠分割现象。深度学习在土壤领域中应用较晚，它能有效地避免噪声干扰，具有较好的输出结果，但其大样本依赖特性使其在实际应用中受限。因此，有必要建立公共土壤图像数据库，实现数据共享，为深度学习在土壤领域的应用提供更广阔的空间。

结构量化是孔隙研究的关键，分割后的二值图像可以直观反映孔隙的几何形态和空间分布状况。由于图像以矩阵的形式存储于计算机之中，孔隙信息提取、量化分析的实质就是计算机运用数学的方法分类、统计矩阵中不同的像素点。二值化图像的灰度值仅有两种状态：一是像素值为0的黑色区域代表孔隙，二是像素值为1的白色区域代表固相基质[36]，通过黑色区域像素值的计算，即可获得孔隙数目、孔隙面积、孔隙率、孔隙大小分布、孔隙成圆率和分形维数等参数，如表 3所示。特征值的大小可直接反映土壤结构状态，如孔隙数目越多，结构越规则，则土壤质地越松散，输水能力越强。

CT技术和MRI技术的发展，实现了土壤孔隙的三维可视化。部分三维指标由二维指标发展而来，其计算不同于二维指标像素值的计算，三维指标是通过体素值确定的，如孔隙体积由孔隙包含的体素之和确定，孔隙面积则定义为边界体素之和[67]。除此之外，三维结构还有一些特有指标，如欧拉数、Γ指标、弯曲度和各向异性等。

欧拉数是评估土壤孔隙连通性最常用的指标之一[68-69]，用于衡量一个结构在分裂成两个独立对象之前可以断开的连接数量。Euler-Poincaré公式如下：

式中：χ为三维对象的欧拉数；β0为孤立对象的数量；β1为冗余连接的数量；β2为闭合腔的数量。计算单个对象的欧拉数时，β0始终等于1，β1和β2的值决定单个对象的欧拉数[70]。β1的值越大，孔隙的分枝越多。而β2的数量相对较少[71]，多数情况下可以忽略[69，71]。因此，χ值越小，孔隙之间的连接越多，连通性越高，反之，孔隙连接越少，独立孔隙越多[72]。

除欧拉数外，还有一个无量纲Γ指标，Γ指标反映了穿过孔隙系统连续路径的概率，χ则反映不考虑路径长度的内部连接数量[73]。当大多数孔隙属于最大连接簇[12，69]时，Γ值就越大。

式中：Np是孔隙p的体素数量；Ni是所有孔隙簇的数量；nk是孔隙簇k中孔隙体素的数量；Γ反映了两个随机选择的孔隙体素属于同一孔隙簇并受到最大孔隙簇nk强烈影响的概率[37]。当所有孔隙体素都属于同一连接簇(linked cluster)时，其范围从0(未连接孔隙簇较多)到1。Γ指标与孔隙度的关系曲线反映了Γ指标随孔隙度的增加而增加，没有明显的渗流阈值[74]。该曲线的斜率反映了相关孔隙簇大小的变化，即数据集样本之间渗流变化的度量[74]。

孔隙弯曲度表示孔隙在某一有效垂直长度上的复杂性，定义为实际孔隙长度除以欧氏长度(Euclidean length)[75-76]，计算公式如下：

式中：T为孔隙弯曲度；Ls为孔隙路径长度之和；L为孔隙垂直长度。其中，T是大于1的无量纲参数。其值越大，表明孔隙流动性越差；其值越小，表明孔隙的分枝越少，越有利于土壤水分流动。

各向异性(A)定义为孔隙投影宽度与投影高度之比[77]，关系式如下：

式中：A为各向异性；Lh为投影宽度；Lv为投影高度。各向异性对传输过程具有重要影响，A值接近1时，土壤常由具有立方或球形孔隙的各向同性介质组成。A > 1时，土壤孔隙多为水平孔隙，反映土壤容易被侵蚀；A < 1时，表明土壤孔隙具有一定方向，存在水平或垂直方向上的延长，低各向异性不利于孔隙水和溶质的保持。

土壤孔隙的研究已由破坏性取样–室内分析逐步过渡到原位监测–可视化研究阶段，研究尺度跨越米到纳米级，结合图像分析可实现孔隙2D、3D重构，能提取更多的孔隙细节信息。但土壤孔隙在时空尺度上受诸多因素的影响，具有较大的异质性，加之目前的研究都是基于小样本取样进行，结果的代表性往往较差。CT技术、MRI技术虽已实现原状孔隙的无损研究，但仍难以实现原位观测同一地点、同一孔隙的动态演变过程，且应用成本较高，很难满足大规模应用推广。因此，今后土壤孔隙定量研究尚有较大的突破空间，如：

1) 在图像获取上，亟需探索土壤孔隙原位动态观测系统。土壤孔隙结构及众多的生态过程在内外界因素作用下存在时间变异性和空间变异性。且目前针对土壤孔隙多属异位研究，仅限于样本的单一扫描。因此，亟需开发一种成本低、能定期自动扫描、可原位动态观测孔隙结构演变的无损检测系统。

2) 在图像处理中，建立各类土壤图像公共数据库，实现神经网络的无监督学习。在图像处理方面，基于神经网络的深度学习依赖大量的数据样本进行训练，但目前土壤孔隙图像的开放性共享资源非常稀缺，因此有必要建立土壤图像的公共数据库，实现数据共享和整合分析。加之高质量标注数据集的缺乏，且人工标注数据集成本较高，如何在不进行标注的情况下训练网络，实现无监督学习是未来研究的重要方面。

3) 多技术耦合，实现土壤多功能探究。如何关联孔隙结构量化数据与土壤生化、物理过程仍然是当前研究的重点和难点。三维X-CT与FIB-SEM结合，可使图像分辨率达到纳米级，与二维近红外光谱技术的结合，可实现在未扰动土壤样本中生物地球化学过程特性与土壤结构之间相互作用的评估。未来还需继续探索多技术之间的耦合，以土壤结构研究为纽带，实现土壤物理过程的原位观测和量化分析，同时链接土壤生物以及生源物质生物地球化学循环过程。