吳彩麗，許 迪，白美健，李益農(1.國家節(jié)水灌溉工程技術研究中心，北京 100048；2.中國水利水電科學研究院 水利研究所，北京 100048)

地面灌溉是目前我國應用最為廣泛的田間灌水方法，但現(xiàn)有地面灌溉田間管理粗放，灌溉過程可控制性差，地面灌溉性能指標值普遍偏低，此問題已成為制約我國地面灌溉技術發(fā)展的瓶頸。加強對灌溉過程的控制和管理，對于提高地面灌溉田間管理水平，改進地面灌水質量具有重要的理論和實際意義。精細地面灌溉智能控制技術的原理是在灌溉過程中通過實時采集田塊前一部分數(shù)據(jù)量的灌溉信息估算土壤入滲參數(shù)和田面糙率系數(shù)，再借助數(shù)值模擬手段對整個灌水過程進行預測，從而優(yōu)化決策出灌溉性能最佳時所需的灌水時間或入畦流量，并根據(jù)優(yōu)化決策方案對灌溉管理參數(shù)(灌水時間或入畦流量)實施控制。由于土壤入滲參數(shù)的空間變異性，若實時采集的灌溉信息數(shù)據(jù)量過少，則估算出的土壤入滲參數(shù)難以代表田塊的整體狀況，會給灌溉模擬和優(yōu)化控制帶來較大誤差；實時采集的灌溉信息數(shù)據(jù)量越多，估算出的土壤入滲參數(shù)越能代表田塊的整體狀況，但若采集的灌溉信息數(shù)據(jù)量過多，則會錯過最佳調節(jié)時機，難以達到最佳灌溉性能。因此，確定需要采集灌溉信息的合理數(shù)據(jù)量是根據(jù)已知灌溉信息準確預測未知灌溉過程，實現(xiàn)對灌溉全過程控制和管理的首要問題。國內外學者對此進行了一些研究，Clemmens等(1992年)通過對灌溉過程中0%、25%、40%、50%、75%和100%的水流推進數(shù)據(jù)量采用貝葉斯方法估算土壤入滲參數(shù)和田面糙率系數(shù)，并將得到的土壤入滲參數(shù)和田面糙率系數(shù)代入地面灌溉模擬模型，對整個水流運動過程進行預測模擬，將水流推進過程模擬值與實測值之間的平均相對誤差進行比較，提出基于水流推進至75%時的數(shù)據(jù)量進行控制相對合理；董孟軍(2010年)對于畦長200 m的田塊選取推進距離60,70,80,…,200 m的15組灌溉信息，采用直接法估算土壤入滲參數(shù)，并將得到的15組土壤入滲參數(shù)代入地面灌溉模擬模型，對整個水流運動過程進行預測模擬，以水流推進過程模擬值與實測值之間的平均相對誤差在15%以內，并且與利用全部水流推進數(shù)據(jù)的平均相對誤差之差在10%以內的最短推進距離為原則，提出基于水流推進至60%時的數(shù)據(jù)量進行控制最為合理。上述學者提出的實施控制需采取的數(shù)據(jù)量均是根據(jù)某一些特定的試驗數(shù)據(jù)得出的經(jīng)驗值，是否具有廣適性不得而知，本文擬基于貝葉斯統(tǒng)計方法，從理論上論證精細地面灌溉智能控制系統(tǒng)需采用的能夠準確估算入滲參數(shù)并預測灌溉全過程的合理數(shù)據(jù)量。

1 方法與試驗

1.1 貝葉斯統(tǒng)計方法

在精細地面灌溉智能控制中，確定不同精度要求所需的合理數(shù)據(jù)量時用到了貝葉斯統(tǒng)計方法?；谪惾~斯統(tǒng)計方法，將先驗信息、樣本信息和似然函數(shù)相結合推導后驗分布及置信區(qū)間，從理論上論證精細地面灌溉智能控制需采用的能夠準確估算入滲參數(shù)并預測灌溉全過程的合理數(shù)據(jù)量。貝葉斯公式包含了歸納推理的一種思想，貝葉斯統(tǒng)計方法的基本觀點是由貝葉斯公式引申而來的。貝葉斯公式是概率論中的一個著名公式，是由英國學者貝葉斯(T.Bayes)在1763年提出來的：設試驗E的樣本空間為S，A為E的事件，B1，B2，…，Bn為S的一個劃分，即某個過程的若干可能的前提，且P(A)>0,P(Bi)>0 (i=1,2,…,n)，則：

(1)

式中：P(Bi)是人們事先對各前提條件出現(xiàn)可能性大小的估計，稱之為先驗概率；P(Bi∣A)既是對以A為前提下Bi的出現(xiàn)概率的重新認識，稱P(Bi∣A)為后驗概率。

