岳海晶,樊贵盛

(太原理工大学 水利科学与工程学院,太原 030024)

考虑土壤结构变形的备耕地入渗参数线性预报模型

岳海晶,樊贵盛

(太原理工大学 水利科学与工程学院,太原 030024)

基于山西省境内大田土壤入渗试验和相关土壤基本理化参数测试数据样本,考虑头水地灌溉时地表土壤骨架变形、容重增加的客观实际,分析了结构变形对Kostiakov三参数入渗模型参数的影响,建立了考虑土壤结构变形的备耕地入渗参数线性预报模型结构；利用MATLAB程序对模型进行了模型系数的线性回归,最终得到了回归系数均显著,模型平均相对误差较小(均能控制在15%以下),预测精度较高的备耕地土壤水分入渗参数的线性预报模型,实现了利用容易获取的土壤基本理化参数进行头水地土壤水分入渗过程预测的有效方法。

土壤结构变形;Kostiakov三参数入渗模型;线性预报模型;土壤基本理化参数;误差分析

灌溉或降雨条件下,水分通过地表进入土壤并在土壤中储存和运动形成土壤水的过程称为土壤水分入渗,土壤水分入渗是田间水循环的一项重要环节,也是地表水转化成土壤水后供植物根系吸收的必要条件,对土壤水分入渗的研究能够为地下水资源的评价和合理指导农田灌溉提供科学的依据。国内外众多研究学者对土壤入渗机理进行大量研究分析后建立了不同类型的土壤水分入渗模型[1],实现了利用经验入渗模型参数来定量表征土壤水分入渗的过程。王维汉[2]总结了入渗模型参数估算的各种方法的适用性和优缺点;冯锦萍,樊贵盛[3]依据大田入渗试验数据研究了土壤入渗参数的线性预测模型;曹崇文[4]详细阐述了建立土壤入渗参数传输函数的三种手段。但是,上述对土壤入渗参数的研究均未考虑灌溉过程中土壤结构变形,因此研究成果有一定的局限性,尤其对备耕土壤局限性更大。备耕土壤是经过翻松后准备播种的土壤,特点是耕作层(0～20 cm)结构疏松,随着灌溉的进行,水分进入土壤导致表层土壤湿陷、崩塌,土壤结构发生了变化。研究表明,随着灌溉的进行,表层土壤的含水量急剧增加,表层土壤的容重会发生显著的变化,黄传琴[5]等通过研究干湿交替过程中土壤的胀缩特性分析了土壤体积变化与土壤含水率的关系;陈亮[6]等阐述了含水量的变化对土壤体积的影响机制;樊贵盛[7]等人明确了脱水过程中土壤密度与含水量存在二次函数的关系。上述研究表明,备耕土壤在灌溉过程中表层土壤颗粒骨架发生变形,导致土壤结构变化,变化趋势是随着含水量增加,表层土壤由疏松变密实,土壤密度由小变大。为使所建土壤入渗参数预报模型更符合实际灌溉过程,笔者以入渗参数物理意义明确的Kostiakov三参数入渗经验模型作为研究对象,利用山西省境内大田土壤入渗试验数据以及试验区土壤的基本理化参数,基于 MATLAB程序语言建立入渗参数的线性预报模型,得到更符合实际灌溉情况下考虑土壤结构变形条件的备耕土壤水分入渗参数的预报模型。

1 材料与方法

1.1 土壤条件

大田入渗试验涉及山西省境内由南到北八市十六县区,包括山区、丘陵、平原等多种地貌形态;试验区土壤类型丰富,有褐土、栗褐土、草甸土、黄绵土、红黏土以及棕壤土等多种类别;土壤质地类型包含有砂土、壤土、粘壤土、粘土等;土壤结构复杂多样,主要有团粒状、网粒状、块状、柱状、片状等多种形态[8]。为减小试验误差给试验结果带来的不利影响,在每个试验区选取5个试验点,土壤密度变化范围为0.774～1.589 g/cm3;体积含水率变化范围为1.75%～42.26%;土壤砂粒的质量分数的变化为19.70%～88.00%,粉粒变化为10.88%～76.10%,粘粒变化为0.92%～20.48%;土壤有机质质量分数为0.416%～6.59%。试验区土壤类型及上述土壤基本理化参数的变化范围覆盖山西省境内所有土壤,试验点分布范围广,数据代表性强。

1.2 试验方案

