盐碱土壤Philip入渗模型参数的非线性预报模型

沈 婧，樊贵盛

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

0 引 言

据粮农组织不完全统计，全世界盐碱地面积为9.544 亿hm2，其中我国约为0.991 亿hm2，是世界第三大盐碱地国家。我国盐碱地广泛分布在华北、东北、西北及滨海地区，其中盐碱化耕地和盐碱荒地总面积达0.333 亿hm2，约占耕地总面积的27.5%，作为一种重要的可开发利用的土地资源，采取合理措施改良盐碱土以提高土地利用率和生产率对缓解我国紧张的人地矛盾意义重大。为此，国内众多学者在盐碱地改良方面进行了诸多探索[1-4]，包括灌水洗盐、平整深耕、选种育苗、间作方式、施用改良剂、田间灌排等，取得了显著成果。

盐碱地的形成受到气候条件、地理条件、土壤条件和地下水等多种因素的综合影响，但其形成的实质是各种可溶性盐分的随水运移，地表集聚。目前，国内外盐碱地的改良多采用灌排工程技术，其原理是利用灌溉水入渗实现对土壤盐分的有效淋洗。盐碱地高含盐量尤其是高钠离子含量使土体结构性差，水力传导度低，水分入渗缓慢，一方面造成了灌溉水利用率低，洗盐效果不佳，同时也加剧了干旱半干旱地区的水资源短缺。因此，合理确定盐碱地水分入渗参数对提高灌水利用率，实现节水型盐碱地改良意义重大。吕殿青等[5]在室内进行了一维积水入渗试验，分析了土壤盐分特征量与水分入渗量之间的关系，最终得出两者的函数关系。张莉等[6]进行了夹砂层土壤入渗过程水盐运动规律的室内模拟试验，观测了在灌水条件下夹砂层和均质土壤的水分入渗过程并对比分析了夹砂层土壤的盐分运动规律，该研究为滨海盐碱地改良和开发利用提供了理论依据。PATEL等[7]通过对印度古吉拉特邦北部的盐碱土壤进行的水分入渗试验发现，入渗率受钠离子吸附、土壤黏粒含量的影响，建议要依据盐碱地的类型选择合理的种植作物。任长江等[8]通过建立水分累计入渗量、入渗率与初始含盐量的多元回归方程，分析了土壤初始含盐量对入渗特性的影响，结果发现：相同时间初始含盐量越高，累计入渗量和入渗率越小。纵观现有研究成果，在盐碱地水盐运移规律及土壤水分入渗机理方面研究较多，成果显著，但试验法获取盐碱地土壤水分入渗参数工作量大，耗时费力，因此，本文以Philip[9]半理论半经验入渗模型参数为研究对象，采用土壤传输函数的方法，对晋北盐碱地区的土壤水分入渗参数进行非线性预测，为间接获取盐碱地土壤水分入渗参数和盐碱地改良提供了一种新思路。

1 材料与方法

1.1 试验区概况及土壤条件

土壤水分原位入渗试验在山西省北部朔州市应县盐碱地进行，地势东南高，西北低，属北温带大陆性季风气候，四季分明，年际温差大，据应县气象局多年观测资料，日温差可达17.1 ℃，年温差高达74.1 ℃，年均气温7 ℃左右，降雨分布极不均匀，主要集中在6-9月，年均降水量为380 mm，年均蒸发量为1 918.6 mm(20 cm蒸发皿)，年均蒸发量远大于年均降雨量。

该县属海河流域和桑干河水系，洪积倾斜平原地貌特征明显，地下水资源丰富，地下水位最高为1.5 m，矿化度落在0.5～1.5 g/L之间。试验区土壤碱化程度相对较轻而盐渍化程度相对较重，pH值介于7.0～8.3之间，电导率范围为0.123 7～1.916 6 S/m2(0～20 cm)和0.125 0～3.158 1 S/m2(20～40 cm)，试验区土壤理化参数情况见表1。

1.2 试验设备与方法

