非线性入渗模型求参中参数初始值确定方法研究

沈洪政，王仰仁，韩娜娜

(天津农学院水利工程学院，天津 300384)

0 引 言

参数拟合初始值的确定是非线性模型参数优化的基本问题。对于非线性模型参数拟合人们提出了多种方法[1-3]。李春友[4]等提出用优化方法计算Van Genuchten模型参数，拟合初始值参考有关文献选取。由于所选初值有时对拟合值有较大的影响,因此在这种情况下就需要多次调初值才可获得较为满意的拟合效果。钱天伟[5]等用单纯形调优法和单纯形加速法作为优化手段，结合计算机程序进行参数优化求解，但其初始值是任意给定的。以上均是从参数优化求解方法角度进行研究，对于优化初始值如何选取，没有给出普遍适用的方法。

在求解非线性模型问题时，通常利用最小二乘法进行参数优化拟合。对于有明确物理意义的参数，可以根据其物理意义进行初始值的确定。对于没有明确物理意义的模型参数，通常采用均匀网格法[6]和外推内插法[7]等进行初始值的选取。对于以上方法，由于需要进行反复寻优，选取的初始点较多，过程复杂，使用不方便。因此，本文以拟合精度较高的非线性土壤入渗Green Ampt模型和Horton公式为例，提出采用泰勒级数法确定初始值，避免了初始值选取的烦琐过程，提高了参数拟合的效率和精度。

1 方法介绍

1.1 目标函数

(1)

1.2 入渗模型介绍

1.2.1 Green-Ampt模型

Green-Ampt模型是Green W H和Ampt G A根据毛管理论提出的近似积水模型[8,9]。在土层均匀、初始土壤干燥的薄层积水入渗条件下，利用达西定律，经过化简与推导，可得到Green-Ampt入渗公式[10]：

(3)

式中：i(t)为入渗率，mm/min；ic为稳定入渗率，mm/min；I(t)为累积入渗量，mm；t为入渗时间，min；b为参数，其中ic和b为待定参数。

1.2.2 霍顿(Horton)公式

霍顿(1940年)提出的入渗公式为：

i(t)=ic+(i0-ic)e-βt

(4)

(5)

式中：i0为初始入渗率，cm/min；β为衰减参数，其余符号意义同前，其中i0、ic、β为待定参数。

Green-Ampt模型和Horton公式具有较高的拟合精度，因而被得到广泛的使用。由公式的形式易知，两者均为非线性模型，其参数的拟合属于非线性拟合问题。在利用最小二乘法(OLS)进行参数优化拟合时，两模型中的参数i0和ic的初始值可根据其物理意义进行确定，而参数b及β由于没有明确的物理意义，用外推内插法等方法确定时过程烦琐。针对此问题，本研究提出利用泰勒级数法对公式进行展开，结合田间实测入渗数据，获得未知参数b和β的初始值，再利用OLS进行拟合，即可得到最优入渗模型参数。

1.3 初始值的确定

1.3.1 参数i0和ic初始值的确定

参数i0和ic分别表示初始入渗率和稳渗率。初始入渗率可选取入渗试验测定的最大入渗率值作为参数i0的优化初始值，稳渗率可选取入渗试验测定的最小入渗率值作为参数ic的优化初始值。

1.3.2 参数b初始值的确定

参数b初始值可采用泰勒级数法确定。为了便于泰勒公式的展开，对于公式(3)将其改写为以t为因变量的显式函数：

(6)

(7)

取式(7)中前3项，代入式(6)中得：

(8)

1.3.3 参数β初始值的确定

(9)

取式(9)中前3项，代入式(5)中得：

(10)

由式(8)和式(10)知，通过泰勒级数法展开后，非线性模型转化为线性模型，结合i0和ic的实际物理意义及田间实测入渗数据即可确定入渗模型参数的初始值，进而可以利用最小二乘法进行参数拟合。现以天津农学院农田水循环试验基地测得的田间入渗数据为例，验证方法的合理性。

2 实例计算

2.1 试验区概况

天津农学院农田水循环试验基地位于天津市杨柳青镇大柳滩村(116°57′E，39°08′N，海拔5.49 m)，全年平均气温11.6 ℃。试验田总面积1 hm2，地下水位变幅在2.06～3.70 m之间，试验区土壤理化性质见表1，土壤水分特征曲线参数[11]见表2和图1。

2.2 测试方法

