基于貝葉斯統計方法的精細地面灌溉智能控制系統合理數據量研究

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

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

1 方法與試驗

1.1 貝葉斯統計方法

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

(1)

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

1.2 灌溉模擬模型

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

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

模型輸入參數為：畦田規格、田面微地形、單寬流量、關口時間、Kostiakov入滲參數及糙率系數。模型輸出參數為：灌水效率、灌水均勻度及儲水效率等。

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

(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為曼寧糙率系數，m/s1/3。

1.3 試驗田塊基本信息

1.3.1先驗信息

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

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

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

1.3.2樣本信息

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

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

2 結果與分析

2.1 確定先驗分布函數

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

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

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

(4)

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

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

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

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

圖1 先驗分布頻率直方圖

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

H0:AREadv的概率密度為：

AREadv<+∞

(5)

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

(6)

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

表4 檢驗計算

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

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

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

(8)

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

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

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

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

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

(9)

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

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

3 結 論

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

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

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

□

[1] 李益農. 我國地面灌溉技術評價及其發展趨勢[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]∥中國農業工程學會2014年學術年會論文集.

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

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

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

[10] 茆詩松，湯銀才. 貝葉斯統計[M]. 北京：中國統計出版社，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] 白美健. 微地形和入滲空間變異及其對畦灌系統影響的二維模擬評價[D]. 北京：中國水利水電科學研究院，2007.

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

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

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