盐碱地土壤原位入渗试验采用双套环入渗仪法，双套环为2个同心圆环(内、外环直径分别为26 、60 cm，两环高度均为25 cm)，内环控制试验的入渗面积，外环起保证内环水体垂直下渗、减少侧向渗流的作用。试验开始前，将双套环打入地表以下20 cm的深度。用自制的水位控制器对内外环水位进行控制，且水位基本保持在地表以上2 cm处。试验过程中，准确记录内环加水的时间和注水量，以便日后拟合入渗参数。试验研究表明，土壤水分入渗在60 min左右基本达到稳定入渗，本试验选定入渗时间为90 min，以确保大田土壤水分入渗都达到稳定入渗状态。

在进行野外试验的同时，还需在室内对土壤理化指标进行测定，包括土壤全盐量、土壤容重、体积含水率、土壤有机质、土壤质地。土壤全盐量采用电导法测定；土壤容重通过蜡封法进行测定；体积含水率测定用传统烘干称重法获得；试验点土壤有机质量的测定是利用化学方法通过重铬酸钾滴定法得到；土壤质地通过比重计法得到筛分曲线，然后分析土壤的颗粒级配，进而确定土壤质地；气温等其他信息通过试验站气象设施观察得到。

表1 试验区土壤理化参数范围Tab.1 Soil Physical and chemical parameters in the test area

1.3 Philip入渗模型

1957年，在系统研究Richard方程的基础上，Philip提出了能较好地描述土壤水分一维垂直入渗情况的入渗模型，该模型是半经验半理论性的模型，具有明确的物理意义，入渗公式中只含有两个入渗参数，形式简单，是被人们认可的、表述土壤水分入渗过程的模型之一，在世界范围内得到较广泛的应用[10]。Philip入渗公式如下：

I=St0.5+At,

(1)

式中：I为累积入渗量：自入渗开始至t时刻的入渗总量，cm；S为吸渗率：入渗开始后第一个单位时段末的累积入渗量与稳定入渗率之差，在数值上等于入渗开始至一分钟末的累积入渗量减去稳渗率，cm/min0.5；t为入渗历时，min；A为稳渗率，在单位势梯度下饱和土壤的入渗速度或非饱和土壤入渗达到相对稳定阶段的入渗速度[11]，cm/min。

1.4 建模样本

本文建立了200组由Philip入渗模型参数和盐碱土壤基本理化参数构成的样本，其中197组作为建模训练样本，3组用于检验所建预报模型的精度。

表2 部分建模样本数据Tab.2Partial sample data

2 输入因子和输出因子的确定

2.1 盐碱地土壤理化特性对入渗参数的影响分析

(1)全盐量对入渗参数的影响：盐碱地土壤的可溶性盐浓度若高于作物的耐盐度，则会影响作物的正常生长，轻则作物出苗率低，重则作物生长受到抑制甚至不生长。大量试验资料表明：在土壤结构、质地、含水率、有机质等条件相近的情况下，盐碱地土壤的入渗能力明显小于非盐碱地土壤的入渗能力，根本原因在于：与非盐碱地的水力传导度相比，盐碱地土壤的水力传导度变小[13]，全盐量越高，土体团聚性越差，土粒越分散，结构性越差[14]，且在入渗过程中，土粒随着水流逐渐下移沉积的可能性越大[15]，可能导致大孔隙孔径减小、小孔隙孔径更小甚至被堵塞，水流下渗受到的阻力增大，吸渗率S越小；在单位势梯度下，单位土壤面积通过的流量减小，水力传导度减小，入渗量小，稳渗率A越小。

(2)有机质对入渗参数的影响：有机质具有很好的胶结作用，比砂粒的黏结力大，促进砂土团粒状结构的形成，比黏粒的黏结力小，降低黏粒间的黏结力，使土壤的大、小孔隙分配合理，从而改善了土壤的通透性，增强孔隙的稳定性，改良土壤结构。在其他条件一定的条件下，有机质量低的土壤，较大、较小孔隙居多，水分进入土壤初期，较大孔隙先充满水，并在重力作用下向下层土壤入渗，S值较大；反之，有机质量高的土壤，中等孔隙发育较好，在水分入渗初期，S值较小。稳渗率A的大小主要由孔隙的大小和多少决定，有机质量高的土壤，中等孔隙居多且土壤团粒结构稳定性较好，A值较大；有机质量低的土壤，大孔隙和小孔隙居多且土壤团粒结构稳定性差，水在入渗过程中，土粒容易随着水流移动，使小孔隙孔径减小，甚者堵塞小孔隙，迫使入渗水流改变方向，下渗路径曲折蜿蜒，A值较小。

