劉 琨，黃冠華

（1.中國農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083；2.中國－以色列國際農(nóng)業(yè)研究培訓(xùn)中心，北京 100083）

1 研究背景

相比于有壓灌溉，溝灌具有低能耗和低投入的優(yōu)點，被廣泛用于寬行作物，例如棉花，玉米和蔬菜。然而，由于溝灌系統(tǒng)的設(shè)計不合理，導(dǎo)致灌溉水利用效率較低，甚至引起次生鹽堿化和地下水污染等問題［1］。通過對溝灌系統(tǒng)的優(yōu)化設(shè)計，可以極大地提高灌水效率［2］。然而，通過田間試驗的方法尋找合理的灌水技術(shù)要素組合費時費力，為此人們提出了很多數(shù)學(xué)模型用于溝灌系統(tǒng)的設(shè)計和管理，其中一些灌溉模型已經(jīng)被開發(fā)成應(yīng)用程序，例如WinSRFR［3］和SIRMOD［4］。數(shù)學(xué)模型模擬可以為用戶采用適宜的灌溉系統(tǒng)設(shè)計方案、評價灌溉系統(tǒng)性能提供有效手段，以達到提高灌溉質(zhì)量的目的［5］。

入滲量估計是溝灌水流模擬研究中的重要問題之一。只有對入滲進行準(zhǔn)確的估計，才能正確地描述灌溉過程中水分的運動，對灌水效果做出合理的評價。現(xiàn)有的溝灌模型大多采用經(jīng)驗公式估計入滲量，經(jīng)驗公式具有計算簡單，參數(shù)較少等優(yōu)點。然而在實際溝灌過程中，溝中水深和濕周不斷變化，影響累積入滲量大?。?］，由于經(jīng)驗公式往往沒考慮土壤的水力特性、水深和濕周變化等因素的影響，導(dǎo)致入滲量的估計出現(xiàn)較大偏差［7］。

隨著對入滲機理的深入認識，人們提出了很多基于物理過程的溝灌入滲模型。例如W?hling等［8］通過一個形狀轉(zhuǎn)換公式將溝灌二維入滲描述成一系列一維入滲的累加，該模型具有較高的精度，但是當(dāng)計算節(jié)點較多時，模型計算效率較低。Warrick等［9］提出一個溝灌近似入滲模型，將溝灌二維入滲看成一維入滲和邊界效應(yīng)之和，并通過數(shù)值試驗驗證了模型，但此模型僅考慮了積水深度恒定的情況，與實際的溝灌過程中積水深度變化的情況有較大差異。Bautista等［10］通過對Warrick模型改進，使其可以應(yīng)用于積水深度變化條件下的入滲計算，該模型只包含一個經(jīng)驗參數(shù)γ，計算較為簡單。Warrick等［9］的研究表明，γ與土壤水力參數(shù)、積水深度和初始有效含水率有關(guān)。Bautista等［10］的研究表明，在定水頭條件下γ值與積水深度有關(guān)，在變水頭條件下γ值可取為常數(shù)。

目前關(guān)于參數(shù)γ的研究還只是定性地分析了其與影響因素之間的關(guān)系，并給出γ的取值范圍，尚未提出參數(shù)γ的定量估計方法。在實際中，土壤的濕潤程度、土壤水力特性和積水深度在灌水過程中隨著時空變化，導(dǎo)致γ的取值具有時空變化特征，若將γ視為常數(shù)，則可能導(dǎo)致入滲計算的不準(zhǔn)確，同時參數(shù)γ時空變化增加了其估計的難度。

近年來，人工智能技術(shù)，如人工神經(jīng)網(wǎng)絡(luò)、遺傳算法和基因表達式編程（GEP）算法等被引入水文與水資源學(xué)領(lǐng)域，并被廣泛應(yīng)用于模型參數(shù)優(yōu)化中。例如，Valipour等［11］使用遺傳算法優(yōu)化修正Kostiakov入滲模型參數(shù)，魯帆等［12］采用多智能體遺傳算法估計了馬斯京根模型參數(shù)，彭昱忠等［13］使用GEP算法獲得了單一重現(xiàn)期暴雨強度計算模型的參數(shù)。GEP算法是一種通用的自適應(yīng)隨機搜索算法，能夠在缺乏先驗知識的情況下，依賴實驗數(shù)據(jù)的挖掘建立較為準(zhǔn)確的預(yù)測預(yù)報公式。與傳統(tǒng)的遺傳算法等相比，GEP算法運算更加靈活，計算效率也更高［14］。