大田土壤入渗试验分别在灌溉前的翻松土壤和灌溉后的密实土壤进行,试验采用双套环入渗仪,将高度均为25.0 cm的外环(直径为64.4 cm)与内环(直径26.0 cm)预埋在深度为20.0 cm的土层,为保证试验的准确性,下环深度达到土壤犁底层。用自制的水位控制器将内外环积水入渗水头控制在2.0 cm,在规定的时间间隔用量筒为内环供水,记录量筒读数;外环水位与内环水位基本持平,无需记录外环加水量。试验研究表明,大田土壤水分入渗到90 min已基本达到稳定入渗状态,因此,以90 min作为入渗试验结束时间。

为表征土壤入渗参数与土壤基本理化参数的定量关系,除需获取土壤入渗试验数据外,还需获取土壤密度、体积含水率、土壤粒径机械组成及有机质质量分数等土壤基本理化参数数据。用100 cm3环刀切割未经扰动的自然状态土样,削除环刀体积外多余的土壤,得到一定体积的土样,将土样烘干称重得到的土样质量与土样体积相除得到土壤密度;将一定质量的土样烘干称重得到土壤质量含水率,结合测得的土壤密度求得土壤体积含水率;利用比重计法测定土壤砂粒、粉粒、粘粒含量得到土壤质地参数;采用重铬酸钾容量法测土壤中有机碳的质量分数,进而得到土壤中有机质质量分数。

1.3 样本数据建立

将大田入渗试验及土壤理化试验数据进行处理,过滤掉有明显错误的奇异值后建立两组数据样本,第一组为灌溉前的翻松土壤数据样本,第二组为灌溉后的密实土壤数据样本,每组样本均包含91组试验数据,试验数据包括:利用MATLAB软件将大田土壤入渗试验结果进行拟合得到的Kostiakov三参数入渗模型的模型参数K,α,f0;灌溉前翻松土壤及灌溉后密实土壤0～10，10～20，0～20，20～40 cm各土层的干容重;0～20，20～40 cm土层范围内的土壤粒径机械组成;0～10，10～20，0～20，20～40 cm各土层范围内的的体积含水率;0～20 cm范围内土壤有机质质量分数。所有试验均在非冻融期非盐碱地土壤上进行,将其中的3组样本数据列于表1。

2 备耕土壤水分入渗参数线性预报模型的建立

2.1 预测模型结构

表1 部分试验区土壤入渗试验样本数据

对未知参数的预测模型有多种形式,如线性预报模型、非线性预报模型、BP模型等,在满足精度要求的条件下,线性模型由于其结构形式简单、建模方法简易、模型应用方便等优点更适合于广大农民群众进行科学灌溉以及基层管理单位进行农业生产活动的管理,线性模型能够在很大程度上简化工作量,提高农业生产工作的效率。因此,本文选取适用范围较为广泛的线性模型进行Kostiakov三参数入渗模型三个参数的预测,预测模型结构采用如下的形式:

Y=β0+β1X1+β2X2+β3X3+…+βnXn.

(1)

式中:Y为预测变量,包括Kostiakov三参数入渗模型的K，α，f0;βi为输入变量的回归系数;Xi为第i个输入变量,此处表示常规土壤理化参数如土壤粒径质量百分比、土壤干容重、含水率和有机质质量分数等;n为变量个数。

2.2 预测模型输出变量

Kostiakov三参数入渗模型参数的物理意义明确,能够很好地模拟土壤水分的入渗过程, 文中将以三参数入渗模型作为研究对象:

i(t)=Kt-α+f0.

(2)

式中:i(t)为t时刻的入渗速率;K为入渗系数,表示入渗开始后第一个单位时间末扣除f0后的累积入渗量;α为入渗指数,表征土壤入渗速度的衰减速度;f0为土壤相对稳定入渗率,是在单位土壤势梯度下饱和土壤的入渗速度或非饱和土壤入渗达到相对稳定阶段的入渗速度。将三参数入渗模型的3个参数K，α，f0作为预测模型的输出变量。

2.3 预测模型输入变量