(3)土壤结构对入渗参数的影响：土壤结构是指土壤的松紧程度、孔隙程度和板结程度，用土壤干容重(土壤相对体积质量)作为反映土壤结构的指标[16]。若土壤干密度大，土壤越密实，土壤孔隙就小，孔隙连通性差，水分入渗路径越弯曲，水流受到的阻力越大，单位势梯度下，水分通量小，土壤的水力传导度就小，则S和A的值越小；反之，越大。

(4)体积含水率对入渗参数的影响：土壤含水率决定着土壤固、液相的比例，常用体积含水率表示。在土壤质地、结构相同或相近的情况下，在入渗开始后的很短时间内，表层土壤很快达到饱和，若土壤含水量大，势梯度就小，单位过水断面的入渗通量就小，S值就小；在土壤结构一定的条件下，含水率越大，孔隙相对减少，稳渗率A就越小。

(5)土壤质地对入渗参数的影响：土壤质地的分类和划分标准有很多种，本文采用国际制对盐碱土壤质地进行划分，并以黏粒含量、粉粒含量和砂粒含量作为数学表征值。在土壤结构、含水率相同的条件下，质地越重，黏粒含量越多，固体相比表面积越大，颗粒表面的吸附能力越大，粒间孔隙孔径越小，在单位势梯度下，入渗水流路径严重曲折，导致水分通量较小，水力传导度低，入渗参数S和A都减小。

2.2 预报模型输出和输入因子的确定

预报模型的输出因子为Philip入渗模型的入渗参数,即吸渗率S和稳渗率A。

预报模型的输入因子根据2.1分析的影响因素初步确定，然后由输入变量的显著性检验决定取舍。在输入因子的初步确定中，考虑了以下三点：

(1)入渗参数受土壤层次的影响。从吸渗率S的概念可知，入渗开始较短时间内，水分还未能到达犁底层，所以影响S的因素仅限于耕作层；而影响A的因素既有耕作层的也有犁底层的。

(2)黏粒、粉粒和砂粒的数学关系：由于黏粒含量(%)+粉粒含量(%)+砂粒含量(%)=100%，且三者对入渗参数的影响高度相关，故本文只选择黏粒含量和粉粒含量作为入渗参数的输入因子[17]。

(3)盐分的集聚性：邹长明等[18]研究认为：盐分主要积聚在0～15 cm土层，尤其是0～5 cm土层。马献发等[19]更是提出：0～10 cm土层的盐分含量远远高于其下各土层。因而，只选择耕作层的全盐量作为输入因子。

那么，初步确定吸渗率S的输入因子为：μ1、μ2、G、ω1、ω2、γ1、γ2、θ1；初步确定稳渗率A的输入因子为：μ1、μ2、G、ω1、ω2、ω*1、ω*2、γ1、γ2、γ3、θ1、θ2。

3 入渗参数预报模型的建立与检验

3.1 建立模型

建立模型的目的是找到一个数学表达式，以测定的土壤常规理化参数为自变量，计算出不易测定的入渗参数。本文基于盐碱地197组田间入渗试验样本，借助MATLAB[12]软件，建立盐碱地Philip入渗模型参数与其土壤理化参数间的关系模型。其具体操作步骤和结果如下：

(1)确定入渗参数与单个输入因子的函数关系。运用控制变量法，选择一个因子作为唯一变量，要求其他因子相同或相近，对数据进行筛选，选出5～10组满足要求的数据，借助MATLAB软件的Cftool功能[12]，输入单个因子及入渗参数的实测值，由离散点的走势，选择或定义函数形式，参考拟合度的大小，最终确定函数的最佳形式。表3即入渗参数与各个输入因子的函数关系。

由表3可知：全盐量μ与入渗参数的关系既包含对数关系又包含线性关系；有机质G与入渗参数呈指数关系；容重γ与入渗参数呈线性关系；含水率θ与入渗参数的关系既包含对数关系又包含线性关系；黏、粉粒含量ω与入渗参数呈线性关系。

