土壤结构是土壤中固体成分（矿物和有机质）和土壤孔隙的空间排列。定量化分析土壤结构对于了解土壤的物理、化学和生物过程有至关重要的作用^[1]。近年来，随着CT扫描技术的发展，实现了在不破坏土壤结构条件下实时获取土壤的孔隙结构^[2]。然而，由于土体内物理、化学和生物过程相互影响以及各种地质过程和人类活动的影响，土壤结构呈现不规则、不稳定和高度复杂性，需要采用合适的理论定量化描述土壤结构。因此，学者将分形理论引入土壤科学中，前期分形理论集中于质量和表面分形维数的研究^[3]。随着分形理论的发展，分形理论逐渐演变为多重分形分析。多重分形分析能够捕捉土壤孔隙结构的尺度规律特征并对微小的局部变化非常敏感，2003年Posadas等^[4]首次基于图像处理研究土壤孔隙结构的多重分形特征。目前，多重分形分析更多应用于二值图像的分析，对于灰度图像自身的研究很少。对于二值图像而言，阈值将影响其的多重分形参数^[5]，然而，目前图像阈值的准确确定仍然是一个悬而未决的问题^[6]。

此外，一些学者指出二值图像并非满足数学多重分形，并且仅在大尺度中存在多重分形特征^[7-8]。当计算二值图像的多重分形特征时，图像的大小限制了可用于研究的尺度。此时仅通过少量尺度获得的多重分形参数是不可靠的，必须使用其他的多尺度分析方法。Bird等^[7]认为在二值图像中应用相对熵的方法能够更准确地提取土壤孔隙尺度变化的信息。与分形方法不同，相对熵不需要任何其他先验假设，因此相对熵可以在所有的研究尺度上进行。该方法已经成功应用于不同的土壤样带数据^[9]，但在定量表征土壤孔隙结构方面鲜有应用。

近年来，生物质炭作为农作物秸秆等有机物质在缺氧或者限氧条件下热裂解得到的产物^[10]，已被广泛应用于土壤改良。研究发现，生物质炭能够提高土壤有机碳含量、改善土壤物理结构^[11]，然而生物质炭对于土壤孔隙结构的具体影响仍然存在争议。An等^[12]认为生物质炭处理后土壤的孔隙结构和水稻生产力显著提高，而Fan等^[13]发现生物质炭处理样品的孔隙率和大孔隙度最低，改良效果不明显。因此，需结合分形理论和图像处理技术从微观角度阐明生物质炭添加对土壤结构的影响。基于此，本研究以江苏沿海围垦区盐渍土为研究对象，基于多重分形理论和相对熵理论，结合Micro-CT扫描技术，对比分析这两种方法在灰度图像和二值图像中的应用前景，并揭示生物质炭施加后土壤孔隙结构的多重分形特征和相对熵特征的变化，从而为生物质炭改良土壤结构提供理论依据。

1 材料与方法 1.1 供试材料与试验设计

为了探究生物质炭对土壤结构的改善效果，在2016—2020年于河海大学江宁校区开展改良盐渍土水稻栽培测坑试验。试验测坑为圆柱形，上部直径为840 mm，高度为850 mm，试验坑的容量为300 L。试验在塑料大棚避雨条件下开展。测坑内装填土壤于2015年取自江苏省如东垦区，该地区土壤质地为粉砂壤土，其土壤基础性质见表 1。供试土壤按1.35 g·cm^–3的容重分层装入测坑。表层土壤（0~20 cm）采用生物质炭改良，设置0%、2%和5%（占表层土壤质量比）3个生物质炭添加梯度。本试验所采用的生物质炭是河南三利公司生产的小麦秸秆生物质炭，其基础性质见表 1。水稻于每年6月移植，10月收获。

考虑到第一年（2016年）土壤的沉降与固结，结构扫描试验从2017年开始。在水稻收割后，在每个测坑中心选择平整的土壤平面，用小铲子清除土壤表面的杂物，采用100 cm³的标准塑料环刀取表层原状土（0~20 cm）进行Micro-CT扫描。CT扫描试验在中国科学院南京土壤研究所进行，采用高分辨率X射线数字岩心分析系统（Phoenix Nanotom S，GE，USA）。Micro-CT扫描通过计算机控制射线源发出射线束，数控扫描平台承载塑料环刀所取的原状土，在计算机控制下以0.5$^\circ \cdot$s^–1的速度旋转，平板探测器负责采集扫描的数据，最后计算机通过采集到的投影数据重建土壤的CT切片图像。本研究中扫描X射线的最大能量为100 kV的电压和100 μA的电流，扫描精度为25 μm。

