劉 博，崔遠(yuǎn)來，尹杰杰，劉方平，鄭傳舉(.武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢 4007；.江西省灌溉試驗中心站，南昌 00；.湖北省樊口電排站管理處，湖北 鄂州 4600)

0 引 言

灌區(qū)農(nóng)業(yè)水資源優(yōu)化配置一般包括單一作物灌溉制度優(yōu)化、作物種植結(jié)構(gòu)優(yōu)化和作物間最優(yōu)水量分配，前者解決灌溉水量在農(nóng)作物生育期間的分配問題，后者解決不同作物之間水量最優(yōu)分配以及作物種植結(jié)構(gòu)優(yōu)化問題。針對以上問題，國內(nèi)外學(xué)者做了大量研究。袁宏源[1]等基于DP方法求解了北方多種旱作物的最優(yōu)灌溉制度。郭宗樓[2]采用非線性規(guī)劃(NLP)，從作物騰發(fā)和根區(qū)土壤水運動規(guī)律出發(fā)，求解了作物最優(yōu)灌溉制度。崔遠(yuǎn)來[3]等基于以上模型，融合稻田水量平衡理論得到非充分灌溉條件下稻田最優(yōu)灌溉制度。Yaron[4]等利用LP-DP模型求解多種作物間灌溉水量的最優(yōu)分配。崔遠(yuǎn)來[5]等運用DP-DP模型對多作物間灌溉水量最優(yōu)分配進(jìn)行求解。近年來很多智能算法在灌溉水量最優(yōu)分配得到廣泛應(yīng)用，付強[6]等將改進(jìn)加速遺傳算法(RAGA)與多維動態(tài)規(guī)劃法(DP)相結(jié)合，潘琳[7]等采用模擬退火遺傳算法求解灌區(qū)多種作物間灌溉水量的最優(yōu)分配問題。楊娜[8]等采用連續(xù)蟻群算法求解水稻灌溉制度。但這些模型中，多將總灌溉水量視為可調(diào)控水量，對于引水灌區(qū)，大多沒有調(diào)蓄能力，各時段可供水量是一定的，以上模型并不適用，本文針對引水灌區(qū)，同時考慮灌溉用水總量約束和灌區(qū)時段供水量約束，以農(nóng)作物產(chǎn)值最大為目標(biāo)，建立了基于DP-PSO灌區(qū)農(nóng)業(yè)水資源優(yōu)化配置模型。

1 灌溉水量最優(yōu)分配的分解協(xié)調(diào)模型

在灌區(qū)總灌溉水量不足條件下，可針對灌區(qū)總效益最大的目標(biāo)，將有限的水量在不同作物之間進(jìn)行最優(yōu)分配，不同作物又將分配得到的水量在不同生育階段之間進(jìn)行最優(yōu)分配，從而確定水源水量最優(yōu)配水過程。將作物子系統(tǒng)作為第一層，灌區(qū)總系統(tǒng)作為第二層，通過分配給每種作物的供水量將兩層聯(lián)系起來，則成為一個具有兩層譜系結(jié)構(gòu)的大系統(tǒng)，適合用大系統(tǒng)分解協(xié)調(diào)模型求解[5]。

1.1 單一作物灌溉水量優(yōu)化模型

以旬為時段，作物水分生長函數(shù)采用Jensen模型，以相對產(chǎn)量最大為目標(biāo)，采用動態(tài)規(guī)劃建模求解，模型參照文獻(xiàn)[5]。

1.2 多種作物間水量優(yōu)化模型

在灌區(qū)總灌水量不足情況下，為達(dá)到最大效益，在不同時期或者同一時期內(nèi)種植有多種農(nóng)作物的灌區(qū)，需要將有限的灌溉水量在不同作物之間進(jìn)行最優(yōu)分配，同時確定各種作物種植面積和灌水量，模型如下：

(1)變量和時段劃分。第k種作物所分配的總凈灌溉水量Qk(萬m3)和灌溉面積Ak(hm2)為自變量，作物生育期按旬劃分成j個時段，Xkj表示第k種作物在j時段所分配凈灌水量，m3/hm2。