表3 入渗参数与单个输入因子的函数形式Tab.3 The function form of infiltration parametersand single input factors

(2)确定回归方程的初步形式。将表2给出的各函数进行机械相加，并添加一个常数项φ，便得出回归方程的初步形式：

S=α1lnμ1+α2μ2+α3eG+α4γ1+α5γ2+

α6lnθ1+β6θ1+α7ω1+α8ω2+φ

(2)

A=α1lnμ1+β1μ1+α2lnμ2+β2μ2+α3eG+

α4γ1+α5γ2+α6γ3+α7lnθ1+α8lnθ2+α9ω1+

α10ω2+α11ω*1+α12ω*2+φ

(3)

(3)回归方程中自变量的显著性检验(T检验)。将样本数据和预报模型的初步函数形式输入到已编好的程序中，运行后可得到初步模型各项输入变量前的回归系数及其进行T检验后的结果，若存在某个|T|≤T0.05/2的值(显著水平α=0.05，T0.05/2=2.014 1)，则剔除掉最小T值所对应的自变量，继续对剩余自变量进行T检验，直至所有自变量对应的|T|值都大于T0.05/2后，停止检验。表4是所有自变量进行T检验的过程和结果。

表4 T检验过程及结果Tab.4 Inspection process and results

注：“___”表示T检验中的最小值。

由表4知，经过4次T检验后，lnμ1、μ2、eG、γ1、lnθ1和ω1对应的|T|都大于T0.05/2，即在入渗参数S的回归方程中显著，那么S最终的输入因子：μ1、μ2、G、γ1、ω1、θ1；经过5次检验后，lnμ1、μ1、lnμ2、eG、γ1、γ2、lnθ1、ω2、lnθ2和ω*1对应的|T|都大于T0.05/2，即在入渗参数A的回归方程中显著，那么A最终的输入因子：μ1、μ2、G、γ1、γ2、θ1、ω2、γ3、θ2、ω*1。

(4)回归方程的显著性检验(F检验)。除了要对方程自变量进行显著性检验，还要对整个回归方程自身进行显著性检验，即F检验，以确保回归方程显著。给定显著水平α=0.05，由样本情况查表可知Fα=0.05值，若计算得到的F值大于Fα=0.05，则方程显著；反之，不显著。F检验结果见表5。

表5 F检验结果Tab.5 F test results

由表5可知：FS和FA都大于Fα=0.05,即入渗参数S和A的回归方程显著，两者的平均误差分别为8.062%、9.043%。则入渗参数的回归方程的最终形式为：

S=0.033 6 lnμ1-0.000 03μ2+0.049eG-

0.637 2γ1-0.038 5 lnθ1-0.009 9ω1+1.265 7

(4)

A=-0.040 6 lnμ1+0.000 03μ1-0.010 4 lnμ2+0.025eG-

0.037 9γ1-0.021 4γ2+0.009 1 lnθ1-0.002 3 lnθ2-

0.000 6ω2-0.002 4ω*1+0.400 3

(5)

3.2 预报模型检验

3.2.1 单参数预测结果检验

表6是3组检验样本(未参与建模的样本)的土壤条件，采用以上所建模型对检验样本的入渗参数进行预报，预报结果见表7。

表6 检验样本的土壤条件Tab.6 Soil conditions for test samples

表7 预报结果Tab.7 Forecast result

由表7可知，预测吸渗率S的相对误差范围是6.977%～8.358%；预测稳渗率A的相对误差范围是6.25%～11.688%，计算I90与实测I90的相对误差介于1.969%～11.331%之间，预报结果的精度在可接受范围。

3.2.2 Philip入渗模型二参数综合检验

将197组建模样本的预测值S和A带入到Phillip入渗公式中，得到对应的90 min累积入渗量I90的计算值，与197组建模样本的实测值I90相比较，两者的拟合结果见图1。

图1 实测值I90与计算值I90拟合效果Fig.1 Fitting effect ofmeasured values and calculated value

由图1所示：197组实测I90与计算I90的拟合程度较高(R2=0.985)，说明由入渗参数S、A计算得到的I90比较接近实测I90，进一步说明所建模型的预报精度较高，模型较可靠。

