首页 > 分享 > 基于分形理论研究土壤结构及其水分特征关系

基于分形理论研究土壤结构及其水分特征关系

花匠小妙招
2024-11-21 09:59

0. 引言

土壤是一种由不规则颗粒组成、具有一定意义下自仿相似结构的分形集(Pachepsky et al., 2000)，分形维数(简称分维)能够定量刻画分形集的几何和统计复杂性，分维的微小变化可以引起形状的急剧改变(Mandelbrot，1983；Chakraborti et al., 2003).

土壤粒度分布(particle size distribution，简称PSD)的分维能够表征土壤质地及其结构特征(Ersahin et al., 2006)，被认为是土壤最基本的物理性质之一，控制着土壤的水力性质.土壤水分特征曲线(soil water characteristic curve，简称SWCC)是表征非饱和土壤含水量或饱和度θ与基质吸力S的定量关系曲线，反映了土壤持水性能及水分运动特征(Nam et al., 2010)，受土壤质地、结构、容重(bulk density)、孔隙率等因素影响(Phoon et al., 2010；Fu et al., 2011；Rajesh et al., 2017).

目前尚不能根据土壤基本性质推导得出土壤水分特征曲线，只能用试验方法测得，然而准确测定这一关系非常困难，已有的研究表明，土壤结构主要影响0~100 kPa范围内的土壤水分特征曲线形状，而曲线其余部分的形状受土壤质地、有机质含量以及粘土矿物影响较大(Brady and Weil, 2008).

分形是描述多孔介质物理特征的有效方法，研究模型包括质量分形、表面边界分形(张超等，2019)、孔隙分形(杜书恒等，2019)、孔隙‒固体分形、毛管束分形和混合分形单元等.根据孔隙空间的迭代次数是否无限，分形理论研究水分特征曲线可分为分形数学模型和准分形物理模型(蔡建超和胡祥云，2015).Wei et al. (2019)基于van Genuchten模型对不同质地石油污染土壤系统的保水性进行分析，认为持水性在某种程度上取决于土壤结构，尤其在低吸力阶段这种作用更明显.Jin et al.(2019)基于分形理论，考虑了薄膜流，提出了一种能够连续描述整个土壤水含量范围内的水分特征曲线模型，模型与实验数据吻合良好.

一方面，受限于实验设备仪器的测量精度与范围，加之土壤结构易受扰动，土壤水力特征参数实验测定复杂，耗时较长，前人研究大多集中于理论模型推导(Sillers et al., 2001)，实验数据验证相对较少，其适用范围和预测精度有待检验.另一方面，目前对于土壤结构特征和土壤水力性质的研究较多，而基于分形理论探讨土壤粒度分布与土壤水力特征之间关系的研究却相对较少，Bird et al.(2000)提出用粒度分布预测土壤水分特征曲线的方法，该方法以质量分布估计样品的孔隙‒固体界面分维，需要先估测孔隙和固相占比.黄冠华和詹卫华(2002)应用Menger海绵结构推导出了土壤水分特征曲线的分形模型，该方法同样基于土壤颗粒质量分形维数，但土壤孔隙结构复杂且不能由一个确定的Menger模型产生.激光散射技术的广泛应用提高了PSD测试速度，为本文探索土壤颗粒体积分形维数与水力性质之间的关系提供了基础.

为定量表征土壤水分运动，本研究选取华北平原子牙河流域为研究对象，通过优选点位野外探坑，采用直径为15 cm、高为10 cm的环刀获取原位原状土样，在实验室用张力计法直接测量土壤样品水分特征曲线，用激光粒度分析仪测定土样粒度分布，基于分形理论，借助统计学原理，将土壤结构以分形形式综合量化表达，采用试验测定与模型验证相结合的方法对水分特征曲线进行分析，寻求华北平原山前倾斜平原‒冲积平原‒滨海平原地区土壤结构与水分特征曲线的联系，研究土壤结构定量指示水力性质.该研究可为量化区域土壤结构、解释土壤水分运动特征提供科学依据和技术支撑.