本文采用GEP算法，建立γ與其敏感因子之間的顯式表達式，為溝灌入滲模型提供了有效的參數(shù)估計方法，并通過比較溝灌入滲模型與HYDRUS-2D計算結(jié)果，檢驗采用GEP算法估計γ方法的精度。

2 理論與方法

2.1 溝灌入滲模型 Bautista等［10］提出的溝灌入滲模型表達如下：

式中：I2D為單位溝長的累積入滲量，cm2；I1D為一維累積入滲量，cm；t為時間，min；W為濕周，cm；θS、θ0分別為飽和含水率和初始含水率，cm3/cm3；S為土壤吸濕系數(shù)，cm/min0.5；γ為經(jīng)驗參數(shù)。

因為積水深度沿濕周變化，計算I1D時使用的積水深度應(yīng)該小于實際的積水深度。Bautista等［15］指出可使用濕周平均水深代替實際的積水深度，濕周平均水深可采用下式計算：

式中：hW為濕周平均水深，cm；h為積水深度，cm；z（x）為濕周節(jié)點上的垂向坐標(biāo)，cm。

hW也用于計算土壤吸濕系數(shù)：

式中：KS為飽和導(dǎo)水率，cm/min；hf為土壤濕潤鋒處的壓力水頭，cm，可以通過Neuman公式計算［16］：

式中：h0為初始含水率對應(yīng)的壓力水頭，cm；K（h）為導(dǎo)水率，cm/min。

整理式（1）可得γ的表達式：

2.2 參數(shù)γ估計的GEP算法

2.2.1 數(shù)據(jù)集合 使用GEP算法時，首先需要建立GEP算法的有效數(shù)據(jù)集合，本研究中通過數(shù)值試驗獲得有效數(shù)據(jù)集合。選擇了6種類型土壤進行數(shù)值試驗，代表不同的土壤質(zhì)地。對于每種土壤，模擬了不同溝截面形狀、積水深度和初始含水率條件下的入滲，根據(jù)式（5）計算出不同模擬情景下的γ值。表1總結(jié)了不同模擬情景的土壤水力特性、積水深度和溝截面形狀。溝截面形狀及尺寸見圖1。

式（5）中的二維和一維入滲量分別通過HY?DRUS-2D［17］和 HYDRUS-1D［19］計算得到。圖 2 給出了HYDRUS-2D的計算區(qū)域。模型的上邊界為定水頭邊界，根據(jù)積水深度設(shè)置；下邊界條件為自由排水邊界；其余邊界為零通量邊界；初始條件為初始含水率對應(yīng)的壓力水頭。HYDRUS-1D計算區(qū)域為200 cm深的土壤剖面，上邊界條件為定水頭邊界，根據(jù)濕周平均水深設(shè)置；下邊界條件為自由排水邊界；初始條件為初始含水率對應(yīng)的壓力水頭。將計算的二維和一維入滲量代入式（5），進而獲得不同模擬情景下的γ值。

表1 不同模擬情景的土壤水力特性、積水深度和溝截面形狀

圖1 溝截面形狀及尺寸（單位：cm）

圖2 HYDRUS-2D計算區(qū)域（單位：cm）

表2 數(shù)據(jù)集合統(tǒng)計指標(biāo)

根據(jù)數(shù)值試驗得到的162數(shù)據(jù)點，以飽和導(dǎo)水率Ks、積水深度h、有效含水率Θ（Θ=（θ0-θr）/（θS-θr））、進氣值參數(shù)α、孔徑分布指數(shù)n、溝深D和溝頂寬B作為輸入變量，γ作為輸出變量，建立GEP算法的有效數(shù)據(jù)集合，數(shù)據(jù)集合的統(tǒng)計指標(biāo)見表2。將有效數(shù)據(jù)集合隨機分成訓(xùn)練集合（70%）和驗證集合（30%），訓(xùn)練集合包含112個數(shù)據(jù)點，用于GEP算法的運算；驗證集合包含50個數(shù)據(jù)點，用于檢驗GEP算法的預(yù)測效果。

