沈洪政，王仰仁，韓娜娜

(天津農(nóng)學(xué)院水利工程學(xué)院，天津 300384)

0 引 言

參數(shù)擬合初始值的確定是非線性模型參數(shù)優(yōu)化的基本問題。對(duì)于非線性模型參數(shù)擬合人們提出了多種方法[1-3]。李春友[4]等提出用優(yōu)化方法計(jì)算Van Genuchten模型參數(shù)，擬合初始值參考有關(guān)文獻(xiàn)選取。由于所選初值有時(shí)對(duì)擬合值有較大的影響,因此在這種情況下就需要多次調(diào)初值才可獲得較為滿意的擬合效果。錢天偉[5]等用單純形調(diào)優(yōu)法和單純形加速法作為優(yōu)化手段，結(jié)合計(jì)算機(jī)程序進(jìn)行參數(shù)優(yōu)化求解，但其初始值是任意給定的。以上均是從參數(shù)優(yōu)化求解方法角度進(jìn)行研究，對(duì)于優(yōu)化初始值如何選取，沒有給出普遍適用的方法。

在求解非線性模型問題時(shí)，通常利用最小二乘法進(jìn)行參數(shù)優(yōu)化擬合。對(duì)于有明確物理意義的參數(shù)，可以根據(jù)其物理意義進(jìn)行初始值的確定。對(duì)于沒有明確物理意義的模型參數(shù)，通常采用均勻網(wǎng)格法[6]和外推內(nèi)插法[7]等進(jìn)行初始值的選取。對(duì)于以上方法，由于需要進(jìn)行反復(fù)尋優(yōu)，選取的初始點(diǎn)較多，過(guò)程復(fù)雜，使用不方便。因此，本文以擬合精度較高的非線性土壤入滲Green Ampt模型和Horton公式為例，提出采用泰勒級(jí)數(shù)法確定初始值，避免了初始值選取的煩瑣過(guò)程，提高了參數(shù)擬合的效率和精度。

1 方法介紹

1.1 目標(biāo)函數(shù)

(1)

1.2 入滲模型介紹

1.2.1 Green-Ampt模型

Green-Ampt模型是Green W H和Ampt G A根據(jù)毛管理論提出的近似積水模型[8,9]。在土層均勻、初始土壤干燥的薄層積水入滲條件下，利用達(dá)西定律，經(jīng)過(guò)化簡(jiǎn)與推導(dǎo)，可得到Green-Ampt入滲公式[10]：

(3)

式中：i(t)為入滲率，mm/min；ic為穩(wěn)定入滲率，mm/min；I(t)為累積入滲量，mm；t為入滲時(shí)間，min；b為參數(shù)，其中ic和b為待定參數(shù)。

1.2.2 霍頓(Horton)公式

霍頓(1940年)提出的入滲公式為：

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

(4)

(5)

式中：i0為初始入滲率，cm/min；β為衰減參數(shù)，其余符號(hào)意義同前，其中i0、ic、β為待定參數(shù)。

Green-Ampt模型和Horton公式具有較高的擬合精度，因而被得到廣泛的使用。由公式的形式易知，兩者均為非線性模型，其參數(shù)的擬合屬于非線性擬合問題。在利用最小二乘法(OLS)進(jìn)行參數(shù)優(yōu)化擬合時(shí)，兩模型中的參數(shù)i0和ic的初始值可根據(jù)其物理意義進(jìn)行確定，而參數(shù)b及β由于沒有明確的物理意義，用外推內(nèi)插法等方法確定時(shí)過(guò)程煩瑣。針對(duì)此問題，本研究提出利用泰勒級(jí)數(shù)法對(duì)公式進(jìn)行展開，結(jié)合田間實(shí)測(cè)入滲數(shù)據(jù)，獲得未知參數(shù)b和β的初始值，再利用OLS進(jìn)行擬合，即可得到最優(yōu)入滲模型參數(shù)。

1.3 初始值的確定

1.3.1 參數(shù)i0和ic初始值的確定

參數(shù)i0和ic分別表示初始入滲率和穩(wěn)滲率。初始入滲率可選取入滲試驗(yàn)測(cè)定的最大入滲率值作為參數(shù)i0的優(yōu)化初始值，穩(wěn)滲率可選取入滲試驗(yàn)測(cè)定的最小入滲率值作為參數(shù)ic的優(yōu)化初始值。

1.3.2 參數(shù)b初始值的確定

參數(shù)b初始值可采用泰勒級(jí)數(shù)法確定。為了便于泰勒公式的展開，對(duì)于公式(3)將其改寫為以t為因變量的顯式函數(shù)：

(6)

(7)

取式(7)中前3項(xiàng)，代入式(6)中得：

(8)

1.3.3 參數(shù)β初始值的確定

(9)