1. 研究区概况

子牙河位于112.00°~117.35°E，36.150°~39.30°N，属于海河水系，西起太行山东麓，东临渤海，南邻南运河，北至大清河，子牙河在天津与西河闸、南运河和大清河相交汇入渤海，流域总面积46 868 km2，河北省内流域面积占比72.6%(Ding et al., 2016)，流域内土壤母质主要由第四纪沉积物所组成.本研究在华北平原子牙河流域6个不同地点、不同深度采集了14个土壤样本，采样点分布在河北省石家庄市、衡水市、沧州市和廊坊市等省辖行政区域，以及天津市滨海新区(图 1)，覆盖了山前倾斜平原、冲积平原和滨海平原地貌，采样位置高程在2~73 m，取样深度在0.43~1.50 m之间.

图 1 采样点分布

Fig. 1. Schematic diagram of sampling sites distribution

正定(ZD)采样点位于太行山中段山前倾斜平原地带，其地貌类型属于滹沱河冲洪积扇的中上部，土层深厚，土壤质地较均匀，0~5 m主要为粉砂质粘壤土、砂壤土和砂土(何雨江等，2013；林丹等，2014).冀州市(JZ)、河间市(HJ)和大城县(DC)采样点位于冲积、泛滥平原地带，地势平坦，地表可见明显的龟裂现象，土色浅淡，质地以壤质为主.滨海新区采样点(一疙瘩村YGD和营房村YF)位于滨海平原地带，地区海拔低于5 m，地表组成物质主要为黏土，受海潮和河流双重影响，地势低平，土壤质地粘重，排水、通气不良(郝翠等，2009；杨晓潇等，2019).采样点地理位置见表 1.

表 1 采样点地理位置

Table Supplementary Table Geographical data of soil sampling locations

地理位置样品编号坐标高程(m) 取样深度(m) 地貌 ZD S1 114°34'36.10"E, 38°08'14.10"N 73 0.43~0.53 山前平原 S2 0.9~1.0 S3 1.40~1.50 JZ S4 115°17'55.42"E, 37°34'58.09"N 27 0.87~0.97 冲积平原 S5 1.1~1.2 HJ S6 116°07'55.00"E, 38°23'53.00"N 13 0.67~0.77 冲积平原 S7 1.25~1.35 DC S8 116°37'37.30"E, 38°39'32.88"N 7 0.7~0.8 冲积平原 S9 0.87~0.97 YGD S10 117°32'24.00"E, 39°00'48.00"N 2 0.5~0.6 滨海平原 S11 0.7~0.8 YF S12 117°34'12.70"E, 38°59'25.43"N 2 0.35~0.45 滨海平原 S13 0.6~0.7 S14 0.96~1.06

2. 材料与方法

2.1 土壤水分特征曲线测定

野外环刀(容积1 767 cm3)取样，室内用张力计法(量程85 kPa)测定样品水分特征曲线.首先将采集到的土样放置在水中浸水饱和，直至土样的上部能看到水或能触摸到水为止，一般重黏土需10 d，砂土需2~10 h.然后，校正张力计，调平天平，连接CR-1000数据监测采集系统(Campbell Scientific公司，美国)、SKYE Mini张力计(英国)和R2000系列电子计重计数天平(渠道科技有限公司，中国北京)组成测定系统，建立通讯，设定采集数据的时间间隔，室内温度保持在25 ℃，开始实验，张力计中的自由水经过陶土头与土壤水建立水力联系，当仪器内外的水势值趋于平衡时，仪器中水的总水势等于土壤中土水势，通过张力计测量土壤水吸力.采集土柱样品在自然蒸发状态下的天平读数和对应的张力计读数.最后，待张力计读数至80 kPa附近时停止实验，取出张力计，拆除测量系统，将样品置于DHG-9146A型电热恒温鼓风干燥箱(精宏实验室设备有限公司，中国上海)烘干.设置温度为105 ℃，砂壤土烘干6~8 h，黏土10 h，烘干至恒重称重(马传明，2013).