本试验采用ImageJ软件分析CT扫描图像。为了避免边缘效应并降低计算机处理时间，选择图像中心512×512×512 voxels部分作为研究区域。然后增加图像的亮度和对比度，以确保图像灰度值的范围为0~255，如图 1（a）所示。接着采用中值滤波（2×2×2 voxels）对灰度图像进行滤波处理，去除图像中的噪点。当灰度图像转换为二值图像时，阈值的确定决定了二值图像能否准确反映原始图像的信息^[14]。本研究通过以下方法确定图像阈值：（1）复制灰度图像并导入ImageJ软件；（2）采用软件默认方法手动调整阈值，对复制的图像进行二值化；（3）仔细比较原始灰度图像和二值图像，通过视觉判断确定适当的阈值，二值图像如图 1b所示。土壤CT孔隙度采用Volume Fraction Plugin计算。为获取土壤孔隙分布，采用LocalThickness插件分割二值图像，分割后孔隙重建图如图 1（c）所示。不同的颜色代表不同的孔径，颜色越亮表示孔径越大。根据分割后图像，采用Histogram得到土壤孔隙分布。

1.2 多重分形分析

多重分形分析是分析系统尺度不变性的标准工具^[15-16]。在本研究中，采用盒计数法计算三维图像的多重分形特征。首先，对于大小为$R \times R \times R$的图像，采用大小为r的盒子去覆盖图像，所需的盒子数量记为$N\left( r \right)$，每个盒子存在孔隙的概率密度函数记为${\mu _i}\left( r \right)$。对于非均匀系统，概率密度函数${\mu _i}\left( r \right)$和计算尺度$\varepsilon$存在以下关系：

${\mu _i}\left( r \right) \sim {r^\alpha }$

(1)

式中，$\alpha$为奇异性指数。此时，根据缩放关系定义$f\left( \alpha \right)$：

$N\left( {\alpha ,r} \right) \sim {r^{ - f\left( \alpha \right)}}$

(2)

式中，$f\left( \alpha \right)$反映点集的分形维数。$\alpha$与$f\left( \alpha \right)$的映射图称为多重分形谱。当多重分形谱显示出倒U的形状时，说明分析的对象显示出多重分形特征。此外，根据度量矩的标度律，也可判断研究对象是否存在多重分析特征，即定义配分函数$\chi \left( {r,q} \right)$在不同维度q下为：

$\chi \left( {r,q} \right) = \mathop \sum \limits_{i = 1}^{N\left( \varepsilon \right)} {\left( {{\mu _i}\left( r \right)} \right)^q}$

(3)

根据式（1）和式（3），配分函数$\chi \left( {r,q} \right)$与计算尺度$\varepsilon$满足：

$\chi \left( {r,q} \right){\text{ \sim }}{r^{ - \tau \left( q \right)}}$

(4)

$\tau \left( q \right) = - \mathop {\lim }\limits_{r \to 0} \frac{{\log \mathop \sum \nolimits_{i = 1}^{n\left( r \right)} {{({\mu _i}\left( r \right))}^q}}}{{\log \left( r \right)}}$

(5)

式中，$\tau \left( q \right)$为q维度的累积量生成函数。当$\chi \left( {r,q} \right)$与r的双对数曲线图满足线性关系时，则研究对象具有多重分形特征。

基于多重分形理论，为了定量化描述对象的多重分形特征，选择以下5个量化参数：（1）${\alpha _0}$；（2）${\alpha _1}$；（3）${\alpha _2}$；（4）振幅$\Delta \alpha$；（5）不对称性系数$\Delta f$。$\Delta \alpha$与$\Delta f$的计算公式如下：

$\Delta \alpha = {\alpha _{\max }} - {\alpha _{\min }}$

(6)

$\Delta f = f\left( {{\alpha _{\min }}} \right) - f\left( {{\alpha _{\max }}} \right)$