试验在2011年分2次测试，分别在灌溉冬水的麦田、不灌冬水的麦田和棉花茬地进行(试验区简图见图2)。测试时间为1月6日(气温为-5.6 ℃，相应麦田测试点次编号分别记为1和2)和3月28日(气温为11.2 ℃，麦田测试点次编号分别记为4和5)；棉花试验田只在3月28日测试，其测试点次编号记为3。

表1 试验区土壤理化性质

表2 试验区土壤水分特征曲线参数

图1 试验区土壤水分特征曲线图

试验采用双环入渗仪测定土壤水分入渗过程[12]。内环横截面积1 000 cm2，内环直径35.7 cm，外环直径60 cm，内外环高度均为30 cm。试验前，将渗透环埋于试验田，埋入土壤的深度为10 cm。内环加水保持5 cm水头，采用手动方式向环内加水，利用量筒计量，使得内外环水位齐平。试验开始后记录每次灌入水量和相应的时间，其内环灌入水量即为渗入水总量，测试直到入渗率基本稳定为止，记录内环水量全部渗完为止的时间。

图2 试验区简图(单位：m)

2.3 结果分析

2.3.1 b值与β值计算结果分析

根据初始值确定方法的介绍，以天津农学院水循环试验基地在2011年1月6日和2011年3月28日测得的入渗数据为基础，经过计算，给出了各测试点b和β平均值的计算结果(表3)以及各测试点b和β值随时间变化的曲线关系图(图3、图4)。

表3 各测试点b与β平均值计算结果表

由表3知，通过平均值法计算得到的Green-Ampt模型中参数b的值在各测试点之间变化相对较大，Horton公式中参数β的值在个测试点之间变化较小。从Green-Ampt模型和Horton公式的实际意义分析来看，当入渗时间t→∞时，入渗速率i趋于稳定入渗率ic，当入渗时间t不断增加时，入渗速率i也趋于稳定入渗率ic。为满足此变化规律，要求b值和β值均为正值，这与计算的结果相符，表明用取平均值的方法计算b值和β值是可行的。

图3 各测试点b值随时间变化曲线

图4 各测试点β值随时间变化曲线

由图3和图4知，在入渗试验初期，由于初始表层土壤含水率及土壤吸力大小等变化较大，入渗速率变化较大，导致b值和β值变化较大。随着入渗试验的进行，入渗速率逐渐稳定，各测试点b值和β值变化也趋于平稳。β值基本上在横轴附近上下变化，b值基本上在79上下变化。说明，采用泰勒级数法确定的b值和β值具有较好的稳定性，且β值较b值更稳定。

2.3.2 Green-Ampt模型及Horton公式参数拟合结果

以b和β求出的值作为初始值，以累积入渗量I(t)的实测值与模拟值误差平方和最小为目标，利用Excel规划求解工具，分别对Green-Ampt模型及Horton公式的参数进行拟合。

表4 Horton公式参数拟合结果表

表5 Green-Ampt公式参数拟合结果表

由表4和表5知：

(1)采用泰勒级数法确定初始值来描述入渗变化具有较高精度，相关系数R2均达到0.85以上，最大值达到了0.998 6，具有较好的相关性。但由于初始入渗率i0和稳渗率ic只是简单选取测试数据的最大值和最小值，同时β和b值均采用测试数据的平均值计算得出，因此，需对参数初始值进行拟合。由拟合后的结果知，相关系数R2均达到了0.99以上，拟合精度得到显著提高。

(2)对比Horton公式和Green-Ampt模型参数拟合结果知，Horton公式拟合参数均为正值，Green-Ampt模型拟合结果中稳定入渗率 出现了负值，与稳定入渗率 的实际物理意义不符。从均方差计算结果来看，Horton公式计算的各测试点之间均方差的值相差不大，均在1.63上下变化。Green-Ampt模型计算的各测试点之间均方差的值相差较大，最大值达到了7.31，最小值为0.76，变化较大。因此，Horton公式的参数比Green-Ampt模型的参数更加稳定，公式更具有适用性。

图5 Green-Ampt模型各测试点累积入渗量实测值与模拟值随时间变化曲线

图6 Horton公式各测试点累积入渗量实测值与模拟值随时间变化曲线

图5和图6给出了各测试点Green-Ampt模型和Horton公式累积入渗量模拟值与对应实测值随时间变化的曲线。由图可见，各测试点Green-Ampt模型和Horton公式累积入渗量的模拟值拟合精度均较高。说明，采用泰勒级数法确定的初始值进行参数拟合具有较高的准确性。另一方面，由于测试时温度的不同，灌冬水的麦田测试点1和测试点5，其累积入渗量之间变化较大，不灌冬水的麦田测试点2和测试点4也有相同变化规律。而对于测试温度相同的麦田测试点1、2(-5.6 ℃)和测试点4、5(11.2 ℃)而言，其累积入渗量实测值与模拟值随时间变化趋势相近，表明，温度对土壤水分入渗过程有很大影响。