计算各个采集时刻对应的土壤体积含水率，结合对应时刻的张力计读数，从而得到该土壤样品的水分特征曲线.

2.2 土壤PSD测定

将测过水分特征曲线的环刀土柱对应层位上的样品烘干，混合均匀后过2 cm筛，用4分法取出粒度分析样，使用QT-2002型激光粒度仪(渠道科技有限公司，北京，中国)测试每个样本粒径累积分布，测量范围为0.1~600.0 μm，准确性误差 < 1%，重复性偏差 < 1%.

3. 结果

3.1 土壤粒度分布

粒度分布分析在判定水动力条件和粒径变化趋势等方面具有显著意义(牛宏等，2016).不同采样点垂向上选取土壤质地变化处环刀取样，研究区内从上游山前平原ZD采样点到下游滨海平原的YF采样点，测定的粒度分布结果见图 2.累积粒度分布曲线反映了颗粒级配范围、颗粒分选及水动力条件特征，显著上升段越陡(大部分在[10 μm，50 μm]的区间内)，表明土壤颗粒分布越均匀.在平面上，ZD采样点最为宽缓，颗粒级配范围广，流域内从上游山前平原到下游滨海平原，水动力条件变弱，土壤颗粒累积粒度分布变陡，粒径分布变窄，粒径变细，分选性增强，累积分布散点图左移.粒径变化趋势在纵向上无明显规律.

图 2 粒径分布

曲线图为粒径频数分布；散点图为粒径累积分布

Fig. 2. The distribution of grain size constitutes

3.2 分维模型与计算结果

本研究进行区域尺度不同深度土壤粒径的分维计算，包括累积分布分维和局部分段分维，其模型可表示为(Tyler and Wheatcraft, 1992；Peng et al., 2015)：

[didmax](3−D)=V(δ＜di)V0," role="presentation">[didmax](3−D)=V(δ＜di)V0,

(1)

式中，di(i=1、2、3、…、n)为土壤颗粒直径(μm)，dmax为土粒最大直径(μm)，D为分形维数，δ为岩土颗粒的粒径变量(μm)，V(δ < di)为粒径小于di的土粒体积(μm3)，V0为土壤样品体积(μm3).该公式适用于已知粒度体积分布求取分维.

对公式(1)两边分别取对数，得到：

(3−D)lg⁡didmax=lg⁡(δ＜di)V0." role="presentation">(3−D)lg⁡didmax=lg⁡(δ＜di)V0.

(2)

计算分维时，首先由实测粒径及其体积数据计算出lg⁡V(δ＜di)V0" role="presentation">lg⁡V(δ＜di)V0，而后以lg⁡didmax" role="presentation">lg⁡didmax为x坐标，lg⁡V(δ＜di)V0" role="presentation">lg⁡V(δ＜di)V0为y坐标，基于最小二乘法进行线性拟合，获得其斜率k，则D=3-k，如图 3a所示，数据点分布很难拟合成一条符合良好的直线，而是呈现出明显的分段特点，表征整体拟合优度的R2只有0.732 6.为更加准确描述土壤粒度分布局部特征，对数据进行多次逐段分析，分别选取激光粒度仪测定区间内相关系数最高的两个区间[5 μm，30 μm]和[10 μm，50 μm]作为特征区间，对于最大粒径小于50 μm的粒度分布，取最大值作为区间右边界，如滨海新区的YGD和YF采样点，按照上述方法拟合，即分段分维数Dn是特征区间In的函数，分别计算该区间内的分段分维，如图 3b、3c所示，两个区间内的分维值分别为2.000 3和2.365 2，明显高于整体分维值1.539 5，两段拟合优度均好于整体，分段分维数较为准确地表征了土壤粒度分布局部特征.