综上所述，无论是预报模型的单参数预测结果检验，还是Philip入渗模型二参数综合结果检验，都表明预报精度较高。因此，所建立的多元非线性预报模型是可行的。

4 结 语

(1)以常规土壤理化参数为输入因子的多元非线性模型对盐碱地Philip入渗模型参数A和S进行预测，两参数平均误差都小于10%，在建模误差范围之内，预报模型的参数检验误差也在可接受范围之内，故所建模型是可行的，可实现对盐碱地入渗参数的准确预报。研究结果也为盐碱地土壤灌水技术参数的合理确定和盐碱地改良提供了理论与技术支撑。

(2)本文选择对入渗参数影响较大的土壤理化特性作为模型的输入因子，其他如地表温度、地下水矿化度以及土壤酸碱度等次要因素均未考虑在内，因此，在后续的研究中应丰富影响因素，以建立精度更高的预报模型。

□

[1] 李庆梅，侯龙鱼，刘 艳，等. 黄河三角洲盐碱地不同利用方式土壤理化性质[J].中国生态农业学报，2009，17(6)：1 132-1 136.

[2] 张小娟，韩 明，郑敏娜. 不同油用向日葵品种在盐碱地的适应性研究[J].山西农业科学，2012，40(4)：332- 335.

[3] QingjieWanga, CaiyunLua, Hongwen Lia,et al. The effects of no-tillage with subsoiling on soil properties and maize yield: 12-Year experiment on alkaline soils of Northeast China[J]. Soil and Tillage Research，2014，137：43-49.

[4] Chun Xia Wang, Quan Jiu Wang, Zhuang Liang. Field Experimental Study on Saline and Alkaline Soil Amelioration through Combining Chemical Method and Under-Mulch Drip Irrigation. Advanced Materials Research, 2012，599：832-837.

[5] 吕殿青，王全九，王文焰，等. 土壤盐分分布特征评价[J].土壤学报，2002，39(5)：720-725.

[6] 张 莉，丁国栋，王翔宇，等. 夹砂层土壤入渗过程水盐运动规律室内模拟试验研究[J]. 安徽农业科学，2010，38(7)：3 605-3 606.

[7] B B PATEL, BHARAT B PATEL, R S DAVE. Studies on infiltration of saline-alkali soils of several parts of mehsana and patan districts of north gujarat[J]. Journal of Applied Technology in Environmental Sanitation.2011，(1)：87-92.

[8] 任长江，白 丹，周 文，等. 土壤初始含盐量对水分入渗特性的影响[J].干旱区研究，2014，31(2)：222-225.

[9] PHILIP J R. The theory of infiltration about sorptivity and algebraic infiltration equations [J]. Soi Sci.，1957，84(4)：257-264.

[10] 原林虎. PHILIP入渗模型参数预报模型研究与应用[D]. 太原：太原理工大学，2013：6-8.

[11] 岳海晶. 樊贵盛.备耕头水地土壤入渗参数的线性预报模型[J]. 中国农村水利水电，2016,(2)：23.

[12] 汪岚，黄彩虹. 基于MATLAB色差预测多元回归模型的研究[J]. 计算机与应用化学，2008，25(8)：1 015-1 018．

[13] 樊贵盛，李 尧，苏冬阳，等. 大田原生盐碱荒地入渗特性的试验[J].农业工程学报. 2012，28(19)：63-70.

[14] So H B,，AYLMORE L A G. How do sodic soils behave? The effects of sadicity on soil physical behavior [J]. Aust. J. Soil Res. 1993，31.

[15] 王 雪，樊贵盛. Na+含量对土壤入渗能力影响的试验研究[J].太原理工大学学报，2009，40(4)：391-393.

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

[17] 郭 华. 土壤入渗模型参数的分阶段非线性预报模型研究[D].太原：太原理工大学，2016：46.

[18] 邹长明，张多姝，张晓红，等．蚌埠地区设施土壤酸化与盐渍化状况测定与评价[J]．安徽农学通报，2006，12(9):54-55．

[19] 马献发，白路平，张继舟.哈尔滨市郊设施土壤积盐规律的研究[J]．黑龙江农业科学，2009，11(1)：50-52．