(2)目標(biāo)函數(shù)：以灌區(qū)總產(chǎn)值最大為目標(biāo)。

(1)

式中：YMk為第k種作物充分供水條件下產(chǎn)量，kg/hm2；PRk為第k種作物產(chǎn)品單價，元/kg；G為各作物效益之和，萬元。

(3)約束條件。

種植面積約束：

(2)

Ak≤αkA

(3)

式中：A為總可灌溉面積，hm2；Ak為第k種作物灌溉面積，hm2；αk為第k種作物最大可灌溉面積占總灌溉面積的百分比。

時段引水量約束：各時段渠首引水量(毛灌溉用水量)不得超過同時段可供水量，即：

?j

(4)

式中：Xkj為第k種作物第j時段的凈灌水量，m3/hm2；Ak為第k種作物灌溉面積，hm2；η為灌溉水利用系數(shù)；Wj為第j時段可供水量，萬m3。

總水量約束：總灌溉需水量不得超過灌區(qū)分配給作物灌溉的總可供水量，即：

(5)

式中：W0為灌區(qū)分配給作物灌溉的總可供水量，萬m3。

非負(fù)約束：各變量非負(fù)，即：

Ak≥0；Qk≥0；Xkj≥0

(6)

2 模型求解

2.1 下層模型優(yōu)化

模型運行時，首先對下層單一作物模型求解，對第k種作物擬定灌溉定額qk，各子系統(tǒng)在給定的qk下，采用DP分別獨立優(yōu)化，得到最優(yōu)效益F(qk)和各階段灌水定額Xkj，具體求解過程見文獻(xiàn)[5]。

2.2 信息反饋

將下層模型的最優(yōu)效益F(qk)和各階段的灌水定額Xkj作為上層模型目標(biāo)函數(shù)和約束條件的變量值，反饋至上層。

2.3 上層模型優(yōu)化

粒子群算法具有參數(shù)少、計算簡便、收斂速度快、計算精度高等優(yōu)點，同時可以通過對參數(shù)的設(shè)置來平衡算法的局部收斂能力和全局收斂能力。作物種植結(jié)構(gòu)及作物間灌溉水量最優(yōu)分配模型采用粒子群算法求解。

2.3.1初始化粒子群群體

對于既存在總量約束，又存在單獨變量上下限約束的問題，設(shè)n維空間上的第i個粒子表示為Xi=(xi1,xi2,…,xin)，則可由以下方法生成初始粒子群體：

(7)

式中：Q為總量約束，即∑nd=1xid≤Q；XL=(xL1,xL2,…,xLn)為變量的上限；XU=(xU1,xU2,…,xUn)為變量的下限；u為[0,1]上的隨機數(shù)。

通過式(7)的約束，將初始粒子控制在約束條件內(nèi)，可以減少無用搜索，提高粒子群算法的收斂速度，更快更準(zhǔn)得找到最優(yōu)解。

2.3.2遞推公式及參數(shù)設(shè)置

粒子速度及位置更新迭代公式如下：

xk+1id=xkid+vk+1id

(8)

vk+1id=wvkid+c1ξ1(xpkid-xkid)+

c2ξ2(xgkd-xkid)

(9)

式中：i=1,2,…,m；i=1,2,…,n，其中m為粒子的個數(shù)，n為變量維數(shù)；c1和c2分別為粒子跟蹤自己歷史最優(yōu)值的權(quán)重系數(shù)和跟蹤群體最優(yōu)值的權(quán)重系數(shù)，一般取c1=c2=2；ξ1和ξ2為[0,1]之間的隨機數(shù)[14]；w為權(quán)重系數(shù)，用來控制前面的速度對當(dāng)前速度的影響，對算法的全局搜索能力和局部搜索能力進(jìn)行平衡調(diào)整。

Y.Shi提出w的線性遞減策略[10]，即隨著迭代的進(jìn)行，線性減少權(quán)重w的值，這種策略能夠兼顧搜索效率和搜索精度，改善優(yōu)化性能，Y.Shi給出的遞減范圍是1.4到0.4。