(7)

${\alpha _0}$、${\alpha _1}$和${\alpha _2}$与${D_0}$、${D_1}$和${D_2}$一一对应，分别反映了研究对象的容量维数、信息维数和关联维数；$\Delta \alpha$反映了对象的复杂程度，$\Delta \alpha$越大，对象的复杂程度越大；$\Delta f$反映了多重分形谱的形状，$\Delta f < 0$，多重分形谱呈现“左钩形状”，研究对象的大概率子集占比更高；$\Delta f > 0$，多重分形谱呈现“右钩形状”，研究对象的小概率子集占比更高。在本研究中，采用ImageJ中MULTIFRAC插件分析灰度图像和二值图像孔隙结构的多重分形特征^[17]，其中q的取值范围为[–5，5]，步长为0.1。

1.3 相对熵分析

熵的尺度特性是量化系统异质性的标准工具^[17]，其揭示系统的异质性如何随着尺度的增加而演变。Tarquis和Bird^[18]指出并非所有系统的熵会随着尺度产生明显的变化，为放大熵随着尺度的变化，提出了相对熵$E\left( r \right)$，相对熵$\left( {E\left( r \right)} \right)$可用式（8）表示：

$E\left( r \right) = - \mathop \sum \limits_i {\mu _i}\left( r \right){\text{log}}{\mu _i}\left( r \right) + \log \frac{r}{R}$

(8)

根据相对熵，计算不同尺度下的相对熵增量能够在较小的尺度上量化系统结构^[16]。本研究中，选择分辨率作为尺度衡量对象，当分辨率增加2倍时，相对熵增量$\Delta E$定义为：

$\Delta E = E\left( {2r} \right) - E\left( r \right)$

(9)

本研究采用以下5个参数量化相对熵结果：（1）$\Delta {E_{\max }}$；（2）$\Delta {E_{\min }}$；（3）$Range\Delta E\left( r \right)$；（4）$r\left( {{E_{\max }}} \right)$；（5）$r\left( {{E_{\min }}} \right)$。$\Delta {E_{\max }}$和$\Delta {E_{\min }}$分别代表相对熵增量的最大值和最小值；$r\left( {{E_{\max }}} \right)$和$r\left( {{E_{\min }}} \right)$分别代表相对熵增量取最大值和最小值时的像素值；$Range\Delta E\left( r \right)$表示熵增量的幅度，可表示为：

${\text{Range}}\Delta E\left( r \right) = \Delta {E_{\max }} - \Delta {E_{\min }}$

(10)

${\text{Range}}\Delta E$越小，系统的混乱程度越低，分形特征越强。在本研究中，采用EXCEL编程计算灰度图像和二值图像孔隙结构的相对熵特征。

2 结果 2.1 添加生物质炭对CT孔隙度和孔径分布的影响

所有处理的CT孔隙度如表 2所示。2%生物质炭处理的CT孔隙度最大，而5%生物质炭处理的CT孔隙度最小。从土壤孔径分布来看，所有处理约有50%的孔隙半径在0~5 pixels之间，其中2017年5%生物质炭处理0~5 pixels孔隙占比最高，达到了67.76%。随着生物质炭施用年限的增加，2020年5%生物质炭处理下0~5 pixels孔隙占比降低了47%，10~25 pixels孔隙占比显著提高。此外，与0%和2%生物质炭处理相比，5%生物质炭最大孔隙半径更小（20~25 pixels）。

2.2 添加生物质炭对CT孔隙多重分形特征的影响

在进行多重分形分析前，需确定多重分形的计算尺度^[19]，即选择$\ln \left( {\chi \left( \varepsilon \right)} \right)$与$\ln \left( \varepsilon \right)$的直线部分进行计算（图 3）。由图 3可知，$\ln \left( {\chi \left( \varepsilon \right)} \right)$随着尺度增加先缓慢降低，然后突然下降，最后在大尺度上呈现线性关系。对于灰度图像，$\ln \left( {\chi \left( \varepsilon \right)} \right)$与$\ln \left( \varepsilon \right)$在8~512 pixels尺度存在线性关系；然而这种关系在二值图像中仅存在于64~512 pixels尺度中。因此，本文选择8~512 pixels尺度计算灰度图像的孔隙多重分形特征；选择64~512 pixels尺度计算二值图像的孔隙多重分形特征。