图 3 正定采样点S1分维拟合图

x轴为

lg⁡didmax" role="presentation">lg⁡didmax

；y轴为

lg⁡V(δ＜di)V0" role="presentation">lg⁡V(δ＜di)V0

Fig. 3. Fractal dimension fitting diagram of ZD sampling (S1)

基于上述模型，各个采样点分维计算结果见表 2.表 2显示土壤累积分布分段分维规律明显，在采样区域和深度范围内，从流域上游到下游，尤其是中下游段(HJ~YF采样点)，随着粒度变细，分段分维D1和D2均呈现出变大趋势，JZ和YGD采样点随着深度增加，黏粒含量变高，分维变大.HJ和DC采样点随着深度增加分维变小，ZD采样点D2随着深度增加分维先变小再变大，这与图 2土壤粒度累积分布一致，说明粒度在[10 μm，50 μm]区间内的分段分维值可以较好地表征土壤粒度累积分布上升段的变化趋势.

表 2 采样点分段分维数

Table Supplementary Table Subsection fractal dimension of sampling point

采样点编号粒径特征区间(μm) 分段分维数 I1段 I2段 D1 D2 ZD S1 5~30 10~50 2.000 3 2.365 2 S2 5~30 10~50 2.039 6 2.194 3 S3 5~30 10~50 2.365 8 2.679 1 JZ S4 5~30 10~50 0.963 3 1.445 2 S5 5~30 10~50 2.262 9 2.739 3 HJ S6 5~30 10~50 1.780 7 2.164 0 S7 5~30 10~50 1.238 4 1.639 6 DC S8 5~30 10~50 2.380 9 2.790 6 S9 5~30 10~50 1.970 0 2.402 4 YGD S10 5~30 10.000~42.603 2.614 2 2.872 8 S11 5~30 10.00~38.93 2.740 3 2.918 0 YF S12 5~30 10.000~35.573 2.757 2 2.922 9 S13 5~30 10.000~35.573 2.771 5 2.926 3 S14 5~30 10.00~38.93 2.783 3 2.925 7 注：对于最大粒径小于50 μm的粒度分布，取粒径最大值作为区间右边界，如滨海新区的S10~S14采样点.3.3 土壤水分特征曲线特征

本文采用试验测定与模型验证相结合的方法对水分特征曲线进行分析，用张力计法测得的实验结果拟合Gardner幂函数(Gardner et al., 1970)经验公式，结果如图 4所示，经验公式为：

式中，S为土壤基质吸力(kPa)，θ为体积含水量(%)，a和b为拟合参数.

由图 4可知，流域内从上游山前平原到下游滨海平原，SWCC由陡峭变为平缓，对照图 2土壤颗粒累积粒度分布曲线，说明黏粒含量越高，含水量随吸力变化越趋缓，在同一吸力条件下土壤的含水率越大.与山前平原和冲积平原相比，滨海平原在同一土壤体积含水率条件下，其吸力值最大(例如在含水率为38%时，山前平原平均吸力8.5 kPa，冲积平原区为24.6 kPa，滨海平原区为62.5 kPa)，这是因为土壤中黏粒含量增加使得细小孔隙发育，黏粒分选好，孔径分布较为均匀，在张力计的量程范围内，随着吸力的变大含水率缓慢减小，因此曲线显得格外平缓.在山前平原和冲积平原区，粒径较粗，孔隙较大，毛管力较弱，大孔隙中的水容易排出，土壤中仅有少量水存留，SWCC呈现出一定吸力以下陡峭(例如山前平原20 kPa以下)，而吸力较大时平缓的特点，曲线形态“上陡下缓”，这是因为中小孔隙毛管力较强，水分较难排出.

图 4 流域内不同位置采样剖面土壤水分特征曲线