2.2.2 輸入因子 分析γ的敏感參數(shù)，確定GEP算法的輸入因子。分析的參數(shù)包括KS、Θ、α、n、h、D和B。以壤土和G2溝截面形狀為基準(zhǔn)，初始參數(shù)設(shè)置如下：KS=0.0173cm/min、α=0.036cm-1、n=1.56、Θ=0.49、h=10 cm、D=40 cm、B=20 cm。將所有的參數(shù)以10%的間隔從-50%增加到+50%，根據(jù)γ的變化情況確定GEP算法的輸入因子。

設(shè)置不同的輸入因子組合，比較不同組合GEP算法的預(yù)測效果，確定最優(yōu)組合，建立γ的顯式表達式。使用決定系數(shù)（R2）、均方根誤差（RMSE）、平均絕對誤差（MAE）和平均相對誤差（MRE）作為評價指標(biāo)，其計算方法如下：

式中：Oi和Pi為實測值和預(yù)測值；和為實測和預(yù)測值的平均值；n為比較變量的數(shù)量。

2.2.3 GEP算法過程 這里對參數(shù)估計的GEP算法過程的步驟進行簡要總結(jié)，詳細的介紹請參考Ferreira［11］。算法的主要步驟如下：（1）初始參數(shù)設(shè)置，包括適應(yīng)度函數(shù)選擇、終點集和函數(shù)集選擇、染色體結(jié)構(gòu)選擇、連接函數(shù)選擇、遺傳算子及其相應(yīng)的發(fā)生概率選擇等；（2）根據(jù)遺傳算子對當(dāng)前種群進行遺傳操作；（3）染色體解碼，評價當(dāng)前種群的適應(yīng)度；（4）根據(jù)終止準(zhǔn)則判斷是否迭代，若不滿足終止準(zhǔn)則，生成新的種群，代數(shù)加一；若滿足終止準(zhǔn)則，停止計算。

重復(fù)（2）—（4）步驟達到預(yù)先指定的代數(shù)或找到一個最優(yōu)解，停止計算。本研究中，使用GeneXpro程序進行GEP算法計算，GEP算法的主要參數(shù)設(shè)置見表3。

2.3 變積水深度條件入滲計算 將建立的γ顯式表達式應(yīng)用于溝灌入滲模型，計算變積水深度條件下的累積入滲量并與HYDRUS-2D計算結(jié)果比較，檢驗GEP算法對γ的估計精度。在計算中，溝中間位置截面的積水深度變化采用零慣量模型和Bautista入滲模型耦合的溝灌水流模擬模型得到，作為入滲模型的上邊界條件。模型的輸入?yún)?shù)見表4，結(jié)果如圖3所示。采用MRE和RMSE評價溝灌入滲模型的效果。

圖3 灌水溝中間位置截面的積水深度變化

3 結(jié)果分析與討論

3.1γ的敏感分析 圖4給出了γ對參數(shù)變化的響應(yīng)結(jié)果，可以看出，γ對積水深度最敏感，壤土條件下γ隨著積水深度的增大而減小。需要注意的是，γ與積水深度之間并不存在確定的正相關(guān)關(guān)系，還受土壤類型等因素的影響。Warrick等［20］的研究表明，砂壤土條件下γ隨積水深度增大而增大，粉質(zhì)壤土條件下γ隨積水深度增大而減小。由圖4可知γ對有效含水率和溝深也很敏感，這與Warrick等［9］的研究結(jié)果一致。γ隨著飽和含水率、進氣值參數(shù)和孔徑分布指數(shù)的變化在-10%～10%范圍內(nèi)波動（圖4），表明γ對土壤水力特性參數(shù)也較為敏感，這與Bautista等［15］的研究結(jié)果一致。

表3 GEP算法參數(shù)設(shè)置

表4 模型輸入?yún)?shù)

圖4 γ對各參數(shù)變化的敏感性分析