三维多重分形的结果如图 4所示。灰度图像和二值图像在各自的计算尺度下均显示出多重分形特征。在灰度图像和二值图像中，2%生物质炭处理的土样多重分形谱均更宽，表明2%生物质炭处理拥有更丰富更复杂的孔隙结构。此外，所有处理的多重分形谱均为短左钩形，这表明在所有处理中，大孔隙更为常见，而小孔隙占比较少，变异性指数更高。此外，灰度图像较二值图像具有更宽的多重分形谱，特别是0%生物质炭处理，其二值图像的多重分形谱几乎呈现一个点，这意味着二值图像孔隙的多重分形特征很弱。由此可知，灰度图像在更小的尺度上存在多重分形特征，且得到的孔隙多重分形特征更强。因此，本文采用灰度图像分析不同处理间孔隙的多重分形特征差异。

不同处理的多重分形参数如表 3所示。由表 3可知，不同处理的${\alpha _0}$、${\alpha _1}$和${\alpha _2}$差异不明显，表明这三个参数不能良好地区分生物质炭处理导致孔隙结构的变化。${{\Delta }}\alpha$能够量化孔隙结构的复杂度。2%生物质炭处理的${{\Delta }}\alpha$最大，较0%和5%生物质炭处理分别增加了14.7%和34.9%。此外，所有处理的${{\Delta }}\alpha$逐年增加，0%、2%和5%生物质炭处理的${{\Delta }}\alpha$年均增长率为16.7%、17.7%和19.2%，生物质炭处理的${{\Delta }}\alpha$年均增长率更高。所有处理在所有年份的${{\Delta }}f$均小于0，表明本试验所有土样的孔隙均以大孔隙为主，小孔隙的占比更低。

2.3 添加生物质炭对CT孔隙相对熵特征影响

灰度图像和二值图像相对熵结果如图 5所示。由图 5可知，灰度图像的相对熵和相对熵增量均接近于0，说明相对熵方法不适合灰度图像土壤孔隙结构的分析。在二值图像中，相对熵并未随着尺度的增加呈现良好的线性关系，而是会随着尺度的增加上下波动。通过相对熵增量的分析，发现相对熵增量随着尺度的减小先迅速增大，再缓慢下降，在16 pixels尺度下相对熵增量达到最大值0.767，这意味着在该尺度下，二值图像的自相似性发生破坏，图像不存在分形特征。

为了更直观化展示二值图像孔隙结构与尺度的关系，选择2017年0%生物质炭处理的一张二值图像开展相对熵分析，如图 6所示。由图 6（b）可知，相对熵增量在64至32 pixels时迅速增加，达到最大值0.81。图 6（c）和图 6（d）展示了在64 pixels和32 pixels网格下的二值切片图像，可以发现当计算网格由64 pixels降低至32 pixels时，部分计算网格几乎完全被孔隙填满，此时这些计算网格的自相似性发生破坏，从而导致图像的分形特征发生破坏。

二值图像三维相对熵参数如表 4所示。0%生物质炭处理的$\Delta {E_{\max }}$和${\text{Range}}({{\Delta }}E)$高于生物质炭处理；而2%生物质炭处理的$\Delta {E_{\max }}$和${\text{Range}}({{\Delta }}E)$较5%生物质炭处理分别降低7.97%和17.90%，说明2%生物质炭处理的土壤孔隙拥有最强的分形特征，而未添加生物质炭处理的分形特征最弱。此外，所有处理的$\Delta E$均在128 pixels处取得最小值；生物质炭添加后，$r\left( {\Delta {E_{\max }}} \right)$有降低趋势。年份同样影响相对熵结果，根据表 4，所有处理的${\text{Range}}({{\Delta }}E)$随着年份逐渐减小，其中2%和5%生物质炭处理的下降速率（8.28%和7.85%）均高于未添加生物质炭处理（7.29%），$r(\Delta {E_{\max }})$随着年份增加有增加的趋势。

3 讨论