取式(9)中前3項(xiàng)，代入式(5)中得：

(10)

由式(8)和式(10)知，通過(guò)泰勒級(jí)數(shù)法展開后，非線性模型轉(zhuǎn)化為線性模型，結(jié)合i0和ic的實(shí)際物理意義及田間實(shí)測(cè)入滲數(shù)據(jù)即可確定入滲模型參數(shù)的初始值，進(jìn)而可以利用最小二乘法進(jìn)行參數(shù)擬合?，F(xiàn)以天津農(nóng)學(xué)院農(nóng)田水循環(huán)試驗(yàn)基地測(cè)得的田間入滲數(shù)據(jù)為例，驗(yàn)證方法的合理性。

2 實(shí)例計(jì)算

2.1 試驗(yàn)區(qū)概況

天津農(nóng)學(xué)院農(nóng)田水循環(huán)試驗(yàn)基地位于天津市楊柳青鎮(zhèn)大柳灘村(116°57′E，39°08′N，海拔5.49 m)，全年平均氣溫11.6 ℃。試驗(yàn)田總面積1 hm2，地下水位變幅在2.06～3.70 m之間，試驗(yàn)區(qū)土壤理化性質(zhì)見表1，土壤水分特征曲線參數(shù)[11]見表2和圖1。

2.2 測(cè)試方法

試驗(yàn)在2011年分2次測(cè)試，分別在灌溉冬水的麥田、不灌冬水的麥田和棉花茬地進(jìn)行(試驗(yàn)區(qū)簡(jiǎn)圖見圖2)。測(cè)試時(shí)間為1月6日(氣溫為-5.6 ℃，相應(yīng)麥田測(cè)試點(diǎn)次編號(hào)分別記為1和2)和3月28日(氣溫為11.2 ℃，麥田測(cè)試點(diǎn)次編號(hào)分別記為4和5)；棉花試驗(yàn)田只在3月28日測(cè)試，其測(cè)試點(diǎn)次編號(hào)記為3。

表1 試驗(yàn)區(qū)土壤理化性質(zhì)

表2 試驗(yàn)區(qū)土壤水分特征曲線參數(shù)

圖1 試驗(yàn)區(qū)土壤水分特征曲線圖

試驗(yàn)采用雙環(huán)入滲儀測(cè)定土壤水分入滲過(guò)程[12]。內(nèi)環(huán)橫截面積1 000 cm2，內(nèi)環(huán)直徑35.7 cm，外環(huán)直徑60 cm，內(nèi)外環(huán)高度均為30 cm。試驗(yàn)前，將滲透環(huán)埋于試驗(yàn)田，埋入土壤的深度為10 cm。內(nèi)環(huán)加水保持5 cm水頭，采用手動(dòng)方式向環(huán)內(nèi)加水，利用量筒計(jì)量，使得內(nèi)外環(huán)水位齊平。試驗(yàn)開始后記錄每次灌入水量和相應(yīng)的時(shí)間，其內(nèi)環(huán)灌入水量即為滲入水總量，測(cè)試直到入滲率基本穩(wěn)定為止，記錄內(nèi)環(huán)水量全部滲完為止的時(shí)間。

圖2 試驗(yàn)區(qū)簡(jiǎn)圖(單位：m)

2.3 結(jié)果分析

2.3.1 b值與β值計(jì)算結(jié)果分析

根據(jù)初始值確定方法的介紹，以天津農(nóng)學(xué)院水循環(huán)試驗(yàn)基地在2011年1月6日和2011年3月28日測(cè)得的入滲數(shù)據(jù)為基礎(chǔ)，經(jīng)過(guò)計(jì)算，給出了各測(cè)試點(diǎn)b和β平均值的計(jì)算結(jié)果(表3)以及各測(cè)試點(diǎn)b和β值隨時(shí)間變化的曲線關(guān)系圖(圖3、圖4)。

表3 各測(cè)試點(diǎn)b與β平均值計(jì)算結(jié)果表

由表3知，通過(guò)平均值法計(jì)算得到的Green-Ampt模型中參數(shù)b的值在各測(cè)試點(diǎn)之間變化相對(duì)較大，Horton公式中參數(shù)β的值在個(gè)測(cè)試點(diǎn)之間變化較小。從Green-Ampt模型和Horton公式的實(shí)際意義分析來(lái)看，當(dāng)入滲時(shí)間t→∞時(shí)，入滲速率i趨于穩(wěn)定入滲率ic，當(dāng)入滲時(shí)間t不斷增加時(shí)，入滲速率i也趨于穩(wěn)定入滲率ic。為滿足此變化規(guī)律，要求b值和β值均為正值，這與計(jì)算的結(jié)果相符，表明用取平均值的方法計(jì)算b值和β值是可行的。