1.2 灌溉模擬模型

采用基于混合數(shù)值解法的一維全水動力學畦灌模型(章少輝，2012年)對畦灌水流運動過程進行模擬。所謂混合數(shù)值解法，即基于隱-顯混合時間格式，利用有限差分法、矢通量分裂與通量差分混合格式(AUSM)有限體積法和有限單元法分別對全水動力學畦灌模型中各矢量項進行空間離散，對形成的3對角矩陣型控制方程代數(shù)方程組進行數(shù)值求解。

基于混合數(shù)值解法的一維全水動力學畦灌模型中采用地表水入滲經(jīng)驗公式中應用最為廣泛的是Kostiakov(1932年)入滲公式(z=kta，z為累積入滲水深，mm；k為入滲參數(shù)，mm/sα；a為入滲指數(shù)，主要取決于土壤質地與結構；t為土壤入滲受水時間，s)描述土壤入滲特性。

模型輸入?yún)?shù)為：畦田規(guī)格、田面微地形、單寬流量、關口時間、Kostiakov入滲參數(shù)及糙率系數(shù)。模型輸出參數(shù)為：灌水效率、灌水均勻度及儲水效率等。

一維全水動力學模型方程：

(3)

式中：x為水平向坐標，m；t為水流運動時間，s；h為地表水深，m；g為重力加速度，m/s2；U為沿x水平坐標向的垂向均布流速，m/s；Q為沿x水平坐標向的單寬流量，m3/(s·m)；ic為地表水入滲率，m3/(s· m2)；b為畦面高程，m；n為曼寧糙率系數(shù)，m/s1/3。

1.3 試驗田塊基本信息

1.3.1先驗信息

設先驗信息來自按畦長不同選取的大興(12個田塊冬灌、春灌)、河北鹿泉(10個田塊冬灌)和山東無棣(8個田塊冬灌、春灌)3個試驗區(qū)的50組典型田塊的歷史灌溉資料，典型田塊基本信息見表1。實測灌溉資料主要包括田塊規(guī)格、田面微地形、灌溉開始和關口時間、灌溉水流推進過程、入畦流量等，田面微地形和水流推進消退時間的測點間距為5 m。

表1 先驗信息田塊灌溉基本資料

續(xù)表1 先驗信息田塊灌溉基本資料

1.3.2樣本信息

取樣本信息為來自新疆222兵團典型田塊的20組灌溉資料，X=(X1,X2,…,X20)，典型田塊基本信息見表2。實測灌溉資料也包括田塊規(guī)格、田面微地形、灌溉開始和關口時間、灌溉水流推進過程、入畦流量等，田面微地形和水流推進消退時間的測點間距為5 m。

表2 樣本信息田塊灌溉基本資料

2 結果與分析

2.1 確定先驗分布函數(shù)

因水流推進時間影響相對較大，選定灌溉水流推進過程模擬值與實測值之間的平均相對誤差(AREadv)為隨機變量。貝葉斯方法的一個主要問題是如何利用先驗信息確定先驗分布，因為先驗分布有很大的主觀性、隨意性。采用常用的直方圖法推斷其先驗分布，估計先驗分布密度函數(shù)。

現(xiàn)以選取60%數(shù)據(jù)量為例，確定水流推進過程模擬值與實測值之間的平均相對誤差的先驗分布,并分析根據(jù)先驗灌溉信息準確預測控制后驗灌溉過程所需的合理數(shù)據(jù)量。對于50組先驗數(shù)據(jù)，根據(jù)60%的實測水流推進數(shù)據(jù)采用反求法計算入滲參數(shù)k和a，并將得到的土壤入滲參數(shù)代入地面灌溉模擬模型(田面糙率系數(shù)取經(jīng)驗值：冬灌取0.08，春灌取0.1，且經(jīng)過驗證不同的經(jīng)驗取值對本文結論影響不大)，對整個水流運動過程進行預測模擬，計算整個水流推進過程模擬值與實測值之間的平均相對誤差。

采用基于混合數(shù)值解法的一維全水動力學畦灌模型對畦灌水流運動過程進行模擬，水流推進過程模擬值與實測值之間的平均相對誤差如下(該誤差值為連續(xù)型隨機變量AREadv，則π(AREadv)為AREadv的先驗分布，π(AREadv|x)為AREadv的后驗分布)。

(4)

式中：t0i-adv是水流推進過程中點i處的實測值；tsi-adv是點i處的水流推進模擬值；n′為測點數(shù)。