為保證每次迭代的群體都滿足總量約束，需要對速度vki進(jìn)行設(shè)置。陳曉楠[11]提出，每次迭代時，若采用式(8)計算的xk+1i不滿足總量約束，則對vki乘以系數(shù)u，u為[0,1]之間的隨機數(shù)，重新計算xk+1i，直到滿足要求為止。

2.3.3適應(yīng)度函數(shù)設(shè)置。

適應(yīng)度函數(shù)一般由目標(biāo)函數(shù)變化而來，即：

f(x)=

(10)

式中：符號意義同前。

2.3.4粒子群算法流程

(1)初始化粒子群的種群規(guī)模、位置和速度，設(shè)置最大的迭代次數(shù)T；

(2)對于每一個粒子i，計算其適應(yīng)值，然后根據(jù)適應(yīng)度值確定粒子i的個體最優(yōu)值pbesti和全局最優(yōu)值gbesti；

(3)對于每一個粒子i，將其當(dāng)前適應(yīng)值與先前個體最優(yōu)值pbesti比較，如果前者優(yōu)于后者，則替換pbesti，否則保留原個體最優(yōu)值；

(4)對于每一個粒子i，用其當(dāng)前個體最優(yōu)值與種群全局最優(yōu)值gbesti比較，如果較好，則將其作為當(dāng)前全局最優(yōu)值，否則保留原全局最優(yōu)值；

(5)根據(jù)公式(8)、式(9)更新粒子的速度和位置；

(6)如果滿足結(jié)束條件，退出，否則轉(zhuǎn)到(2)。

3 實例研究

3.1 灌區(qū)概況

贛撫平原灌區(qū)地處江西省贛中偏北部，跨南昌縣、進(jìn)賢縣、豐城市、青山湖區(qū)等地域，總灌溉面積68 980 hm2，主要灌溉作物種類為早、中、晚稻，相關(guān)信息見表1。除水稻外，灌區(qū)還種植豆類、油菜、蔬菜等旱作物，由于旱作物種植面積和需水量均較少，且品種多樣、經(jīng)濟效益較高，假定灌區(qū)旱作物的種植面積和灌溉用水需求完全滿足。

2012年《江西省人民政府關(guān)于實行最嚴(yán)格水資源管理制度的實施意見》“贛府發(fā)[2012]29號”文規(guī)定，贛撫平原灌區(qū)從2015年起用水總量控制為10.08億m3?，F(xiàn)狀灌區(qū)引水量25億m3遠(yuǎn)大于該分配定額，因此未來限額供水條件下研究灌區(qū)農(nóng)業(yè)灌溉用水優(yōu)化配置十分必要。

表1 贛撫平原灌區(qū)水稻有關(guān)基本資料

3.2 作物生育階段劃分

以旬為時段，贛撫平原灌區(qū)水稻生育階段及DP模型求解最優(yōu)灌溉制度的相關(guān)參數(shù)見表2。

表2 水稻各生育階段相關(guān)參數(shù)

3.3 灌區(qū)可供水量計算及降雨頻率選取

灌區(qū)引水工程有兩處，分別為焦石進(jìn)水閘和柴埠口進(jìn)水閘，兩進(jìn)水閘最大引水流量為249 m3/s。結(jié)合1979-2014年撫河天然徑流資料，分旬計算，可得1979-2014年各旬可供水量。

灌區(qū)降雨和渠道可供水量是影響灌區(qū)農(nóng)作物灌溉用水的兩大主要因素，由于灌區(qū)降雨和可供水量存在不同步現(xiàn)象，為充分考慮兩方面的影響，根據(jù)南昌站多年實測資料統(tǒng)計，按水稻生育期降雨量總量排頻，選取90%、75%和50% 3個降雨頻率，針對每個降雨頻率，分別計算1979-2014年渠首不同可供水條件下灌區(qū)農(nóng)業(yè)配水方案。

3.4 灌區(qū)農(nóng)業(yè)水資源優(yōu)化配置結(jié)果

