基于DP-PSO算法的灌區(qū)農(nóng)業(yè)水資源優(yōu)化配置

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

0 引 言

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

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

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

1.1 單一作物灌溉水量?jī)?yōu)化模型

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

1.2 多種作物間水量?jī)?yōu)化模型

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

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

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

(1)

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

(3)約束條件。

種植面積約束：

(2)

Ak≤αkA

(3)

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

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

?j

(4)

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

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

(5)

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

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

Ak≥0；Qk≥0；Xkj≥0

(6)

2 模型求解

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

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

2.2 信息反饋

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

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

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

2.3.1初始化粒子群群體

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

(7)

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

通過(guò)式(7)的約束，將初始粒子控制在約束條件內(nèi)，可以減少無(wú)用搜索，提高粒子群算法的收斂速度，更快更準(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為粒子的個(gè)數(shù)，n為變量維數(shù)；c1和c2分別為粒子跟蹤自己歷史最優(yōu)值的權(quán)重系數(shù)和跟蹤群體最優(yōu)值的權(quán)重系數(shù)，一般取c1=c2=2；ξ1和ξ2為[0,1]之間的隨機(jī)數(shù)[14]；w為權(quán)重系數(shù)，用來(lái)控制前面的速度對(duì)當(dāng)前速度的影響，對(duì)算法的全局搜索能力和局部搜索能力進(jìn)行平衡調(diào)整。

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

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

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

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

f(x)=

(10)

式中：符號(hào)意義同前。

2.3.4粒子群算法流程

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

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

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

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

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

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

3 實(shí)例研究

3.1 灌區(qū)概況

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

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

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

3.2 作物生育階段劃分

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 有無(wú)時(shí)段約束供水過(guò)程(90%降雨，2009年可供水)

4 結(jié) 語(yǔ)

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

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

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

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

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

[3] 崔遠(yuǎn)來(lái). 非充分灌溉條件下稻田優(yōu)化灌溉制度的研究[J]. 水利學(xué)報(bào), 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)來(lái)，李遠(yuǎn)華. 作物缺水條件下灌溉供水量最優(yōu)分配[J]. 水利學(xué)報(bào), 1997,(3):38-43.

[6] 付 強(qiáng). 推求水稻非充分灌溉下優(yōu)化灌溉制度的新方法——基于實(shí)碼加速遺傳算法的多維動(dòng)態(tài)規(guī)劃法[J]. 水利學(xué)報(bào), 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é)報(bào), 2010,(S1):134-138.

[9] 高玉芳,張展羽. 混沌人工魚群算法及其在灌區(qū)優(yōu)化配水中的應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報(bào), 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é)報(bào), 2008,(3):103-106.

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

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

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

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

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

[17] 張展羽,李壽聲,何俊生. 非充分灌溉制度設(shè)計(jì)優(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.