圖3 各測(cè)試點(diǎn)b值隨時(shí)間變化曲線

圖4 各測(cè)試點(diǎn)β值隨時(shí)間變化曲線

由圖3和圖4知，在入滲試驗(yàn)初期，由于初始表層土壤含水率及土壤吸力大小等變化較大，入滲速率變化較大，導(dǎo)致b值和β值變化較大。隨著入滲試驗(yàn)的進(jìn)行，入滲速率逐漸穩(wěn)定，各測(cè)試點(diǎn)b值和β值變化也趨于平穩(wěn)。β值基本上在橫軸附近上下變化，b值基本上在79上下變化。說(shuō)明，采用泰勒級(jí)數(shù)法確定的b值和β值具有較好的穩(wěn)定性，且β值較b值更穩(wěn)定。

2.3.2 Green-Ampt模型及Horton公式參數(shù)擬合結(jié)果

以b和β求出的值作為初始值，以累積入滲量I(t)的實(shí)測(cè)值與模擬值誤差平方和最小為目標(biāo)，利用Excel規(guī)劃求解工具，分別對(duì)Green-Ampt模型及Horton公式的參數(shù)進(jìn)行擬合。

表4 Horton公式參數(shù)擬合結(jié)果表

表5 Green-Ampt公式參數(shù)擬合結(jié)果表

由表4和表5知：

(1)采用泰勒級(jí)數(shù)法確定初始值來(lái)描述入滲變化具有較高精度，相關(guān)系數(shù)R2均達(dá)到0.85以上，最大值達(dá)到了0.998 6，具有較好的相關(guān)性。但由于初始入滲率i0和穩(wěn)滲率ic只是簡(jiǎn)單選取測(cè)試數(shù)據(jù)的最大值和最小值，同時(shí)β和b值均采用測(cè)試數(shù)據(jù)的平均值計(jì)算得出，因此，需對(duì)參數(shù)初始值進(jìn)行擬合。由擬合后的結(jié)果知，相關(guān)系數(shù)R2均達(dá)到了0.99以上，擬合精度得到顯著提高。

(2)對(duì)比Horton公式和Green-Ampt模型參數(shù)擬合結(jié)果知，Horton公式擬合參數(shù)均為正值，Green-Ampt模型擬合結(jié)果中穩(wěn)定入滲率 出現(xiàn)了負(fù)值，與穩(wěn)定入滲率 的實(shí)際物理意義不符。從均方差計(jì)算結(jié)果來(lái)看，Horton公式計(jì)算的各測(cè)試點(diǎn)之間均方差的值相差不大，均在1.63上下變化。Green-Ampt模型計(jì)算的各測(cè)試點(diǎn)之間均方差的值相差較大，最大值達(dá)到了7.31，最小值為0.76，變化較大。因此，Horton公式的參數(shù)比Green-Ampt模型的參數(shù)更加穩(wěn)定，公式更具有適用性。

圖5 Green-Ampt模型各測(cè)試點(diǎn)累積入滲量實(shí)測(cè)值與模擬值隨時(shí)間變化曲線

圖6 Horton公式各測(cè)試點(diǎn)累積入滲量實(shí)測(cè)值與模擬值隨時(shí)間變化曲線

圖5和圖6給出了各測(cè)試點(diǎn)Green-Ampt模型和Horton公式累積入滲量模擬值與對(duì)應(yīng)實(shí)測(cè)值隨時(shí)間變化的曲線。由圖可見，各測(cè)試點(diǎn)Green-Ampt模型和Horton公式累積入滲量的模擬值擬合精度均較高。說(shuō)明，采用泰勒級(jí)數(shù)法確定的初始值進(jìn)行參數(shù)擬合具有較高的準(zhǔn)確性。另一方面，由于測(cè)試時(shí)溫度的不同，灌冬水的麥田測(cè)試點(diǎn)1和測(cè)試點(diǎn)5，其累積入滲量之間變化較大，不灌冬水的麥田測(cè)試點(diǎn)2和測(cè)試點(diǎn)4也有相同變化規(guī)律。而對(duì)于測(cè)試溫度相同的麥田測(cè)試點(diǎn)1、2(-5.6 ℃)和測(cè)試點(diǎn)4、5(11.2 ℃)而言，其累積入滲量實(shí)測(cè)值與模擬值隨時(shí)間變化趨勢(shì)相近，表明，溫度對(duì)土壤水分入滲過(guò)程有很大影響。

2.3.3 溫度對(duì)入滲結(jié)果的影響

通過(guò)圖5和圖6的分析結(jié)果可知，溫度對(duì)入滲結(jié)果有很大影響，溫度越高，入滲速率越大，反之越小。茲以參數(shù)穩(wěn)定性較好的Horton公式為例，考慮溫度對(duì)入滲結(jié)果的影響，在式(4)和式(5)中引入溫度修正系數(shù)，

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