參數(shù)γ主要用于重力對入滲邊界影響的修正，以及土壤剖面含水量分布和吸濕系數(shù)近似處理可能產(chǎn)生誤差的修正［22］。由于積水深度和溝截面形狀決定入滲邊界的范圍，因此γ對積水深度和溝截面形狀最敏感。土壤有效含水率、導(dǎo)水率、進氣值參數(shù)和孔徑分布指數(shù)影響土壤剖面含水量分布情況和吸濕系數(shù)的大小，因此γ對土壤有效含水率和土壤水力特性參數(shù)也較為敏感。

3.2 最優(yōu)輸入因子組合 根據(jù)γ的敏感性分析結(jié)果選擇h、Θ、D、KS、α和n作為GEP算法的輸入因子，并設(shè)置5組輸入因子組合方案。表5給出了不同方案的GEP算法對γ預(yù)測效果的評價結(jié)果。由表5可知，方案1中僅使用γ的三個最敏感參數(shù)Θ、h和D，γ的預(yù)測精度較低。方案2、3、4在方案1的基礎(chǔ)上分別增加了參數(shù)KS、n和α，與方案1相比，方案2、3、4的γ預(yù)測精度都有所提高，這表明土壤水力特性參數(shù)對γ有較大的影響。Bautista等［15］指出，與積水深度相比，土壤類型對γ的影響更加顯著，因此將土壤水力特性參數(shù)作為輸入因子可以有效地提高GEP算法對γ的預(yù)測精度。與方案2相比，方案3、4的預(yù)測精度更高，訓(xùn)練組的R2值分別為0.825和0.775，這表明與飽和導(dǎo)水率相比，孔徑分布指數(shù)和進氣值參數(shù)對γ的貢獻更大。

由表5可知，方案3在訓(xùn)練組中的預(yù)測效果較好，R2值最高，但是在驗證組預(yù)測精度較低，R2值僅為0.576。方案5在訓(xùn)練和驗證組的預(yù)測效果均較好，MAE值最小，分別為0.053 cm2和0.049 cm2。由于方案5考慮的影響因素最全面，預(yù)測效果較好（見圖5），因此確定方案5為最優(yōu)的輸入因子組合。最優(yōu)輸入因子組合GEP算法得到的γ表達式為：

表5 不同輸入因子組合的GEP算法對γ預(yù)測效果評價指標(biāo)

圖5 方案5實測值和GEP算法預(yù)測值之間的回歸分析

3.3 變積水深度條件下γ的估計值 圖6給出了變積水深度條件下不同土壤類型采用GEP算法估計的γ值?？梢钥闯靓秒S著積水深度增加，最終穩(wěn)定在一個恒定值附近。這是由于在一定時間后積水深度變化不大（圖3），此時的入滲過程可視為定積水深度條件下的入滲，因此γ值穩(wěn)定在恒定值附近。對于不同類型的土壤，γ隨著土壤導(dǎo)水率增加而增加，黏壤土的γ值為0.82，砂土的γ值為1.13。這與Bautista等的結(jié)果一致，其研究結(jié)果表明質(zhì)地較粗的土壤的γ值更大。Bautista等給出的黏壤土、砂黏壤土和壤砂土的γ值分別為0.8、0.9和1.0［10］，本研究中估計的3種土的γ穩(wěn)定值分別為0.82、0.88和1.12，與Bautista等的結(jié)果相近，表明GEP算法對γ的估計效果較好。

3.4 變積水深度條件下累積入滲量 將式（10）應(yīng)用于溝灌入滲模型計算累積入滲量，并與HYDRUS-2D計算的“精確”結(jié)果進行比較，結(jié)果如圖7所示。從圖中可以看出，采用GEP算法估計γ值的溝灌入滲模型計算精度較高，MRE值均小于5%，滿足入滲計算的精度要求，這表明所建立的基因表達式（式（10））在變積水深條件下對γ的估計效果較好。