計算出50組水流推進過程模擬值與實測值之間的平均相對誤差最小值、最大值分別為2.18% 、27.46%，即所有數(shù)據(jù)落在區(qū)間[2.18，27.46]上，現(xiàn)取能覆蓋區(qū)間[2.18，27.46]的區(qū)間[0,28]，將區(qū)間[0,28]等分為7個小區(qū)間，小區(qū)間的長度記為Δ，Δ=(上限-下限)/7=4。數(shù)出落在每個小區(qū)間內的數(shù)據(jù)的頻數(shù)fi，計算出頻率fi/n(n=50,i=1,2,…,7)見表3。

表3 平均相對誤差頻率計算

根據(jù)表3中數(shù)據(jù)自左向右依次在各個小區(qū)間上作以fi/n/Δ為高的小矩形，見直方圖1。小矩形的面積就等于數(shù)據(jù)落在該小區(qū)間的頻率fi/n。從圖1中看出，直方圖的外廓曲線中間高，兩頭低，有一個峰，類似正態(tài)分布曲線。

圖1 先驗分布頻率直方圖

為進一步確認先驗分布是否服從正態(tài)分布，作χ2擬合檢驗如下。

H0:AREadv的概率密度為：

AREadv<+∞

(5)

由極大似然估計法得μ，τ2的估計值分別為14.7和5.02，在H0下將θ可能取值的區(qū)間(-∞,+∞)分為7個小區(qū)間，并取事件Ai(見表3中第1列所示)。若H0為真，則θ的概率密度函數(shù)的估計為：

(6)

按式(5)查正態(tài)分布的分布函數(shù)表可得概率P(Ai)的估計，計算見表4。

表4 檢驗計算

現(xiàn)χ2=52.33-50=2.33, 因為χ20.1(k-r-1)=χ20.1(7-2-1)=7.779>2.33，故在0.1水平下接受H0，由此可從理論上證明水流推進過程模擬值與實測值之間的平均相對誤差的先驗分布服從正態(tài)分布。正態(tài)分布(Normal distribution)又名高斯分布(Gaussian distribution)，是一個在數(shù)學物理及工程領域都非常重要的概率分布，在統(tǒng)計學的許多方面有著重大的影響力。在自然現(xiàn)象和社會現(xiàn)象中，大量隨機變量都服從或近似服從正態(tài)分布。

2.2 根據(jù)先驗地表水流推進數(shù)據(jù)控制畦灌過程的適宜數(shù)據(jù)量

由水流推進過程模擬值與實測值之間的平均相對誤差的先驗分布π(AREadv)～N(14.7,52),結合樣本信息推斷后驗分布π(AREadv|x)及置信區(qū)間。

(8)

由此可以獲得θ的1-α置信區(qū)間。令：

同理，可以根據(jù)灌溉信息采集數(shù)據(jù)量為40%，50%，60%，70%，80%的計算數(shù)據(jù)，分析部分數(shù)據(jù)推求整個水流運動過程的模擬值與實測值的平均相對誤差(AREadv)隨灌溉信息采集數(shù)據(jù)量x的變化規(guī)律，作散點圖并畫出趨勢線見圖2。

圖2 AREadv隨灌溉信息采集數(shù)據(jù)量x的變化趨勢散點與擬合曲線

由圖2可知，AREadv隨灌溉信息采集數(shù)據(jù)量x呈冪函數(shù)變化趨勢，其冪指數(shù)為負數(shù)的減函數(shù)，擬合方程為：

AREadv(x)=0.041 8x-2.50≤AREadv≤1

(9)

由此可得，針對由先驗灌溉信息準確預測控制后驗灌溉過程所需的合理數(shù)據(jù)量研究，根據(jù)貝葉斯統(tǒng)計方法可論證出如下結論：對于一般的田塊(即與本文試驗田塊的灌溉技術要素相近)，當采用灌溉信息采集數(shù)據(jù)量分別為40%、50%、60%、70%和80%的實測水流推進數(shù)據(jù)時，水流推進過程模擬值與實測值之間的平均相對誤差為40%、25%、15%、10%和7%，在精細地面灌溉智能控制技術實際應用中，可以根據(jù)需要的精度要求選擇相應的灌溉信息數(shù)據(jù)量。

前人曾通過實測數(shù)據(jù)對灌水技術要素與灌溉信息數(shù)據(jù)量的影響進行研究，結果表明灌溉信息數(shù)據(jù)量受入畦單寬流量、田面縱坡、畦長等灌溉要素的影響，隨入畦單寬流量、田面縱坡、畦長的增大而減小，隨微地形的增大而增大(董孟軍，2010年)。因此，對于灌水技術要素比較特殊的田塊，需要根據(jù)實際情況選擇適宜的灌溉信息數(shù)據(jù)量。

