Philip入渗模型参数的非线性预报模型

郭 华，樊贵盛

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

土壤水分入渗过程[1]是指降水或灌溉水通过地表进入土壤的过程，是田间水循环的重要组成部分，也是灌水技术参数优化、解决地下水补给等问题的重要依据。为了反映土壤水分入渗过程，国、内外学者提出了许多入渗模型[2-4]，如Horton模型、Smith模型、蒋定生公式等，其中Horton模型和Kostiakov模型是早期提出的经验型模型，但没有明确的物理意义；Green-Ampt模型有明确的物理基础，但对参数精度要求低；而Philip入渗模型不仅有明确的物理基础，对参数精度要求高，而且能较好的反映土壤水分入渗过程。邵明安等人[5]曾用Philip模型预测了土壤水分运动参数，武世量[6]研究了Philip公式标定因子与土壤理化参数的关系，任尚岗等人[7]基于Philip公式分析了单环入渗的三维动态过程，原林虎[8]采用Philip模型描述了水分入渗过程，但目前而言，对Philip入渗模型参数的预报尚鲜见。

土壤水分入渗参数[9]是表征土壤入渗能力的指标，直接影响着大田土壤的灌水质量和灌溉效果，然而由于土壤条件空间变异性和其他限制因素的存在，很难通过实测获得大量的入渗参数。基于此，学者们通过土壤传递函数法[10](Pedo-Transfer Functions，PTFs法)，即建立试验较易获得的土壤理化参数(如土壤含水率、干密度、质地等)与入渗参数之间的函数关系，来间接获得入渗参数，以解决获取入渗参数困难这一问题。其中冯锦萍[11]用线性回归分析法对入渗参数进行了预报，其方法简单，但预测精度低；原林虎[8]采用人工神经网络法对入渗参数进行了预测，其精度较高，但易出现局部极小等不良现象，不能保证全局最优。因此，本研究基于土壤传递函数法，以黄土高原区大田耕作土壤为试验对象，将试验较易获得的土壤体积含水率、干密度、粉粒含量、黏粒含量和有机质含量等土壤理化参数作为输入变量，建立土壤水分入渗参数为输出变量的土壤传递函数，以实现对入渗参数的多元非线性预报，从而提高预测精度，达到全局最优。

1 材料与方法

1.1 试验条件

试验对黄土高原区大田耕作土壤进行，主要在大同、吕梁、应县等地进行。大同位于山西最北端，平均海拔700～1 400 m，干寒多风，温差较大，由于多期地貌构造的变动，形成了山地、丘陵、平川等地形；吕梁地处山西中部西侧，平均海拔1 000～2 000 m，四季分明，降雨量分布不均，由于吕梁山脉的存在，形成了由南到北，由高山逐渐变为平川的地势形态；应县位于山西北部的桑干河中游，海拔在1 000～2 300 m之间，全年气候偏低，其黄土丘陵区主要集中于东北、西北地带，其余为平川地区。试验区均位于山西境内，是被黄土广泛覆盖的山地形高原，属大陆性半干旱季风气候，在自然作用和长期农业活动的共同影响下，土壤条件在空间上发生很大变化，形成了复杂多样的土壤类型。试验区有壤土、沙壤土、黏壤土、粉沙质黏壤土、粉沙壤土、黏沙壤土、轻沙壤土等土壤类型，其中耕作层土壤的沙粒含量为9.7%～88%，粉粒含量为10.88%～76%，黏粒含量为0.84%～23.56%；犁底层土壤沙粒、粉粒、黏粒含量分别变化在9.66%～93%、5.41%～85%、0.91%～33.15%。耕作层土壤体积含水率为2.08%～70.67%，有机质含量为0.416～37.54 g/kg，其中地表以下0～10 cm土层的干密度的范围是0.854～1.426 g/cm3，地表以下10～20 cm土壤干密度变化在0.891～1.523 g/cm3；犁底层土壤体积含水率为3.03%～75.87%，其土壤干密度为0.774～1.762 g/cm3。

1.2 试验设备与方法

