基于改進的Kostiakov模型的地面灌溉入滲參數與田間糙率系數優化識別研究

孫 玥

(遼寧河庫管服務中心(遼寧省水文局)，遼寧 沈陽 110003)

農田灌溉過程中管理制度及灌水質量評價的重要指標為土壤入滲系數及田面糙率系數，農田灌溉設計基礎數據也包括土壤入滲系數及田面糙率系數[1]。農業是我國國民經濟的基礎，地面灌溉方式是國內90%以上農田采用的灌溉方式[2]。地面灌溉效果改善需要對土壤入滲系數及田面糙率系數進行準確評估，眾多研究成果表明[3- 5]，土壤入滲特性是地面灌溉水流推進的重要影響因素。地面灌溉主要研究內容是建立合適的土壤入滲模型以及對入滲參數進行準確確定[6]。土壤入滲系數確定主要有2種方法，第一種方法是通過田間試驗測定，第二種方法是地面灌溉試驗測定數據進行間接反推計算[7]。田間試驗測定方法由于區域土壤特性有所差異，使得其很難具有代表性。間接反推算法在當前國內應用較為成熟，通過參數優化識別，大大提高土壤入滲參數確定的精度，其中Kostiakov公式當前在國內土壤入滲參數間接反推中應用較為成熟[8- 9]，但傳統Kostiakov公式由于不能考慮時間變量對土壤入滲特性的影響，在土壤入滲參數確定存在一定局限，為此本文對該公式引入時間變量進行改進，建立改進的Kostiakov公式對土壤入滲參數進行確定，結合均勻設計試驗確定地面灌溉水流推進數據，選用地面灌溉水流模擬WinSRFR模型[10]并結合多元回歸方法，建立土壤入滲參數和田面糙率系數優化模型，結合粒子群算法對模型進行求解確定土壤入滲系數及田面糙率系數。此方法土壤入滲系數及田面糙率系數的確定僅需地面灌溉推進水流數據即可。研究成果對于灌區設計中灌溉入滲參數及田面糙率系數確定方法提供借鑒和參考。

1 研究方案

土壤入滲參數和田面糙率系數的取值范圍，結合田間試驗數據進行確定，試驗方案采用均勻試驗進行設計，不同參數組合條件下的地面灌溉水流推進采用WinSRFR模型進行模擬，與實測水流觀測數據進行比較，將地面灌溉水流模擬值和實測值之間差值作為目標函數進行模擬試驗，模擬值和實測值之間差值的平方和與參數之間通過多元回歸方法建立土壤入滲參數和田面糙率系數的參數優化設別模型，確定地面灌溉最優參數組合。

2 地面灌溉水流模型構建

2.1 模型選取

(1)WinSRFR模型

不同參數組合下的地面灌溉水流模型選用WinSRFR模型，模型可以對地面灌溉水流推進過程進行動態模擬，模型計算方程為：

(1)

(2)

其中：

(3)

式中，Q—任意時刻t進入農田斷面流量，L/s;I—單位長度內(x，t)水量下滲量，mm/m;T—入滲時間凈值，min;S0—灌溉田地的坡度，%;St—坡降阻力，%。q—單寬流量，L/s;n—曼寧糙率;y—任意時刻t進行田面的水深，mm。

(2)改進的Kostiakov模型

針對傳統Kostiakov模型不能考慮土壤入滲量隨時間變化的而局限，引入時間變量t對傳統模型進行改進，建立改進的Kostiakov模型對土壤入滲量進行模擬，改進方程為：

Z=Kta+f0t

(4)

式中，Z—水深入滲量，mm;t—地面灌溉水流時間，h;K、a、n—土壤入滲參數。

因此，在進行均勻方案設計時需要考慮的參數分別為K、a、n。

2.2 模型參數確定

地面灌溉水流推進及消退曲線采用WinSRFR模型進行模擬時，需要對以下參數進行確定，其中包括固定和可變2類參數。

(1)固定參數。包括試驗田的尺寸、地形、灌水流量、灌水時間作為輸入固定參數變量，其中沿著試驗田田面相對高程實測值為試驗田微地形數據，結合遼寧地區某試驗田實例，進行固定參數的輸入，具體參數見表1—3。