土壤密度可以衡量土壤结构,土壤结构表征在内外因素共同作用下形成的团聚体的形态、特征、数量及分布,土壤结构通过影响土壤的密实程度和孔隙的大小及分布来改变土壤水分运动的驱动力即土水势来影响田间土壤水分的入渗;土壤体积含水率作为衡量土壤水分含量的指标通过影响土壤层与层之间的吸力梯度来影响土壤水分的运动,对入渗能力产生较大的影响;土壤质地可用土壤中不同粒径土壤颗粒的百分比表示,砂粒、粉粒、粘粒的相对含量不同意味着土壤颗粒的表面能不同,孔隙大小、分布及连通性不同,影响了土壤水分运动的通道及水分子所受的吸力从而影响土壤水分的入渗;土壤有机质主要指土壤中含有无机碳的腐殖质等,无机碳含量高的物质能够增加土壤的粘结性有利于团聚体的形成,而团聚体的数量和分布又能够影响土壤良好孔隙的数量和分布从而增加土壤结构的稳定性,而且有机质丰富的土壤其毛管孔隙较多,能够将水分存储在土壤中,影响土壤水分的下渗[9]。

对于入渗系数K,表示的是第一个单位时间末的土壤入渗情况,历史较短,土壤水分最多能下渗到犁底层以上的土层范围,因此对于K而言表层土壤的理化参数对其影响较大,除此之外,在第一个单位时间内地表疏松土壤密度在灌水后还未来得及发生变化土壤入渗已完成,因此对K的预测要采用灌溉前的样本数据;对于入渗指数α,表征土壤入渗速度的衰减能力,土壤入渗途经土层(0～20 cm)范围内的土壤理化性质对其均有一定的影响;对于相对稳定入渗率f0,它主要表征入渗达到相对稳定时段的入渗率,0～20 cm范围内的土壤理化性质对其影响较明显,但是当土壤达到相对稳定阶段时地表土壤含水率已基本饱和,此时初始含水率对参数f0的影响已基本不存在,因此理化参数含水率对于f0而言已不能作为一项输入变量,当土壤入渗到相对稳定阶段时,地表土壤的容重变化也已基本稳定,地表容重对f0的影响很大。

土壤质地作为一项重要的土壤基本理化参数,其作为输入变量时通常以砂粒、粉粒、粘粒的相对质量分数来表征,三种土壤颗粒质量分数之和基本为1,为使模型结构的每一个输入变量独立,一般以两种土壤粒径作为输入变量,对不同的入渗参数,不同的土壤粒径组合对其产生的影响不同,因此需要通过数学手段定量分析每一个参数的具体输入变量。

对预测模型输入变量进行定性分析可知影响每个入渗参数的输入变量的数量、类型以及哪些土层范围的变量对其影响最明显,初步确定如表2所示。

2.4 预测模型结构

2.4.1 入渗参数K结构的确定

以砂粒和粉粒表征土壤质地的输入变量作为Ⅰ,以砂粒和粘粒表征土壤质地的输入变量作为Ⅱ,以粉粒和粘粒作为表征土壤质地的输入变量作为Ⅲ,对三种模型形式进行T检验,得到入渗系数K对的最佳预报模型结构,建立过程如下表3所示。

由检验结果可以看出,入渗指数K的三种模型形式均能满足数理回归分析的要求,回归系数T值均大于t0.025(2.00),如何选取最佳的组合形式需进行误差分析判断其平均相对误差选择精度最高的形式作为K的最终结构。

2.4.2 入渗指数α结构的确定

同样以表征质地的不同的输入参数对模型结构划分为A，B，C,T检验结果如表4所示。

表2 入渗参数影响因素

表3 入渗系数K结构建立T检验

表4 入渗系数α结构建立T检验

由上表可以看出,无论是砂粒和粘粒的组合B还是粉粒和粘粒的组合C模型回归系数的显著性均很差,且有机质含量的不显著与定性的理论分析相悖,只有砂粒和粉粒的组合A既能满足数理回归分析的要求又能满足定性理论分析的内容,因此入渗指数α的最佳模型为A。

2.4.3 相对稳定入渗率f0结构的确定

对f0的模型结构ⅰ、ⅱ、ⅲ进行同样的T检验,检验结果如表5所示。

表5 相对稳定入渗率f0结构建立T检验

由检验结果可以看出,砂粒和粘粒的组合ⅱ以及粉粒和粘粒的组合ⅲ均能满足数理回归分析的要求,回归系数T值均大于t0.025(2.00),说明ⅱ和ⅲ的变量组合均可作为相对稳定入渗率f0结构的输入变量,但是哪种模型形式最佳需要进行进一步的误差分析，判断其平均相对误差选择精度最高的形式作为f0的最终结构。

综上分析,确定了Kostiakov三参数入渗模型的预测模型结构,将模型结果列于表6.表中:γ0为0～10 cm土层的容重,g/cm3;γ1为0～20 cm土层的容重,g/cm3;θ0为0～10 cm土层的体积含水率,%;θ1为0～20 cm土层的体积含水率,%;w1为0～20 cm土层的砂粒质量分数,%;w2为0～20 cm土层的粉粒质量分数,%;w3为0～20 cm土层的粘粒质量分数,%;w为0～20 cm土层的有机质质量分数,%。