3 結 論

本文基于貝葉斯統(tǒng)計方法，根據(jù)先驗信息、樣本信息和似然函數(shù)，推導后驗分布及置信區(qū)間，從理論上證明精細地面灌溉智能控制需采用的能夠準確估算入滲參數(shù)并預測灌溉全過程的合理數(shù)據(jù)量。綜合分析得出以下結論：

(1)采用60%的實測水流推進數(shù)據(jù)估算土壤入滲參數(shù)并預測灌溉全過程產(chǎn)生的平均相對誤差落入[12.96,17.38]區(qū)間的概率為95%，即采用60%已經(jīng)數(shù)據(jù)量推斷未知灌溉過程是相對準確與合理的。

(2) 要使部分數(shù)據(jù)推求整個水流運動過程的模擬值與實測值的平均相對誤差控制在40%、25%、15%、10%和7%左右，所需的合理數(shù)據(jù)量為40%、50%、60%、70%和80%。在精細地面灌溉智能控制技術實際應用中，可以根據(jù)需要的精度要求選擇相應的灌溉信息數(shù)據(jù)量。本研究成果為精細地面灌溉智能控制系統(tǒng)中適用條件和合理數(shù)據(jù)量的確定提供了理論依據(jù)。

□

[1] 李益農. 我國地面灌溉技術評價及其發(fā)展趨勢[C]∥ 中國水利學會首屆青年科技論壇論文集. 北京：中國水利水電出版社，2004:201-207.

[2] Clemmens A J, Keats J B. Bayesian inference for feedback controlⅠ:theory[J]. Journal of Irrigation and Drainage Engineering,1992,118(3): 397-415.

[3] Clemmens A J, Keats J B. Bayesian inference for feedback control.Ⅱ:surface irrigation example[J]. Journal of Irrigation and Drainage Engineering,1992,18(3):416-432.

[4] Clemmens A J. Feedback control of basin irrigation system [J], Journal of Irrigation and Drainage Engineering, 1992,118(3):480-496.

[5] 董孟軍. 基于推進距離的一維畦灌反饋控制方法研究[D].北京：中國水利水電科學研究院，2010.

[6] 白美健，李益農，許 迪，等. 精細地面灌溉實時反饋控制技術適應性分析[C]∥中國農業(yè)工程學會2014年學術年會論文集.

[7] 章少輝，許 迪，李益農. 基于混合數(shù)值解法的一維全水動力學畦灌模型[J]. 農業(yè)工程學報，2009,(9)：7-14.

[8] 盛 驟，謝式千，潘承毅. 概率論與數(shù)理統(tǒng)計[M]. 北京：高等教育出版社，2001.

[9] 韋來生，張偉平. 貝葉斯分析[M]. 合肥：中國科學技術大學出版社，2013.

[10] 茆詩松，湯銀才. 貝葉斯統(tǒng)計[M]. 北京：中國統(tǒng)計出版社，2012.

[11] Clemmens A J. Roughness or infiltration from surface irrigation advance [C]∥ ASCE. Proc. Irrigation and Drainage, 1990:46-53.

[12] Walker W R, Busman J D. Real-time estimation of furrow infiltration [J]. Journal of Irrigation and Drainage Engineering, 1990,116(3):299-318.

[13] Burt C M,Clemmens A J, Strelkoff T S, etal. Irrigation performance measure: efficiency and uniformity[J]. Journal of Irrigation and Drainage Engineering, 1997,123(6):423-442.

[14] Kanya L, Khatri, R J Smith. Toward a simple real-time control system for efficient management of furrow irrigation[J]. Journal of Irrigation and Drainage Engineering, 2007,56(4):463-475.

[15] R K Koech, R J Smith, M H Gillies. A real-time optimisation system for automation of furrow irrigation[J]. Journal of Irrigation and Science, 2014,32:319-327.

[16] 許 迪, 李益農. 精細地面灌溉技術體系及其研究的進展[J]. 水利學報, 2007,38(5):529-536.

[17] 白美健. 微地形和入滲空間變異及其對畦灌系統(tǒng)影響的二維模擬評價[D]. 北京：中國水利水電科學研究院，2007.

[18] 許 迪，蔡林根，茆 智，等. 引黃灌區(qū)節(jié)水決策技術應用研究[M]. 北京：中國農業(yè)出版社，2004.

[19] 章少輝. 土壤入滲參數(shù)和糙率系數(shù)的優(yōu)化反求方法及應用[D]. 北京：中國水利水電科學研究院，2005.

[20] 王占軍. 地面灌溉參數(shù)與灌水技術要素優(yōu)化研究[D]. 西安：西安理工大學，2009.