Fig. 4. SWCC curves at different sampling positions and soil depths

Gardner幂函数经验公式中，参数a、b的拟合结果见表 3，系数(R2)表征模拟值与实测值的拟合效果，值越大拟合程度越好.幂函数经验公式对该流域土壤水分特征曲线拟合优度好，除JZ采样点外，其余采样点R2均大于0.9.对照图 4和表 3可知，公式(3)中的参数a决定了曲线的斜率，即土壤含水量随土壤吸力增加而减少的速度，a值越大，曲线越平直，持水能力越强，变化越慢，如冲积平原区(S6~S9).b值决定曲线的位置，b值越大，曲线越远离x轴，如山前平原区和滨海平原区(S1~S3、S10~S14).参数a和b的大小主要受控于土壤结构(不同粒径级配和黏粒含量).

表 3 不同采样点下幂函数经验公式拟合结果

Table Supplementary Table The fitting results of empirical formula in different sampling points

采样点编号幂函数参数拟合结果 a b R2 ZD S1 56.455 3 -0.162 38 0.916 6 S2 70.747 4 -0.315 4 0.922 7 S3 55.502 7 -0.330 1 0.950 0 JZ S4 180.803 6 -0.529 1 0.767 5 S5 77.985 9 -0.213 4 0.850 5 HJ S6 114.335 7 -0.412 3 0.979 2 S7 213.844 1 -0.534 9 0.987 7 DC S8 44.891 58 -0.147 88 0.973 3 S9 52.983 8 -0.150 55 0.949 8 YGD S10 42.066 2 -0.047 56 0.945 7 S11 47.767 4 -0.038 6 0.937 3 YF S12 39.325 6 -0.116 05 0.966 9 S13 40.956 4 -0.048 6 0.981 8 S14 49.236 3 -0.060 2 0.947 8

4. 讨论

4.1 土壤PSD分维与SWCC拟合参数的相关关系

前人的研究结果表明，土壤分形结构与其水力性质关系密切(程冬兵等，2009；王艳艳和何雨江，2019).本文利用SPSS 26软件分析公式(2)拟合参数与分维Pearson相关性，衡量其是否存在线性关系，置信水平P＝0.05，即P < 0.05为显著相关，P < 0.01为极显著相关，由于两两变量之间的关联程度未知(正相关、负相关或不相关)，因此本研究应用双尾检验(该检验一般被用于没有强烈方向性期望的实验中)，研究土壤累积PSD分维、分段分维与其Gardner水分特征曲线模型参数的相关关系，判断变量之间的相关性，结果见表 4.

表 4 土壤结构分维与土-水特征曲线拟合参数的Pearson相关系数

Table Supplementary Table Pearson correlation coefficient between FDs and SWCC fitting parameters

项目 D D1 D2 a -0.494 -0.89** -0.901** b 0.16 0.901** 0.906** 注：**表明在p < 0.01水平上极显著相关.

由表 4可知，土壤水分特征曲线的拟合参数(a、b)与土壤累积PSD分维不具有显著相关关系，然而拟合参数与特征区间的分段分维具有极显著相关关系，尤其是与粒度分布在[10 μm，50 μm]区间内的分段分维值相关关系更强，a与分段分维呈极显著负相关，分维值越小a越大；b与分段分维呈极显著正相关，分维值越大b越大，该结果为应用回归分析寻找粒度分布分维与水分特征曲线参数之间的定量关系提供了基础.这说明土壤粒度在[10 μm，50 μm]区间内的颗粒分布强烈地控制着0~80 kPa范围内的土壤水分特征曲线形状，即前文所示累积粒度分布图的显著上升段(图 2)主要控制着土壤水分特征曲线特征(图 4).

4.2 土壤PSD分维表征土-水特征曲线拟合参数