本试验结合多重分形理论与相对熵理论探究生物质炭改良土壤微观孔隙的复杂行为，结果表明，灰度图像较二值图像在更小的尺度上存在更强的多重分形特征。这是因为当灰度图像转化为二值图像时，图像像素完全被孔隙或土壤基质填满，灰度图像中部分信息发生丢失^[20]，因此二值图像微观孔隙仅在更大尺度上存在多重分形特征（64~512 pixels），此时可用于计算孔隙多重分形特征的尺度仅有4个，使得到的孔隙多重分形特征很弱。而相对熵方法则通过放大二值图像的孔隙结构随着研究尺度变化的差异性，能够研究所有尺度下二值图像孔隙结构，更好地表征弱多重分形对象的多尺度异质性。综上所述，为分析土壤微观孔隙的复杂尺度行为，多重分形方法适合于灰度图像分析，相对熵方法适合于二值图像。

由灰度图像的多重分形结果（表 3）可知，2%生物质炭处理孔隙的复杂度最高，不对称性最强，这说明该处理中小孔隙占比更低，变异性更强。这可能是因为生物质炭的添加有助于促进水稳性大团聚体的形成，大团聚体能够显著增加孔隙的连通性^[21]，使得盐渍土中孤立小孔隙的数量显著降低，连通的大孔隙增加，从而使得土壤孔隙的复杂度与不对称性均提高。高生物质炭添加量（5%生物质炭）降低了盐渍土的孔隙复杂度与不对称性，这可能是因为生物质炭和土壤互相填充，导致大孔径的孔隙被占据而减少所致（图 2）^[22]。此外，不同处理之间的${\alpha _0}$、${\alpha _1}$和${\alpha _2}$差异不明显，表明这三个参数无法反映生物质炭处理对于土壤孔隙结构的影响，这与本团队之前的研究一致^[23]。相对熵的方法给出了类似的结果。由二值图像的相对熵参数可知（表 4），生物质炭处理的$\Delta {E_{\max }}$和${\text{Range}}({{\Delta }}E)$更低，说明生物质炭添加后土壤孔隙的多重分形特征增加。此外，随着生物质炭用量的增加，特别是5%生物质炭的处理，$r\left( {\Delta {E_{\max }}} \right)$有降低的趋势，这意味着土壤中孔隙的最大孔径有降低的趋势（图 2），这一结果证实了生物质炭会与土壤相互作用，破坏土壤中原有的大孔隙。

此外，生物质炭对于土壤孔隙结构有持续性影响。随着生物质炭施用年限的增加，所有处理的$\Delta \alpha$逐年增加而${\text{Range }}\Delta E$逐年降低，这意味着孔隙的复杂程度随着年份有增加的趋势。本研究在水稻种植前均会对测坑中的土壤开展翻耕，在此过程中改变土壤孔隙的界面形态特征使得土壤通透性增加，土壤中生物活动强度增大，生物孔隙与植物根系交错形成更为复杂的孔隙结构^[24]。生物质炭的添加则有助于增强土壤中生物和微生物的活动^[25]，因而孔隙结构复杂度拥有更快的变化速率。Soto-Gómez等^[26]指出更复杂的孔隙结构改善孔隙形态、增加孔隙数量和连通性，这意味着土壤导水透气性得到改善，生物质炭的添加能够提高滨海盐渍土孔隙复杂度并加快土壤孔隙结构的发育，从而实现土壤改良。

4 结论

应采用不同的理论方法分析土壤微观孔隙的复杂尺度行为，多重分形方法适合于强分形特征对象（例如灰度图像），相对熵方法适合于弱多重分形特征对象（例如二值图像）。2%生物质炭处理的孔隙结构复杂程度最高，土壤的孔隙结构改善效果最好。与0%和5%生物质炭处理相比，$\Delta \alpha$分别提高了14.7%和34.9%，${\text{Range }}\Delta E$分别降低了7.97%和17.90%。高生物质炭的添加量（5%生物质炭处理）减少孔隙结构的复杂度，降低生物质炭的改良效果。所有处理的孔隙复杂度随着年份逐年增加。生物质炭的添加能够加快土壤孔隙结构的发育，生物质炭处理的孔隙复杂度增加速率更快。2%和5%生物质炭处理的$\Delta \alpha$年均增长率为17.7%和19.2%，${\text{Range }}\Delta E$年均降低率为8.28%和7.85%，均高于0%生物质炭处理（16.7%和7.29%）。