表1 農田灌溉試驗田尺寸及灌溉輸入參數

(2)按照均勻設計方案對土壤入滲參數及田面糙率系數的初始值進行確定。固定參數對于田面灌溉水流模擬均為輸入的固定值，可變參數主要按照均勻設計方案進行K、a、n的調整。

3 參數優化識別模型建立

未知參數通過實測數據進行反推是參數識別的實質，函數關系式通過實測值和模擬值進行建立，二者之間的誤差最小值通過模型優化進行尋求。本文采用地面灌溉實測水流數據，不同參數組合下的地面灌溉水流推進數據采用WinSRFR模型進行模擬，結合均勻設計在地面入滲參數及田面糙率系數范圍確定基礎上進行方案優化試驗，通過對參數的不斷優化調整，使得實測值和模擬值之間的誤差值最低，再通過地面灌溉水流實測數據進行驗證，參數優化識別目標函數構建方程為：

(5)

式中，y(K，a，n)—差值平方和的目標函數；Xi—灌溉水流推進到i觀測點的距離，m;N—灌溉水流觀測點總個數。

表2 農田灌溉試驗田相對高程 單位：m

表3 農田灌溉試驗灌溉水流推進及消退測定數據

K，a，n各參數取值范圍為目標函數的約束條件：

Kmin≤K≤Kmax

(6)

amin≤a≤amax

(7)

nmin≤n≤nmax

(8)

式中，Kmin、amin、nmin—K，a，n各參數的最小取值；Kmax、amax、nmax—K，a，n各參數的最大取值。

4 均勻試驗設計

4.1 均勻試驗設計表

(1)試驗指標確定。多個試驗指標需要在方案中進行設計時，要進行綜合分析。

(2)因素的選取。影響試驗指標較大的因素需要結合專業知識或實際經驗進行挑選。均勻設計方案主要采用均勻設計試驗進行分析，地面灌溉最優參數組合通過模擬試驗進行確定。地面土壤入滲參數和田面糙率系數是土壤入滲模型主要影響因素，地面灌溉水量及灌水時間則是地面灌溉技術優化的主要影響因素。

(4)明確均勻設計表試驗方案，進行試驗。按照選取的均勻設計表及試驗因素明確均勻設計試驗方案，對相關變化的設計因素在設計表中進行設計。

(5)統計試驗分析。對試驗數據進行多元回歸統計試驗分析。

結合國內地面灌溉水流參數相關文獻數據，確定K，a，n的參數取值范圍：

15≤K≤350

(9)

0.15≤a≤0.85

(10)

0.02≤n≤0.45

(11)

表4 水平因素編碼賦值

對不同因素進行如下變換：

(12)

(13)

(14)

在因素變換的基礎上，對其均勻試驗進行方案設計，均勻設計方案見表5。

表5 均勻試驗設計方案

4.2 多元回歸分析

在均勻試驗設計方案的基礎上，采用EXCEL中多元回歸公式對均勻試驗數據進行回歸統計分析，結果見表6—7。

表6 回歸顯著檢驗及統計結果

采用均勻試驗設計方案下地面灌溉差值平方和的擬合曲線經過多元回歸分析其相關系數高于0.9，Fa臨界值低于0.05擬合效果，具有顯著的回歸效果。從表7中各回歸變量顯著性檢驗結果可看出，各變量參數t檢驗值除a2的檢驗值高于0.05外，其他變量P檢驗值均低于0.05，呈現顯著相關，將不顯著相關回歸變量a2進行剔除，得到目標函數差值平方和y與K，a，n各參數之間的回歸方程：

表7 回歸方程各變量回歸系數顯著性檢驗值

y=796.1521+180.6697K-31615.8a+60913.82n-26.6674Ka-78.2701Kn+23707.28an-0.15357K2-66825.7n2

(15)

5 參數優化識別

5.1 目標函數及約束條件

按照方程(15)建立決策變量為土壤入滲參數K、a及田面糙率系數n的參數優化模型，采用粒子群算法對參數進行識別，目標函數方程為：

miny(K，a，n)=796.1521+180.6697K-31615.8a+60913.82n-26.6674Ka-78.2701Kn+23707.28an-0.15357K2-66825.7n2