假定灌區(qū)生活、環(huán)境、工業(yè)、經(jīng)濟作物及水稻秧田泡田期水量完全滿足，計算得到灌區(qū)可用于水稻本田生育期的凈灌溉水量為2.5億m3，采用以上DP-PSO模型求解3種不同降雨頻率下多年平均產(chǎn)值及水稻種植面積優(yōu)化結(jié)果，見表3，以90%降雨頻率為例，各年產(chǎn)值及水稻種植面積優(yōu)化結(jié)果見圖1。若不考慮各時段供水量約束，3個降雨頻率下多年平均產(chǎn)值及水稻種植面積見表4。

表3 考慮時段供水量約束時多年平均產(chǎn)值及水稻種植面積優(yōu)化結(jié)果

表4 不考慮時段供水量約束時多年平均產(chǎn)值及水稻種植面積優(yōu)化結(jié)果

(1)由表4可見，不考慮時段供水約束時，優(yōu)化配水后，早稻和晚稻種植面積都達(dá)到最大可種植面積，而中稻種植面積由現(xiàn)狀的7 333 hm2減至0，這是因為中稻生育期耗水量大，而豐產(chǎn)產(chǎn)量及單價僅略高于早晚稻，水的投入產(chǎn)出比低，因此會盡量削減中稻面積而改種早晚稻。3種降雨頻率年下，隨降雨量增加，缺水程度降低，因此雖然早晚稻種植面積相同，但產(chǎn)值呈增加趨勢。

(2)考慮時段供水約束后，由圖1可見，部分年份無法使早晚稻達(dá)到最大可種植面積，這是因為這些年份早稻或晚稻生育期間存在時段供水約束，無法滿足最大面積的灌溉需求。以1986年可供水為例，5月下旬、8月中下旬和9月中旬至10月中旬渠首可供水量均不足1 000萬m3，造成早晚稻種植面積的削減，在90%降雨頻率下早晚稻種植面積分別為4.51和5.11萬hm2，而在中稻敏感指數(shù)最大的9月上旬可供水量為6 000萬m3，因此會考慮種植部分中稻以提高整體效益，中稻種植面積為0.49萬hm2。

圖1 90%降雨頻率下各年產(chǎn)值及水稻種植面積優(yōu)化結(jié)果

(3)從表3可以看到，90%、75%和50%降雨頻率下，隨降雨的增加多年平均總種植面積、早晚稻種植面積均呈增加趨勢，而中稻面積正好相反，因為在降雨量增加情況下，從渠道取用的灌溉需水量減少，時段可供水量約束對早晚稻面積限制削減，會盡量種植水分生產(chǎn)率高的早晚稻，而減少中稻種植面積；3種頻率年下多年平均總灌溉需水量均未達(dá)到2.5億m3且呈現(xiàn)依次遞減趨勢，表明有很多時段灌溉需水與渠道可供水不同步。

(4)以90%降雨頻率下，2009年可供水為例，分析有無考慮時段供水量約束情況下的實際灌溉需水過程，見圖2。不考慮時段供水約束時，9月份灌溉需水量最多，但9月上、中旬均超出水源可供水能力，此外還有6月中旬超出該旬供水能力?？紤]時段供水約束后，9月份需水量明顯減少，而8月份需水量增加，特別是8月中旬，6月上旬需水量同樣增加。計算結(jié)果表明，若按不考慮時段供水約束實際配水，則到6月中旬、9月上中旬會存在供水不足的情況，導(dǎo)致灌區(qū)總體產(chǎn)值降低，為23.23億元，而考慮供水約束后，改變了灌溉需水分配過程，使需水與可供水更匹配，灌區(qū)總體產(chǎn)值為23.60億元，提高了水資源利用效率。

圖2 有無時段約束供水過程(90%降雨，2009年可供水)

4 結(jié) 語

(1)采用DP-PSO的農(nóng)業(yè)水資源優(yōu)化配置模型，充分考慮了灌溉用水總量約束和時段可供水量約束，適用于水源可供水過程有約束條件下農(nóng)業(yè)水資源優(yōu)化配置研究。

(2)對于多重約束問題，提出了一種提高收斂速度的初始化粒子群方法，并通過粒子速度的動態(tài)變化來保證每代粒子都滿足約束。這些方法減少了無用搜索，提高粒子群算法的收斂速度及精度。