土壤颗粒分布和孔隙结构均具有分形特征，常用的分形结构有Sierpinski垫片、Sierpinski地毯、Menger海绵以及Koch曲线等(Tyler and Wheatcraft, 1990；Baliarda et al., 2000；Bird et al., 2000)，已有的SWCC分形模型需要已知饱和含水率(θs)、残余含水率(θr)或进气值等参数(de Gennes，1985；Rieu and Sposito, 1991；Tyler and Wheatcraft, 1992；Bird et al., 2000；Xu and Dong, 2004)，对于张力计来说无法测得黏土脱湿后的残余含水率(远大于张力计量程).

为此，本文结合4.1讨论结果，选取土壤粒度值[10 μm，50 μm]作为特征区间，以该区间内分段分维作为特征值，建立矩阵拟合各采样点的土壤分段分维与土壤水分特征曲线参数的关系，以D2为自变量，a、b为函数，置信度设为95%，a、b与D2的关系式如图 5所示，参考Jin et al.(2019)文中的公式：

ln⁡θ=(3−D)(−ln⁡S)+(3−D)ln⁡Sa+ln⁡θs," role="presentation">ln⁡θ=(3−D)(−ln⁡S)+(3−D)ln⁡Sa+ln⁡θs,

(4)

图 5 a和b分别与D2的拟合关系

Fig. 5. a & b are fitted with D2, respectively

式中，Sa为进气值(kPa)，θs为饱和体积含水率(%).

本文公式(3)中，a=k1(3-D)，b=k2(D-3)，通过线性回归分析得出k1≈-100.78，k2≈1/3，获取适用于子牙河流域的土壤PSD和SWCC的分形关系式(5)：

θ=100.78×(3−D)S(D−3)/3," role="presentation">θ=100.78×(3−D)S(D−3)/3,

(5)

方程两边取对数，可得公式(6)：

lg⁡θ=2.0+lg⁡(3−D)+D−33lg⁡S," role="presentation">lg⁡θ=2.0+lg⁡(3−D)+D−33lg⁡S,

(6)

式中，θ为体积含水量(%)，D为分形维数，S为土壤基质吸力(kPa).

本研究的模型所需参数少，测试简单，利用土壤粒度分布分形特征即可指示水力性质.

5. 结论

(1) 土壤粒度分布具有分形特点，从子牙河流域上游到下游，随着粒度变细，分段分维D1和D2均呈现出变大趋势，尤其是分维值D2能够表征土壤粒度累积分布上升段的变化趋势.

(2) 张力计法试验测定的土壤水分特征曲线结果表明，流域内从上游山前平原到下游滨海平原，SWCC由陡峭变为平缓，对照土壤颗粒累积分布曲线，说明黏粒含量越高，含水量随吸力变化趋缓，在同一吸力条件下土壤持水能力越强.

(3) 本文用幂函数模型拟合试验测定结果，所需参数个数少，对该流域土壤水分特征曲线拟合优度好，除JZ采样点外，其余采样点R2均大于0.9.模型中参数a决定了曲线的斜率，即土壤含水量随土壤吸力增加而减少的速度；b值决定曲线的位置，b值越大，曲线越远离x轴.参数a和b的大小，主要受控于土壤结构(不同粒级土壤排列组合及黏粒含量).

(4) 土壤PSD分维与土壤水分特征曲线拟合参数的相关关系分析结果表明，子牙河流域土壤粒度在[10 μm，50 μm]区间内的分布分维强烈地控制着0~80 kPa范围内的土壤水分特征曲线形状，以该区间的分段分维为纽带，获得0~80 kPa吸力范围内土壤水分特征曲线的关系式：.

致谢:本次研究野外取样工作到了华北水利水电大学王艳艳和河南理工大学徐流洋的帮助，成文过程中得到中国地质科学院水文地质环境地质研究所陆川研究员的指导，审稿专家和责任编辑老师提出了宝贵的修改意见，在此一并表示诚挚的谢意！