(16)

約束條件為：

15≤K≤350

(17)

0.15≤a≤0.85

(18)

0.02≤n≤0.45

(19)

5.2 參數識別步驟

(1)首先初始化參數K，a，n，搜索點速度及方位坐標進行初始化，在參數約束條件范圍內進行隨機生產，各粒子當前位置為個體極值點坐標，且對極值點粒子的適應度進行計算。各粒子個體中極值最好的為全局極值，粒子序號通過全局極值點進行記錄，該粒子當前位置坐標為全局極致點坐標。

(2)對于每個粒子的K，a，n適應度進行計算，當前粒子極值好于當前粒子值，則該粒子位置坐標位置為個體極致點坐標位置，且對個體極值進行更新。

(3)將尋求的粒子群全局最優適應度和各粒子適應度進行對比，若其當前適應度好于全局最優適應度，則該粒子位置為全局極值點坐標位置，對該粒子所處位置進行序號標記，其對全局極值進行更新，全局極值為當前粒子坐標位置。

(4)結合粒子群算法對粒子位置及收斂速度進行迭代更新。

(5)對算法是否符合收斂準則要求，若符合則輸出尋優值，否則重新回到步驟(2)循環計算。

5.3 參數優化識別結果

按照粒子群優化算法求解步驟對方程(16)進行求解計算，確定K，a，n參數的優化識別值分別為：K=84.025(mm/ha)、a=0.7571、n=0.075。

5.4 結果驗證

將參數優化識別地面灌溉3個參數K=84.025(mm/ha)、a=0.7571、n=0.075分別代入到WinSRFR模型進行地面灌溉水流模擬，并與試驗農田灌溉實測水流推進數據進行誤差對比，見表8。

表8 試驗田灌溉水流實測值和模擬值差值分析

通過WinSRFR模型設定參數優化識別的土壤入滲參數及田面糙率系數后，灌溉水流及灌溉時間等數據采用表1中數據作為模型固定參數輸入，實現不同推進時間過程下的灌溉水流推進距離的模擬，從田間觀測的灌溉水流不同推進時間下的灌溉水流模擬數據可看出，灌溉水流推進模擬值和實測值之間的差值在±2.5m范圍內，通過統計差值平方和y為55.93，均小于均勻設計試驗中各方案下的差值平方和，此外還將實測推進距離和模擬推進距離進行相關性分析，其相關系數可達到0.9以上，通過實例驗證，采用的模型方法切實、可行。

5.5 不同方法參數優化識別精度對比

為對不不同方法對土壤水入滲參數及田面糙率系數參數優化識別的精度，結合遼寧省臺安水田灌溉試驗站測定數據，分別采用5種方法進行參數優化識別，并統計各方法下的模擬值和實測值之間的相對誤差均值ARE，見表9。

表9 不同方法參數優化識別精度對比

對比5種算法確定K，a，n后進行地面灌溉水流推進距離的模擬，通過模擬值和灌溉試驗站測定數據對其相對誤差均值ARE進行統計分析，從分析結果可看出，模型下相對誤差均值ARE為13.5%，相對誤差均值好于其他4種模型算法，Kostiakov-Lewis入滲模型及Kinematic-wave模型2種模型算法由于未對識別的參數進行優化，總體相對誤差均超達到20%以上。此外本文模型相比于傳統Kostiakov模型其相對誤差均值ARE提高6.2%。傳統Kostiakov模型和YSM模型相對誤差均值ARE總體較為接近。

6 結語

(1)相比于傳統如正交試驗方式，本文采用的均勻設計試驗方法可有效降低試驗次數和周期，明顯提高模型尋優求解的收斂度，對于試驗因子較多，且希望試驗次數減少，試驗要素復雜的尋優計算較為適合。

(2)以灌溉水流推進距離差值平方和最小為目標函數，并提出了相關參數的上下約束條件，可有效降低參數輸入隨機性，有效提高參數識別的精度。

(3)對均勻設計試驗水平因素下限未進行明確，存在不足，在后續研究中還應對均勻設計水平下限進行確定，提高均勻設計試驗方案設計的合理性。