表6 三个模型参数预报模型结构

2.5 预测模型结构回归结果

2.5.1 回归方程系数

基于数理统计学最小二乘法原理编制MATLAB线性回归分析程序语言,将前面选定的91组样本数据代入已建模型结构进行回归分析,得到如表7所示回归结果。

表7 参数回归方程系数

2.5.2 模型显著性检验

预测模型线性回归结构建立以后需要采用方差分析法进行回归模型的显著性检验,将由91组样本数据进行计算得到的样本F值与给定的显著水平(α=0.05)对应的Fα(m,n-m-1)值比较,n为样本长度,m为输入变量个数,显著性结果列于表8。

表8 模型回归显著性检验结果

由上表可知,3个参数模型的F值(10.549 2～52.238 2)均大于F0.05(2.29～2.45),说明所建立的模型结构是显著的。

3 误差分析

3.1 参数预报模型平均相对误差分析

参数模型结构建立以后需对建立模型结构的样本数据的平均相对误差进行分析,判断其精度是否能够满足农业生产活动的需求。计算91组样本数据实测值和拟合值的绝对值误差,将全部误差进行平均得到参数预报模型的平均相对误差,列于表9。

表9 参数预报模型平均相对误差

由上述结果可知所建各参数预报模型的平均误差分别为11.66%～14.47%,均小于15%,说明预报模型精度较高,能够满足要求。对于入渗指数K,砂粒和粘粒的变量组合Ⅱ平均相对误差最小仅为11.66%,比组合Ⅰ低1.75%,选为K的最佳预报模型;对于相对稳定入渗率f0的两种模型变量组合形式,ⅱ的平均相对误差14.25%比ⅲ的平均相对误差14.47%低0.22%,说明ⅱ的精度相对较高,选择变量组合形式ⅱ作为f0的最终模型结构。

根据输入变量T检验建立回归系数均显著的模型结构,然后根据平均相对误差分析在所有显著性均显著的模型结构中选择精度最高的组合形式作为土壤结构变形条件下的备耕土壤水分入渗参数线性预报模型,结果如下:

K=3.822 4-1.388 1γ0-2.309 2θ0-

0.876 8w1+2.709 6w3+0.219 3w,

(3)

α=-0.111 3+0.140 8γ1-0.264 1θ1+

0.220 9w1+0.286 3w2+0.011 4w,

(4)

f0=0.115 1-0.028 9γ0+0.048 9w2-

0.234 0w3+0.019 0w.

(5)

3.2 预报参数下的Kostiakov三参数入渗模型的平均相对误差

为分析上述模型参数预报结构下的Kostiakov三参数入渗模型的预测精度,将每组样本数据的参数预测值代入三参数入渗模型,求得给定时间下的土壤入渗速率。实际灌溉或降雨条件下,60min时土壤入渗已基本达到稳定阶段,但一般情况下为了安全起见,常以90min时的土壤入渗速率来衡量土壤水分入渗能力[10]。由参数预测值计算得出的91组样本数据的90min土壤入渗速率与实测90min土壤入渗速率进行比较,得到90min土壤入渗速率的平均相对误差为11.79%,精度能够满足农业生产活动的需要。

3.3 备耕土壤水分入渗参数预报模型实例验证

选取大同新平镇、朔州边耀、晋城郎庄村、运城南方平、长治郜村进行所建模型结构的实例预测。方法如下:对5个地方的各层土壤进行基本理化试验,得到灌溉前后的土壤密度、体积含水率、砂粒、粉粒、粘粒质量分数以及有机质质量分数等基本理化参数,代入入渗系数K、入渗指数α以及相对稳定入渗率f0的预测模型结构,得到相应的参数预测值即可得到试验地区的Kostiakov三参数入渗模型,由下表10结果可看出参数相对误差范围为2.63%～16.74%,在允许误差范围内,预测结果可对农业生产活动进行指导。

表10 参数预报模型应用实例结果

4 结论与建议