(3)針對贛撫平原灌區(qū)，考慮灌區(qū)普遍存在的降雨和水源可供水量不同步現(xiàn)象，計算了3種降雨頻率與不同年可供水組合情況的農(nóng)業(yè)水資源優(yōu)化配置方案，分析了有無時段供水約束下作物種植結(jié)構(gòu)及水資源配置方案，結(jié)果表明有時段供水約束下的優(yōu)化結(jié)果更符合引水灌區(qū)實際。

[1] 袁宏源，劉肇祎. 高產(chǎn)省水灌溉制度優(yōu)化模型研究[J]. 水利學(xué)報, 1990,(11):1-7.

[2] 郭宗樓. 灌溉水資源最優(yōu)分配的DP-DP法[J]. 水科學(xué)進(jìn)展, 1994,(4):303-308.

[3] 崔遠(yuǎn)來. 非充分灌溉條件下稻田優(yōu)化灌溉制度的研究[J]. 水利學(xué)報, 1995,(10):29-34.

[4] Yaron D, A Dinar, S Meyers. Irrigation scheduling-Theoretical approach and application problems[J]. Water Resources Management, 1987,1(1):17-31.

[5] 崔遠(yuǎn)來，李遠(yuǎn)華. 作物缺水條件下灌溉供水量最優(yōu)分配[J]. 水利學(xué)報, 1997,(3):38-43.

[6] 付 強. 推求水稻非充分灌溉下優(yōu)化灌溉制度的新方法——基于實碼加速遺傳算法的多維動態(tài)規(guī)劃法[J]. 水利學(xué)報, 2003,(1):123-128.

[7] 潘 林, 肖新棉,雷永富. 模擬退火遺傳算法在灌溉水量最優(yōu)分配中的應(yīng)用[J]. 節(jié)水灌溉, 2006,(6):45-47,50.

[8] 楊 娜. 連續(xù)蟻群算法在水稻灌溉制度優(yōu)化中的應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報, 2010,(S1):134-138.

[9] 高玉芳,張展羽. 混沌人工魚群算法及其在灌區(qū)優(yōu)化配水中的應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報, 2007,(6):7-11.

[10] ShiY, EberhartR. A Modified Partiele Swann OPtimizer[C]∥ IEEE World Congresson Computational Intelligenee, Anehorage, AK, USA,1998:69-73.

[11] 陳曉楠. 基于粒子群的大系統(tǒng)優(yōu)化模型在灌區(qū)水資源優(yōu)化配置中的應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報, 2008,(3):103-106.

[12] 李 智,鄭 曉. 粒子群算法在農(nóng)業(yè)工程優(yōu)化設(shè)計中的應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報, 2004,(3):15-18.

[13] 劉建華. 粒子群算法的基本理論及其改進(jìn)研究[D]. 長沙：中南大學(xué), 2009.

[14] 崔遠(yuǎn)來,袁宏源,李遠(yuǎn)華.考慮隨機降雨時稻田高效節(jié)水灌溉制度[J].水利學(xué)報, 1999,(7):40-45.

[15] 榮豐濤,王仰仁. 山西省主要農(nóng)作物水分生產(chǎn)函數(shù)中參數(shù)的試驗研究[J].水利學(xué)報,1997,(1):78-83.

[16] 周祖昊,袁宏源. 有限供水條件下灌區(qū)優(yōu)化配水[J].中國農(nóng)村水利水電, 2002,(5):5-7.

[17] 張展羽,李壽聲,何俊生. 非充分灌溉制度設(shè)計優(yōu)化模型[J].水科學(xué)進(jìn)展,1993,(3):207-214.

[18] Kumar R, khepar S D. Decision model for optimal cropping patterns in irrigation based on crop water production function[J]. Agric. Water Manage., 1980,(3):65-76.

[19] Raju K S, Stanley Lee E, Biere A W, et al. Irrigation scheduling based on a dynamic crop response model[J]. Advances in Irrigation, 1983,(2):257-271.

[20] Tsakiris G P. A method for applying crop sensitivity factors in irrigation scheduling[J]. Agricultural Water Management, 1982,5(4):335-343.

[21] HuX, ShiY, Eberhart R. Recent advancesin Particle swarm[C]∥ Congresson Evolutionary Computation, 2004:19-23.