估算畦灌土壤入渗参数的线性回归法

估算畦灌土壤入渗参数的线性回归法王维汉1,2，缴锡云1，彭世彰1，马海燕1,2（1. 河海大学水文水资源与水利工程科学国家重点实验室，南京 210098；2. 河海大学现代农业工程系，南京 210098)摘要：将畦灌地表水深与水流推进距离按指数函数进行非线性最小二乘拟合，并根据水量平衡原理提出了估算土壤入渗参数的线性回归法。

实例计算表明：地表水深与水流推进距离之间呈现较好的指数函数关系；线性回归法估算土壤入渗参数计算工作量较小，计算精度较高。

关键词：入渗参数；畦灌；地表水深；线性回归法土壤入渗特性是影响地面灌溉过程的一个十分重要的因素。

土壤入渗参数的估算是地面灌溉研究中的一个重要内容。

早在1956年，Haise就提出了利用筒测仪来测量土壤入渗参数，但由于入渗过程的时空变异性很大[1-2]，土壤入渗参数常常难以准确估算。

国内外学者对估算土壤入渗参数进行了大量研究，提出了多种计算方法。

Elliott和Walker[3](1982)针对沟灌提出了估算土壤入渗参数的两点法，需要分别观测水流前锋推进到沟长中点和沟末端的时间及沟首过流断面面积，在水量平衡的基础上估算入渗参数。

两点法需要观测的数据少，计算简单，但计算结果精度往往不够。

Maheshwari[4](1988)首次将优化技术应用到土壤入渗参数的估算中，通过测量水流推进过程和畦首地表水深的变化过程来计算土壤入渗参数。

这种方法(以下简称M法)适用于任何入渗模型，精度较高，使用较为广泛，但计算量偏大。

Shepard[5](1993)提出了估算土壤入渗参数的一点法，该方法需要测量水流推进到沟末端的时间和沟中的平均过水面积来计算土壤入渗参数，但它只适用于Philip入渗模型，从而其应用受到一定的限制。

Esfandiari[6](1997)对M法进行改进，提出了利用水流推进资料和沿沟长若干点地表水深资料采用模式搜索技术来估算土壤入渗参数的方法，但同样是存在计算量较大的问题。

国内学者也提出了许多估算土壤入渗参数的计算方法。

王文焰[7](1993)提出了利用两个畦田的水流推进消退过程来估算土壤入渗参数的方法。

费良军[8](1999)提出了利用畦灌水流的地表水深资料及水流推进过程来估算土壤入渗参数的方法。

这两种方法都需要至少观测两个畦田的灌水资料才能估算土壤入渗参数，精度较高，但试验工作量较大，且不便于评价入渗模型的合理性。

缴锡云[9](2001)对M法的计算方法进行了改进(以下简称M-J法)，张新民[10](2005)又对M法的计算方法做了改进(以下简称M-Z法),这两种改进减少了一定的试验工作量，但计算量较大的问题仍算存在。

综上所述，在以上估算土壤入渗参数的方法中，两点法、M法等采用地表储水形状系数计算地表储水量，但当畦田长度较长时(大于100m)，地表储水形状系数变化较大，会给计算带来较大误差[11]。

本文依据水位传感器观测的地表水深资料，规避地表储水形状系数，提出计算精度较高的线性回归法，来估算土壤入渗参数。

1理论分析1.1水量平衡方程式建立x，地表水面线与入渗水量分布曲线如图1所示。

在畦灌地表水流推进过程中，对应于推进距离a基金项目：河海大学水文水资源与水利工程科学国家重点实验室开放研究基金资助项目（2005405811）；河海大学科技创新基金资助项目（2084－40401105）作者简介：王维汉（1981—），男，河南南阳市人，在读博士研究生，从事节水灌溉理论与技术研究。

通讯作者：缴锡云（1962—），男，河北文安县人，博士，教授，主要从事节水灌溉理论与技术研究。

根据水量平衡原理有z h V V qt += (1)式中：q 为入畦的单宽流量，L 2·T -1；t 为灌水历时，T ；h V 为水流前锋推进至a x 处时，单位畦宽上的地表储水量，L 2；z V 为对应的单位畦宽上的入渗水量, L 2。