為了進一步研究變積水深度對γ的影響，通過擬合溝灌入滲模型計算的累積入滲量與HYDRUS-2D計算的“精確”結(jié)果，反求出變積水深度的γ值，并與定積水深度的γ值進行比較（見表6），結(jié)果表明變積水深度的γ值與定積水深度的γ值取值接近，這與Bautista等的結(jié)果一致，其研究發(fā)現(xiàn)變積水深度和定積水深度的γ值相近［10］，這表明γ對積水深度變化條件不敏感。雖然基因表達式（式（10））是基于定積水深度條件建立的，在變積水深度條件下對γ的估計效果也較好。GEP算法估計的γ穩(wěn)定值與反求的變積水深度的γ值接近（見表6），因此使用估計的γ穩(wěn)定值計算溝灌累積入滲量，在滿足計算精度要求的同時，可以使溝灌入滲模型的應(yīng)用更加方便。

圖6 變積水深度條件下不同類型土壤采用GEP算法估計的γ值

表6 變積水深度和定積水深度條件的γ值

3.5 土壤類型對邊界效應(yīng)的影響 為了研究土壤類型對邊界效應(yīng)的影響，本文模擬了不同土壤在定積水深度（10 cm）和相同初始土壤含水率（0.25）條件下的入滲。入滲過程中與一維累積入滲量增加速率相比，總?cè)霛B量增加速率更大，因此入滲過程中邊界效應(yīng)對總?cè)霛B量的貢獻逐漸增大（見圖8）。這是在入滲邊界中考慮了水流重力作用影響的結(jié)果［21］，在入滲開始階段，土壤基質(zhì)勢起主導(dǎo)作用，隨著入滲時間的推進，土壤含水率增加，毛管基質(zhì)勢作用逐漸減小，重力勢的作用逐漸增大，因此邊界效應(yīng)對總?cè)霛B量貢獻逐漸增大。由圖9可知對于質(zhì)地較粗的土壤，邊界效應(yīng)對總?cè)霛B量貢獻更大，黏壤土中邊界效應(yīng)的最終貢獻約為25%，砂土中約為37%。這是由于粗質(zhì)土壤的砂粒含量較多，對水分的吸附能力小，重力對水分運動的影響大，進而導(dǎo)致邊界效應(yīng)對總?cè)霛B量的貢獻較大。

圖7 采用GEP算法估計γ值的溝灌入滲模型計算的入滲量與HYDRUS-2D計算結(jié)果比較

圖8 砂壤土總?cè)霛B量和一維累積入滲量模擬結(jié)果

圖9 入滲過程中邊界效應(yīng)對總?cè)霛B量貢獻

4 結(jié)論

本文采用基因表達式編程算法，建立了溝灌入滲模型參數(shù)γ與最優(yōu)組合因子之間的基因表達式，提出了參數(shù)γ的估計方法，并分析了變積水深度對γ值和累積入滲量影響以及土壤類型對邊界效應(yīng)的影響。主要結(jié)論如下：

（1）γ與積水深度、溝形狀和土壤水力特性等因素有關(guān)，其中γ對積水深度、初始有效含水率和溝深度最敏感，建立的基因表達式為γ與積水深度、初始有效含水率、溝深、飽和導(dǎo)水率和進氣值參數(shù)之間的定量關(guān)系。

（2）變積水深度的γ值與定積水深度的γ值取值接近，基于定積水深度條件建立的基因表達式在變積水深度條件下對γ的估計效果較好，與利用HYDRUS-2D模型計算得到的“精確”入滲量相比，應(yīng)用基于γ估計值的溝灌入滲模型計算的累積入滲量誤差小于5%，滿足計算精度要求。對于質(zhì)地較粗的土壤，γ值更大，黏壤土γ值約為0.8，砂土的γ值約為1.2。

（3）入滲過程中受重力作用的影響，邊界效應(yīng)對總?cè)霛B量的貢獻逐漸增大。與細質(zhì)土壤相比，粗質(zhì)土壤中邊界效應(yīng)對總?cè)霛B量的貢獻更大，本研究條件下，黏壤土中邊界效應(yīng)最終貢獻約為25%，砂土中約為37%。

參 考 文 獻：

［1］ ESFANDIARI M，MAHESHWARI B L.Field evaluation of furrow irrigation models［J］.Journal of Agricultural Engineering Research，2001，79（4）：459-479.

［2］ WALKER W R，SKOGERBOE G V.Surface Irrigation：Theory and Practice［M］.Prentice-Hall Inc，Englewood Cliffs，N.J.，1987.