2.3.3 温度对入渗结果的影响

通过图5和图6的分析结果可知，温度对入渗结果有很大影响，温度越高，入渗速率越大，反之越小。兹以参数稳定性较好的Horton公式为例，考虑温度对入渗结果的影响，在式(4)和式(5)中引入温度修正系数，

i(t)=KT[ic+(i0-ic)e-βt]

(12)

其中：

KT=ea(T-15)

(13)

式中：KT为温度修正系数；T为气温，℃；a为参数，15为参考温度，℃。对式(11)，将各测试点(1和5、2和4)在同一时间的实测入渗量数值之比与不同温度(-5.6 ℃，11.2 ℃)下的温度修正系数KT值之比建立等式关系[式(14)]，由此便可求得参数a的值(见表6)，由测试点1和5求得的a值为0.044 1，由测试点2和4求得a的值为0.050 7。参数a的值确定以后，再将不同温度(-5.6 ℃，11.2 ℃)值代入式(11)中，即可求得温度修正系数值。

(14)

式中：T1、T2分别为不同灌水地块的气温；I1(t)、I2(t)分别为不同测试点的实测入渗量值。图7和图8分别表示不考虑温度修正的入渗曲线及考虑温度修正的入渗曲线。由图7和图8知，引入温度修正系数以后，不同温度的入渗曲线很好的重合到一起(参考温度条件下)，显著地改善了入渗曲线拟合精度。

表6 a值计算结果表

图7 不考虑温度修正的入渗曲线

图8 考虑温度修正的入渗曲线

3 结 语

本文以非线性土壤入渗Green-Ampt模型和Horton公式为例，提出采用泰勒级数法确定其模型参数初始值，同时考虑了温度对入渗结果的影响。结果表明，在初始值拟合基础上，相关系数均达到0.85以上，经最小二乘法拟合后，实测值与模拟值相关系数提高到0.99以上。考虑温度修正系数以后，不同温度的入渗曲线很好的重合到一起(参考温度条件下)，显著地改善了入渗曲线拟合精度。

泰勒级数法具有更好的理论基础，避免了非线性规划求参中需要设置多个初始值的不足，可显著提高计算效率，而且可以提高参数寻优的准确度，适用于无法线性化的非线性模型参数初始值的确定。但是，在拟合计算过程中Green-Ampt模型参数 值出现了负值情况，与其实际物理意义相矛盾，尚需进一步研究。

□

[1] 朱良君，张光辉，任宗萍. 4种土壤入渗测定方法的比较[J]. 水土保持学报，2012，32(6):163-168.

[2] Ravindra V. Kale & Bhabagrahi Sahoo. Green-Ampt Infiltration Models for Varied Field Conditions: A Revisit[J]. Water Resource Manage，2011,25:3 505-3 536.

[3] 张新民，王根绪，胡想全，等. 优化技术在土壤入渗参数计算中的应用[J]. 灌溉排水学报，2004，(6): 52-54.

[4] 李春友，任 理，李保国. 利用优化方法求算Van Genuchten方程参数[J]. 水科学进展，2001，(12): 473-476.

[5] 钱天伟，陈繁荣，杜晓丽，等. 一种推求van Genuchten方程参数的高性能优化方法[J]. 土壤学报，2004， (11): 973-975.

[6] 尚松浩. 水资源系统分析方法及应用[M]. 北京：清华大学出版社，2012:179-180.

[7] 张 莹. 运筹学基础[M]. 北京：清华大学出版社，2004:135-137.

[8] Green W H，Ampt G A. Studies on soil physics:I. Flow of air and water through soils[J]. AgricSci，1911，(4):1-24.

[9] 李 毅，王全九，邵明安，等. Green-Ampt入渗模型及其应用[J]. 西北农林科技大学学报，2007，(2):225-230.

[10] 雷志栋，杨诗秀，谢森传. 土壤水动力学[M]. 北京：清华大学出版社，1989:130-131.

[11] 车 政，王仰仁，王永红，等. 农田土壤水分特征曲线参数拟合及其剖面变异特性研究[J]. 灌溉排水学报，2016:1-5.

[12] 周鑫洋，王仰仁，刘群昌. 四个入渗公式对农田潮土的适应性研究[J]. 灌溉排水学报，2015，(8):96-100.