1.1.1 地表储水量在畦灌水流推进过程中, 假定地表水深h 与水流推进距离a x 符合如下的指数函数关系式)1(0s h a b e h h h −+= (2)式中：h 为地表水深，L；0h 为畦首地表水深，L ；a h 、b h 为拟合参数; s 为沿畦长方向上距畦首的距离，L。

则地表储水量可表示为h V =∫ax hds 0ds e h h ab x s h a ])1([00∫−+=bx h a a a h e h x h h a b )1()(0−++= (3)式中：h V 为单位畦田宽度的地表储水量，L 2；a x 为水流推进距离，L。

1.1.2 地下储水量地下储水量的计算涉及到土壤入渗模型，常用的入渗模型有Kostiakov 模型、Kostiakov-Lewis 模型、Green-Ampt 模型、Philip 模型、Horton 模型等。

在实际计算中Kostiakov 经验模型应用较为广泛[12]。

本文采用Kostiakov 入渗模型ατK Z = (4)式中：Z 为单位畦田面积上的累积入渗量，以水深表示，L ；K 为入渗系数，L·T -α；α为入渗指数，无量纲；τ为入渗历时，s t t −=τ，T ；s t 为水流前锋推进时间，T ；假定水流前锋推进过程符合幂函数规律[13]，即rs ps t = ，p 和r 为拟合参数。

于是，地下储水量可表示为ds ps t K V ax r z ∫−=)(α (5)根据Fok 和Bishop 对式(5)的推导[13]，上式可表示为z V =a x Kt r r αααα)1)(1()1(11[++−++ (6) 1.1.3 水量平衡方程式综合式(1)、式(3)和式(6)有Fig .1 Sketch of water surface profiles and infiltration waterdistribution for border irrigation地表水深ha b x h a a a x Kt r r h e h x h h qt a b αααα])1)(1()1(11[)1()(0++−+++−++= (7) 1.2 公式推导及计算方法对式(7)移项，令bx h a a a z h e h x h h qt V a b )1()(0−−+−= (8) =)(αf )1)(1()1(11r r ++−++ααα (9) 则式(7)可简化为ααKt x f V az=)( (10)式中：z V 表示地下入渗水量的观测值，L 2。

对式(10)两边取对数得t K x f V azln ln )(ln[αα+= (11)为了便于说明计算过程，令])(ln[azx f V Y α= t X ln = （12）则式(11)可转化为X K Y α+=ln (13)于是可利用线性回归法来估算入渗参数K 和α，具体计算步骤如下：(1)线性回归。

对灌水过程中的不同时刻i t （=i 1，2，3…），计算i X ；设入渗指数初始值5.0)0(=α，将)0(α代入到)(αf 中并计算得一系列i Y ，视等式(13)右边的K ln 和α分别为线性回归的截距和斜率，对i Y 和i X 进行线性回归，便可得到入渗参数的第一次回归值，记为)1(K 和)1(α。

(2)迭代。

用)1(α代替)0(α，重复步骤(1)的计算过程，可得到入渗参数的第二次回归值，记为)2(K 和)2(α，如此进行m 次迭代便可得到入渗参数的回归值）（、m m Kα)(。

(3) 误差控制。

对第1−m 次和第m 次线性回归得到的入渗参数）（、1)1(−−m m K α和）（、m m K α)(，利用式(14)计算两次回归求得参数之间的误差。

)1(1−−−=m m m K K K δ，)1()(1−−−=m m m a ααδ (14)式中：1−m Kδ,1−m aδ分别为K 和α的第m 次回归和第1−m 次回归之间的误差。

如果两次回归得到的入渗参数同时满足式(15)，则停止迭代，得到最终参数为)(m K K =，)(m αα=。

K m Kεδ≤−1，ααεδ≤−1m (15)式中：εK 和αε分别为K 和α的允许误差。

(4)精度计算。

精度用残差平方和e S 表示，计算式如下2)(zj zj e V V S −= (16) 式中：zj V 为由试验测得的水流推进到第j 个水深控制点时的单宽入渗水量，L 2；zj V 为利用入渗参数计算得到的第j 个水深控制点时的单宽入渗水量，L 2。

对)0(α在0到1之间任意取值，计算发现，)0(α的不同取值不影响迭代的稳定性，为了加快迭代的收敛速度，建议)0(α在0.3到0.7之间取值。

采用上述线性回归及迭代计算入渗参数，比目前使用广泛的M 法大大减少了计算工作量，而且在计算地表储水量时，避开了难以定值的地表储水形状系数，因而具有更高的计算精度。

2 试验验证2.1 试验区概况及试验方法试验于2005年在中科院南皮生态农业试验站内进行。

试验区位于河北省沧州市南约40km ，地处北纬38。

06′，东经116。

40′。

该地区多年平均降水量550mm ，平均日照时间2318h ，水面蒸发1025~1218mm ，地下水埋深5~7m ，耕地土壤为轻质壤土，灌前测定试验地土层1m 内的土壤含水量15.85％~19.63%，土壤干容重1.42g/cm 3。

在2条畦田（编号分别为B 1和B 2）中分别进行了2次灌水试验，畦田规格及入畦流量见表1。

表1 畦田规格及入畦流量 Table 1 Border data and inflow discharge 畦宽畦长平均纵坡单宽流量畦田编号 /m /m% /m 3/min·mB 1 2.0 95 0.043 0.140 B 22.0 95 0.0640.081为了方便观测地表水流的推进过程和地表水深资料，在畦田两侧从畦首起每隔5m 布置一个标杆，并在畦田中心线上每隔10m 布置一个水位传感器。

灌水过程中用水表和秒表分别记录灌水流量和水流前锋推进时间；用奥得赛电容式水位传感器(Odyssey™4.5，精度0.1mm, 量程0.5m)记录距畦首0m 、10m 、20m 、30m 、40m 、50m 、60m 、70m 和80m 处的地表水深，直至灌水结束试验。

用水来自机井抽取的地下水。

2.2 水流推进过程及地表水深对畦田B 1的地表水流推进资料，用幂函数拟合得到水流推进过程函数关系式如下7721.10151.0st s = 9999.02=R (17)式中：s t 、s 的意义同前。

对畦田B 1的水深数据采用式(2)建立不同时刻的地表水深与水流推进距离的函数关系式，拟合参数及其相关系数见表2。

畦田B 1、B 2地表水深与水流推进距离拟合曲线见图2，图3。