(12)

其中：

KT=ea(T-15)

(13)

式中：KT為溫度修正系數(shù)；T為氣溫，℃；a為參數(shù)，15為參考溫度，℃。對(duì)式(11)，將各測(cè)試點(diǎn)(1和5、2和4)在同一時(shí)間的實(shí)測(cè)入滲量數(shù)值之比與不同溫度(-5.6 ℃，11.2 ℃)下的溫度修正系數(shù)KT值之比建立等式關(guān)系[式(14)]，由此便可求得參數(shù)a的值(見表6)，由測(cè)試點(diǎn)1和5求得的a值為0.044 1，由測(cè)試點(diǎn)2和4求得a的值為0.050 7。參數(shù)a的值確定以后，再將不同溫度(-5.6 ℃，11.2 ℃)值代入式(11)中，即可求得溫度修正系數(shù)值。

(14)

式中：T1、T2分別為不同灌水地塊的氣溫；I1(t)、I2(t)分別為不同測(cè)試點(diǎn)的實(shí)測(cè)入滲量值。圖7和圖8分別表示不考慮溫度修正的入滲曲線及考慮溫度修正的入滲曲線。由圖7和圖8知，引入溫度修正系數(shù)以后，不同溫度的入滲曲線很好的重合到一起(參考溫度條件下)，顯著地改善了入滲曲線擬合精度。

表6 a值計(jì)算結(jié)果表

圖7 不考慮溫度修正的入滲曲線

圖8 考慮溫度修正的入滲曲線

3 結(jié) 語(yǔ)

本文以非線性土壤入滲Green-Ampt模型和Horton公式為例，提出采用泰勒級(jí)數(shù)法確定其模型參數(shù)初始值，同時(shí)考慮了溫度對(duì)入滲結(jié)果的影響。結(jié)果表明，在初始值擬合基礎(chǔ)上，相關(guān)系數(shù)均達(dá)到0.85以上，經(jīng)最小二乘法擬合后，實(shí)測(cè)值與模擬值相關(guān)系數(shù)提高到0.99以上?？紤]溫度修正系數(shù)以后，不同溫度的入滲曲線很好的重合到一起(參考溫度條件下)，顯著地改善了入滲曲線擬合精度。

泰勒級(jí)數(shù)法具有更好的理論基礎(chǔ)，避免了非線性規(guī)劃求參中需要設(shè)置多個(gè)初始值的不足，可顯著提高計(jì)算效率，而且可以提高參數(shù)尋優(yōu)的準(zhǔn)確度，適用于無(wú)法線性化的非線性模型參數(shù)初始值的確定。但是，在擬合計(jì)算過(guò)程中Green-Ampt模型參數(shù) 值出現(xiàn)了負(fù)值情況，與其實(shí)際物理意義相矛盾，尚需進(jìn)一步研究。

□

[1] 朱良君，張光輝，任宗萍. 4種土壤入滲測(cè)定方法的比較[J]. 水土保持學(xué)報(bào)，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] 張新民，王根緒，胡想全，等. 優(yōu)化技術(shù)在土壤入滲參數(shù)計(jì)算中的應(yīng)用[J]. 灌溉排水學(xué)報(bào)，2004，(6): 52-54.

[4] 李春友，任 理，李保國(guó). 利用優(yōu)化方法求算Van Genuchten方程參數(shù)[J]. 水科學(xué)進(jìn)展，2001，(12): 473-476.

[5] 錢天偉，陳繁榮，杜曉麗，等. 一種推求van Genuchten方程參數(shù)的高性能優(yōu)化方法[J]. 土壤學(xué)報(bào)，2004， (11): 973-975.

[6] 尚松浩. 水資源系統(tǒng)分析方法及應(yīng)用[M]. 北京：清華大學(xué)出版社，2012:179-180.

[7] 張 瑩. 運(yùn)籌學(xué)基礎(chǔ)[M]. 北京：清華大學(xué)出版社，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入滲模型及其應(yīng)用[J]. 西北農(nóng)林科技大學(xué)學(xué)報(bào)，2007，(2):225-230.

[10] 雷志棟，楊詩(shī)秀，謝森傳. 土壤水動(dòng)力學(xué)[M]. 北京：清華大學(xué)出版社，1989:130-131.

[11] 車 政，王仰仁，王永紅，等. 農(nóng)田土壤水分特征曲線參數(shù)擬合及其剖面變異特性研究[J]. 灌溉排水學(xué)報(bào)，2016:1-5.

[12] 周鑫洋，王仰仁，劉群昌. 四個(gè)入滲公式對(duì)農(nóng)田潮土的適應(yīng)性研究[J]. 灌溉排水學(xué)報(bào)，2015，(8):96-100.