试验在非冻融期进行。试验前将内、外环直径分别为26、64.4 cm，内、外环高度均为25 cm的双套环入渗仪埋于试验地块，试验时对内、外环均进行灌水，并用自制的水位控制装置将内、外环水位控制在2～3 cm，以消除积水水头对势梯度的影响；根据Philip入渗公式，通过定量描述累积入渗水量与灌水时间的关系，便可间接获得土壤水分入渗参数。大量试验资料表明，土壤水分达到稳渗状态的时间在60 min，考虑到试验资料的可靠性和一致性，将水分稳定入渗时间定为90 min。

在大田土壤水分入渗试验进行的同时，还需对试验地块的土壤理化参数进行测定。土壤理化参数包括土壤体积含水率、干密度、质地以及有机质，其中土壤体积含水率采用烘干、称重的方法测定；土壤干密度用环刀法获得；通过筛分曲线法分析各个土层的颗粒级配，从而获得各层土壤质地；土壤有机质主要存在于耕作层(0～20 cm)土壤中，通过重铬酸钾滴定法获得。

1.3 Philip入渗模型

在对Richard方程系统研究的基础上，Philip[12]提出了适合土壤一维垂直入渗的入渗模型，该模型形式简单，计算方便，参数容易确定，具有明确的物理意义，且能较好反映土壤水分入渗过程，在多年的土壤入渗研究中应用广泛，其具体表达式为：

I=St0.5+At

(1)

式中：I为单位面积上的累积入渗量，cm；t为入渗历时，min；S为吸湿率，在数值上等于第一个单位时段末的累积入渗量，决定着入渗初期入渗率的大小，cm/min0.5；A为稳定入渗率,接近土壤的水力传导度， cm/min。

2 土壤入渗参数的预报模型

2.1 输出因子与输入因子

随着农业用水的浪费、灌溉水利用率低、土壤环境退化等问题的出现，如何合理灌溉成为需要急切解决的问题。为此，学者们提出了土壤水分入渗模型，其中土壤水分入渗参数作为地面灌溉水流运动模拟模型的输入参数，其取值的合适与否，直接决定着灌水质量和灌溉效果，但通过实测获取大量的入渗参数相当困难，若建立以试验较易获得的土壤理化参数为输入因子，以入渗参数为输出因子的土壤传递函数，便可方便获得入渗参数S和A。因此，本文着力研究Philip入渗模型参数S和A的预测问题，即模型的输出因子为入渗参数S和A。

S在数值上等于第一个单位时段末的累积入渗量，只与表层土壤的理化性质有关；而A主要取决于土壤结构状态，即土壤孔隙的数量和连通性，受到整个湿润土层理化参数的影响。相关研究[8]表明，影响入渗参数的主要因素为土壤体积含水率、干密度、质地和有机质含量，土壤体积含水率越高的土壤，其水吸力越小，土水势越大，即水势差越小，土壤达到饱和的速度越快，S值越小；土壤结构以干密度来表征，能够反映土壤的孔隙率和密实度，土壤越密实，干密度越大，土壤孔隙越小，连通性越差，入渗通量也越小，入渗参数S和A随之变小；土壤质地反映不同大小矿物颗粒的组成状况，其中以黏粒含量作为划分土壤类型的标准，土壤质地越重，黏粒含量越高，土壤孔隙度越小，孔隙越弯曲，连通性越差，对水分入渗的阻力越大，导致水力传导度越小，并且固体颗粒表面积相应增大、吸附能力增强，S和A随之变小；土壤有机质是由动、植物残体分解合成的产物，可促进土壤团聚体的形成，增强土壤的稳定性，有机质含量少的土壤，其结构不稳定，孔隙不均匀，大孔隙相对较多，使得入渗初期的入渗通量较大，随着水分的入渗，裂隙两边的土壤发生崩塌、挤压，堵塞土壤中的大孔隙，入渗路径严重弯曲，导致相对稳定入渗率减小，因而，有机质含量少的土壤，吸湿率S值较大，相对稳定入渗率A较小。

由上述分析可知，S主要受耕作层土壤物理性质的影响，A主要受耕作层和犁底层土壤理化性质的共同影响。土壤体积含水率以θ0(0～20 cm)、θ1(20～40 cm)作为指标；土壤结构以干密度γ0(0～10 cm)、γ1(10～20 cm)、γ2(20～40 cm)作为指标；考虑到土壤沙粒含量、粉粒含量和黏粒含量之间的线性关系，将耕作层(0～20 cm)土壤的粉粒含量w1(0.05～0.002 mm)、黏粒含量w2(<0.002 mm)以及犁底层(20～40 cm)土壤的粉粒含量w3、黏粒含量w4作为土壤质地指标；土壤有机质以耕作层土壤的有机质含量G作为指标。因此，S的回归模型以土壤理化参数θ0、γ0、γ1、w1、w2、G作为输入因子，A的回归模型以土壤理化参数θ0、θ1、γ0、γ1、γ2、w1、w2、w3、w4、G作为输入因子。试验区部分土壤条件如表1所示。