［3］ BAUTISTA E，CLEMMENS A J，STRELKOFF T S，et al.Modern analysis of surface irrigation systems with WINSRFR［J］.Agricultural Water Management，2009，96：1146-1154.

［4］ WALKER W R.SIRMOD III—Surface Irrigation Simulation，Evaluation and Design：User’s Guide and Techni?cal Documentation［M］.Dept.of Biological and Irrigation Engineering，Utah State Univ.，Logan，Utah.2003.

［5］ 許迪，李益農(nóng).精細地面灌溉技術(shù)體系及其研究的進展［J］.水利學(xué)報，2007，38（5）：529-537.

［6］ 劉亶仁，路京選.溝灌二維入滲條件下累計入滲量變化規(guī)律的研究［J］.水利學(xué)報，1989（4）：11-21.

［7］ ZERIHUN D，F(xiàn)URMAN A，WARRICK A W，et al.Coupled surface-subsurface flow model for improved basin irrigation management［J］.Journal of Irrigation and Drainage Engineering，2005，131（2）：111-128.

［8］ W?HLING T H，SCHMITZ G H，MAIHOL J C.Modeling two-dimensional infiltration from irrigation furrows［J］.Journal of Irrigation and Drainage Engineering，2004，130（4）：296-303.

［9］ WARRICK A W，LAZAROVITCH N，F(xiàn)URMAN A，et al.Explicit infiltration function for furrows［J］.Journal of Irrigation and Drainage Engineering，2007，133（4）：307-313.

［10］ BAUTISTA E，WARRICK A W，SCHLEGEL J L，et al.Approximate furrow infiltration model for time-variable ponding depth［J］.Journal of Irrigation and Drainage Engineering，2016，04016045.

［11］ VALIPOUR M，MONTAZAR A A.Optimize of all effective infiltration parameters in furrow irrigation using visual basic and genetic algorithm programming［J］.Australian Journal of Basic and Applied Sciences，2012，6（6）：132-137.

［12］ 魯帆，蔣云鐘，王浩，等.多智能體遺傳算法用于馬斯京根模型參數(shù)估計［J］.水利學(xué)報，2007，38（3）：289-294.

［13］ 彭昱忠，元昌安，林開平，等.暴雨強度計算模型參數(shù)擬合優(yōu)化的新進化方法［J］.廣西大學(xué)學(xué)報：自然科學(xué)版，2013，38（5）：1173-1178.

［14］ FERREIRA C.Gene expression programming：a new adaptive algorithm for solving problems［J］.Complex Sys?tems，2001，13（2）：87-129.

［15］ BAUTISTA E，WARRICK A W，STRELKOFF T S.New results for an approximate method for calculating two-di?mensional furrow infiltration［J］.Journal of Irrigation and Drainage Engineering，2014，140（10）：349-356.

［16］NEUMAN S P.Wetting front pressure head in the infiltration model of Green and Ampt［J］.Water Resources Re?search，1976，12（3）：564-566.

［17］ ?IM?NEK J，?EJNA M，van GENUCHTEN M Th.The Hydrus-2D software package for simulating water flow and solute transport in two-dimensional variably saturated media，Version 1.0，IGWMC-TPS-53［Z］.Interna?tional Ground Water Modeling Center，Colorado School of Mines，Golden，Colo.1996.

［18］ HILLS R G，PORRO I，HUDSON D B，et al.Modeling one-dimensional infiltration into very dry soils.1.Model development and evaluation［J］.Water Resources Research，1989，25（6）：1259-1269.

［19］ ?IM?NEK J，?EJNA M，van GENUCHTEN M Th.HYDRUS 1D software package for simulating the one-dimen?sional movement of water heat and multiple solutes in variably saturated media，version 2.0［Z］.International Ground Water Modeling Center，Colorado School of Mines，Golden，Colo.1998.

［20］ WARRICK A W，LAZAROVITCH N.Infiltration from a strip source［J］.Water Resources Research，2007，43，W03420.

［21］ HAVERKAMP R，ROSS P J，SMETTEM K R J，et al.Three-dimensional analysis of infiltration from the disc in?filtrometer 2.Physically based infiltration equation［J］.Water Resources Research，1994，30（11）：2931-2935.