考虑灌溉时土壤骨架变形的情况,对Kostiakov三参数入渗模型参数K，α，f0的输入变量进行定性的理论分析以及定量的数理统计T检验，建立了预测参数与土壤基本理化性质土壤干容重、体积含水率、土壤质地及有机质间的线性结构,采用MATLAB程序语言对建立的模型结构进行线性回归分析,得到了回归系数均显著、模型精度较高的备耕土壤水分入渗参数的线性预报模型。入渗系数K的平均相对误差为11.66%,入渗指数α的平均相对误差为13.09%,相对稳定入渗率f0的平均相对误差为14.25%,90min土壤入渗率的平均相对误差为11.79%,均能将误差控制在15%以下,并对预测模型结构进行实例预测,验证其可行性,结果表明预测结果能够满足精度要求,可用来指导农业生产活动。

所建土壤结构变形条件下的备耕土壤水分入渗参数的线性预报模型精度能够满足需要,但是仅考虑了影响土壤入渗模型参数的主要影响因素,没有考虑如地表结皮、水温、地温等次要影响因素,在以后的研究过程中可以适当将次要影响因素考虑在内,进一步提高预测精度。

[1] 寇小华,王文,郑国权.土壤水分入渗模型的研究方法综述[J].亚热带水土保持,2013(3)：53-55.

[2] 王维汉,陈晓东,等.土壤入渗参数的估算方法及其变异性研究进展[J].中国农学报,2011,27(6):272-275.

[3] 冯锦萍,樊贵盛.土壤入渗参数的线性传输函数研究[J].中国农村水利水电,2014(9):8-11.

[4] 曹崇文.利用土壤传输函数确定入渗参数的方法研究[D].太原：太原理工大学,2007.

[5] 黄传琴,邵明安.干湿交替过程中土壤胀缩特征的实验研究[J].土壤通报,2008,39(6):1243-1247.

[6] 陈亮,卢亮.土体干湿循环过程中的体积变形特性研究[J].地下空间与工程学报,2013,9(2):229-235.

[7] 樊贵盛,孔令超,等.土壤脱水过程中结构与含水率间的定量关系[J].灌溉排水学报,2012,31(5):40-43.

[8] 韩永鸿，樊贵盛.土壤结构与田间持水率间的定量关系研究[J].太原理工大学学报,2012，43(5):615-619.

[9] 刘耀宗,张经元.山西土壤[M].北京：科学出版社,1992:45-294.

[10] 李卓,刘永红,等.土壤水分入渗影响机制研究综述[J].灌溉排水学报,2011,30(5):124-130.

[11] 樊贵盛,李雪转,李红星.非饱和土壤介质水分入渗问题的试验研究[M].北京:中国水利水电出版社,2012:63-81.

[12] 雷志栋,杨诗秀,谢森传.土壤水动力学[M].北京:清华大学出版社,1988.

[13]WostenJHM,PachepskyYA,RawlsWJ.Pedotransferfunctions:bridgingthegapbetweenavailablebasicsoildataandmissingsoilhydrauliccharacteristic[J].JournalofHydrology,2001,251:123-150.

(编辑：贾丽红)

The Linear Prediction Model of Soil Infiltration Parameters for Preparation Land Considering of Deformation of Soil Structure

YUE Haijing,FAN Guisheng

(CollegeofHydroscienceandEngineering,TaiyuanUniversityTechnology,Taiyuan030024,China)

On the basis of the field infiltration and soil physicochemical parameters test data of Shanxi province,the influence of structural deformation on the infiltration parameters of Kostiakov model in three parameters was analyzed by taking into account the actual irrigation situation of surface soil skeleton deformation and soil volume-weight increase when watering the preparation land at first irrigation,and the linear prediction model structure of infiltration parameter was established considering the deformation of soil structure. Finally, the linear prediction model of soil infiltration parameters for preparation land was selected by using MATLAB program to conduct the linear regression of model coefficient.The regression coefficients were significant, the average relative error of the model was relatively small (under 15%) and the prediction accuracy was high.This effective method realized the prediction of soil water infiltration process of preparation land at first irrigation with readily available soil physicochemical parameters.

deformation of soil structure;Kostiakov infiltration model in three parameters;linear prediction model;soil physicochemical parameters;error analysis

1007-9432(2015)05-0616-07

2014-05-20

国家自然科学基金项目：区域尺度上土壤入渗参数多元非线性传输函数研究(40671081);山西省科技攻关项目(20100311124)

岳海晶(1989-),女,山西清徐人，硕士生,主要从事灌排理论与技术的研究，(E-mail)yuehaijing8811@163.com

樊贵盛(1955-)，博士生导师,教授,主要从事土壤物理、灌排理论与技术的研究，(E-mail)fanguis5507@263.net

S152;TV93

A

10.16355/j.cnki.issn1007-9432tyut.2015.05.026