表1 试验区部分土壤条件Tab.1 Some soil condition of test area

2.2 模型的建立与分析

回归分析是建立实测数据间的相互依赖关系，分析数据的内在规律，以用于预报、控制等问题，其中非线性回归分析难度相对较大，为此，MATLAB语言得到广泛应用。MATLAB软件[13]是高级的数值分析、处理与计算软件，具有优秀的计算能力和卓越的数据化可视功能，在矩阵运算方面有其他语言难以比拟的优越性，可节省大量的编程时间。土壤入渗参数与各理化参数之间的关系本质上是非线性的，由于入渗参数受到多个理化参数的共同影响，用多个理化参数的最优组合共同来预测入渗参数相比只用一个理化参数更符合实际。因此，本研究以MATLAB软件为计算工具，建立土壤理化参数与入渗参数之间的多元非线性传递函数，对试验区100组样本进行回归，以更加精准的预测入渗参数，使预测模型更符合实际 。

(1)单因素函数形式的确定。基于大量的实测资料，在其他因素相同时，选取多组单因素试验样本，采用MATLAB中的Cftool工具，绘制散点图和回归曲线图，根据曲线走势，确定单因素最优回归方程，如表2所示。

由单因素最优回归方程的形式可知：在其他因子相同或接近相同时，入渗参数与体积含水率既存在对数关系，又存在线性关系，干密度与入渗参数呈线性关系，耕作层和犁底层土壤的粉粒、黏粒含量与入渗参数均呈对数关系，耕作层土壤的有机质含量与入渗参数为对数关系。

表2 入渗参数与各土壤理化参数间的单一定量关系Tab.2 Single quantitative relationship between each infiltration parameter and soil physicochemical parameters

(2)多元非线性函数的确定。基于对各影响因素与入渗参数内在关系的研究，确定单因素回归方程，以各单因素回归方程中的因子作为回归模型的输入因子，建立入渗参数的多元非线性回归模型，并对各因子进行T检验，剔除每次检验中T值最小对应的因子，直至所有因子∣T∣≥T0.05/2为止(T0.05/2=2.014 4)，检验结果如表3所示。

表3 S回归方程和A回归方程自变量的确定Tab.3 Independent variables' determination of S and A regression equations

注：*表示每次进行T检验时的最小值，将其对应的因子剔除，不作为回归方程的自变量。

对回归方程输入因子的检验结果可知，lnw2在S的回归方程中不显著，lnθ0、θ1、lnw3在A的回归方程中不显著。最终确定入渗参数回归方程的自变量，如下所示：S回归方程的自变量为lnθ0、θ0、γ0、γ1、lnw1、lnG；A回归方程的自变量为θ0、lnθ1、γ0、γ1、γ2、lnw1、lnw2、lnw4、lnG。

基于土壤传递函数----非线性回归分析法，以上述自变量为输入参数，建立以S和A为输出参数的多元非线性回归方程，如下所示：

S=1.519 7-0.092 9 lnθ0+0.001 1θ0-0.635γ0+

0.040 8γ1+0.050 9 lnω1+0.140 7 lnG

(2)

A=0.321 5+0.000 2θ0-0.029 2 lnθ1-0.103 5γ0-

0.024 4γ1+0.036 4γ2-0.013 9 lnω1+0.015 4 lnω2-

0.012 lnω4+0.018 2 lnG

(3)

(3)灰色关联分析。为反映自变量与预测变量的关联程度，以Philip模型参数S和A为参考序列，各自变量为比较序列，依据灰色关联理论[14,15]，对自变量进行排序，排序结果见表4。

表4 自变量的排序结果Tab.4 Sorted Results of independent variables

由表4可知，不同回归模型中，自变量与预测值之间的关联度不同，S回归模型中自变量的关联度顺序为θ0>lnG>γ1>lnw1>γ0> lnθ0，A回归模型中自变量的关联度顺序为lnw2>γ1> lnG>θ0>lnw1>γ0>γ2> lnθ1>lnw4。

(4)模型检验。在自变量检验的基础上，对回归方程进行显著性、误差以及拟合度分析，以保证方程整体的可靠性，分析结果如表5所示。

模型分析结果显示：S和A回归模型的F值均大于F0.05、平均相对误差分别为7.53%和7.64%、负相关系数分别为0.882、0.933，累积入渗量I的综合平均相对误差为7.81%，则入渗参数S和A的多元非线性预测模型显著性强、精度高、拟合度高。

表5 模型分析Tab.5 Model analysis

注：R为复相关系数，R越接近1，表明方程的拟合度越高。

2.3 预报实例

利用以下2组样本，采用上述100组样本建立的非线性预测模型对Philip模型参数进行预报，预报样本的土壤条件如表6所示，预报结果如表7所示。

表6 预报样本的土壤条件Tab.6 Soil condition of forecasting samples

表7 非线性预报结果Tab.7 Nonlinear forecasting results

由非线性预报结果可知，入渗参数S和A的预测相对误差分别变化在5%～4.63%、5.46%～7.81%，累积入渗量I的综合影响相对误差为3.42%～4.32%，以上所得精度均可满足非冻融条件下Philip模型参数的预报要求。预报结果表明用土壤体积含水率、干密度、粉粒含量、黏粒含量以及有机质含量对Philip模型参数进行非线性预报具有一定的可靠性。

3 结 语

(1)不同土壤条件下，用土壤理化参数对Philip入渗模型参数进行预报是可行的，即采用土壤传递函数----非线性回归分析法，以土壤理化参数作为输入变量，对Philip入渗模型参数S、A进行预报是可行的，其综合相对误差可控制在的8%以下，实现了对入渗参数的高度预测，解决了获取入渗参数困难的问题。

(2)不同入渗参数的输入变量不同，入渗参数S回归模型的最终输入变量为耕作层土壤(0～20 cm)的土壤体积含水率、干密度、地表以下0～10 cm土壤的干密度、粉粒含量和有机质含量；入渗参数A回归模型的最终输入变量为 耕作层土壤的土壤体积含水率、干密度、地表以下0～10cm土壤的干密度、粉粒含量、黏粒含量、有机质含量以及犁底层土壤的土壤体积含水率、干密度、黏粒含量。

本文所建立的Philip入渗模型参数的非线性预报模型，体现了非冻融期土壤理化参数与入渗参数之间的非线性本质，较好地反映了土壤水分入渗过程，克服了预测精度低、不能保证全局最优等缺陷，但Philip模型仅适用均质一维垂直入渗的情况，还需进一步研究非均质土壤的入渗过程。

[1] 邵明安，王全九，黄明斌. 土壤物理学[M]. 北京：高等教育出版社，2006.

[2] 王全九，来剑斌，李 毅. Green-Ampt模型与Philip入渗模型的对比分析[J]. 农业工程学报， 2002,18(2):13-16.

[3] Green W H, Ampt G A.Studies on soil physics, flow of air and water through soils [J].J. Agr. Sci.，1911,76(4)：1-24.

[4] Horton RE．An approach toward a physical interpretation of infiltration-capacity[J]．Soil Sci., 1957,84(4):257-264．

[5] 邵明安，王全九. 推求土壤水分运动参数的简单入渗法[J]. 土壤学报， 2000, 37(1):1-8.

[6] 武世亮，聂卫波，马孝义. 土壤特性与Philip入渗公式标定因子的空间变异性关系[J]. 排灌机械工程学报， 2014,32(8):730-736.

[7] 任尚岗，张振华，杨润亚，等. Philip公式在三维入渗及参数测算中的应用[J]. 干旱地区农业研究， 2011,29(5):192-196.

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

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

[10] 朱安宁，张佳宝，陈效民，等．封丘地区土壤传递函数的研究[J]. 土壤学报， 2003,40(1):54-59.

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

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

[13] 彭 勇，牛俊峰. 用MATLAB回归非线性模型参数[J]. 化学通报， 2007,(11):880-883.

[14] 邓聚龙.灰色系统控制系统[M].武汉：华中理工大学出版社， 1993.

[15] 刘思峰，党耀国，方志耕.灰色系统理论及其应用[M].北京：科学出版社， 2004.