基于CNSGA-Ⅱ算法的新疆棉花灌溉制度多目標(biāo)優(yōu)化研究

唐曉宇，何 英，彭 亮，王 俊

(1.新疆農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，烏魯木齊 830052；2.新疆水利水電勘測設(shè)計研究院，烏魯木齊 830000)

0 引 言

在內(nèi)陸干旱地區(qū)，水資源十分短缺，降雨稀少蒸發(fā)強(qiáng)烈，灌溉是實(shí)現(xiàn)作物穩(wěn)產(chǎn)、高產(chǎn)的重要措施[1]。但由于水資源短缺及不同產(chǎn)業(yè)的用水競爭，灌溉可用水量日益減少[2]，難以保證作物充分灌溉以實(shí)現(xiàn)高產(chǎn)，人們開始從根本上尋求灌溉用水最合理的利用方式，提高灌溉用水利用率[3]。在這種情況下，學(xué)者們提出非充分灌溉，其實(shí)質(zhì)是允許作物一定程度的減產(chǎn)，從而使作物灌溉用水量和產(chǎn)量達(dá)到新的平衡，是應(yīng)對水資源短缺的優(yōu)化措施。

采用傳統(tǒng)算法對作物非充分灌溉制度進(jìn)行優(yōu)化計算時，所得到的解易陷入局部最優(yōu)而不能保證達(dá)到全局最優(yōu)[4]。因此，學(xué)者們提出了解決此缺點(diǎn)的優(yōu)化算法。付強(qiáng)等[5]建立遺傳動態(tài)規(guī)劃模型(RA-GA-DP)對三江平原井灌水稻的灌溉制度進(jìn)行優(yōu)化計算，取得了滿意的結(jié)果。王斌等[6]在對灌溉制度優(yōu)化研究中，分別采用自由搜索優(yōu)化算法(Free Search，F(xiàn)S)與傳統(tǒng)算法對優(yōu)化模型進(jìn)行求解，其優(yōu)化結(jié)果表明FS算法的性能更好，且提高了尋優(yōu)精度。楊習(xí)清等[7]利用混合粒子群算法 (SELPSO)對景泰川一期灌區(qū)春小麥的灌溉制度進(jìn)行優(yōu)化研究，結(jié)果表明該算法相較于遺傳算法尋優(yōu)結(jié)果更佳。以上3種算法都沒有涉及將灌溉日期考慮在決策變量中，難以指導(dǎo)農(nóng)田灌溉實(shí)踐。尚松浩[8]采用單純形搜索法，以灌水日期為決策變量，很好地解決了農(nóng)田灌溉實(shí)踐指導(dǎo)問題，但單純形搜索法作為一種局部搜索方法，計算所得解為局部最優(yōu)解，存在一定的缺陷性?；糗娷姷萚9]以灌水日期作為決策變量，將模擬技術(shù)及遺傳算法應(yīng)用于冬小麥灌溉制度優(yōu)化中，結(jié)果表明相較于傳統(tǒng)算法，遺傳算法能得到全局最優(yōu)解，且結(jié)果穩(wěn)定。郄志紅等[10]基于改進(jìn)的NSGA-Ⅱ算法，以灌水日期和灌溉水量作為決策變量，建立了能夠同時對灌水日期和灌溉水量進(jìn)行優(yōu)化的多目標(biāo)模型，優(yōu)化目標(biāo)由以往單目標(biāo)優(yōu)化模型轉(zhuǎn)變成多目標(biāo)優(yōu)化模型，但其考慮的約束條件過少，易導(dǎo)致灌水日期的選擇對作物生長產(chǎn)生不利影響。

本文以新疆阿瓦提豐收灌區(qū)棉花灌溉制度為研究對象，以農(nóng)田水量平衡模型和作物水分生產(chǎn)函數(shù)為理論依據(jù)，以日為時間步長，以作物的灌水時間和灌水定額作為自變量，將作物的相對產(chǎn)量和全生育期總灌水量作為目標(biāo)函數(shù)，構(gòu)建內(nèi)陸干旱區(qū)作物灌溉制度多目標(biāo)優(yōu)化模型，采用CNSGA-Ⅱ算法對優(yōu)化模型求解，最終確定了該地區(qū)棉花最優(yōu)的灌水方案。

1 材料與方法

1.1 研究區(qū)概況

豐收灌區(qū)位于新疆阿克蘇地區(qū)阿瓦提縣的西南部，北緯40°20′～40°38′，東經(jīng)80°10′～80°20′，處于亞歐大陸腹地的塔里木盆地西北邊緣[11]。該地區(qū)屬暖溫帶大陸性荒漠氣候，夏季炎熱，冬季寒冷，晝夜溫差較大，光照充足，降雨稀少，蒸發(fā)強(qiáng)烈。年蒸發(fā)量1 905.2 mm，年平均氣溫10.5 ℃，年平均降水量47.8 mm，主要集中在6-8月，約占全年降水量的60%[12]。灌區(qū)主要種植作物為棉花、玉米、小麥等，灌區(qū)棉花生育期的灌溉制度如表1所示。

表1 豐收灌區(qū)棉花的灌溉制度Tab.1 Irrigation system of cotton in Fengshou area

1.2 農(nóng)田水量平衡模型

農(nóng)田水量平衡模型是用于描述旱作物在生育期任一時段內(nèi)，根系層(土壤計劃濕潤層)內(nèi)的水量平衡關(guān)系[13]，其基本表達(dá)式如下：

Wt-W0=WT+P0+M+K-E-T

(1)

式中：Wt、W0分別為根系層在時段末與時段初的儲水量；WT為隨計劃濕潤層增加而增加的水量；P0為根系層內(nèi)保存的有效降水量；K為時段t內(nèi)的地下水補(bǔ)給量；M為時段t內(nèi)的灌溉水量；E、T為時段t內(nèi)作物實(shí)際的蒸發(fā)量和蒸騰量。

農(nóng)田中主要的水分消耗項(xiàng)是作物需水量，其主要受大氣蒸發(fā)能力(通常以參考作物蒸發(fā)蒸騰量ET0表示)、生長狀態(tài)(以作物系數(shù)Kc表示)及土壤水分條件(由土壤水分脅迫系數(shù)Ks表示)等影響。作物蒸發(fā)蒸騰量(作物需水量)可以采用作物系數(shù)法計算，即：

ET=KsETm

(2)

ETm=KcET0

(3)

式中：ETm為供水充足時的田間最大蒸發(fā)蒸騰量；K2為土壤水分脅迫系數(shù)；ET0為參考作物蒸發(fā)蒸騰量；Kc為作物系數(shù)。

Kc與作物種類、生育期和作物的群體葉面積指數(shù)等因素有關(guān)，可以將其概化為生育期內(nèi)時間t的函數(shù)[14]，即：

Kc=Kcmexp{-[(t-tm)2/c2]}

(4)

式中：Kcm為生育期內(nèi)最大作物系數(shù)；tm為最大作物系數(shù)所對應(yīng)的時間；c為形狀參數(shù)。

采用FAO推薦的方法計算Ks，根區(qū)土壤含水量用根區(qū)耗水量Dr表示[15]。Ks計算式如下：

(5)

式中：TAW為土壤總有效水量；Dr為作物根區(qū)消耗水量；RAW為土壤易利用有效水分。

TAW采用下式計算：

TAW=1 000 (θFC-θWP)Zr

(6)

式中：θFC為田間持水量；θWP為永久凋萎點(diǎn)含水量；Zr為根區(qū)深度。

RAW的計算式如下：

RAW=ρTAW

(7)

式中：ρ為水分脅迫發(fā)生前土壤易利用有效水分與土壤總有效水量的比值。

1.3 水分生產(chǎn)函數(shù)模型

1.3.1 模型的選取

作物水分生產(chǎn)函數(shù)表示田間水分消耗與作物產(chǎn)量之間的關(guān)系[16]。本研究采用Jensen乘法模型，具體表達(dá)式如下：

(8)

式中：Y為作物實(shí)際產(chǎn)量；Ym為作物最高產(chǎn)量；EM、TM為作物潛在蒸發(fā)量和蒸騰量；n為作物生育期的總階數(shù)；i為作物生育階段劃分序號；λi為作物在第i生育階段的缺水敏感指數(shù)。

1.3.2 參數(shù)求解

傳統(tǒng)的回歸分析法在求解水分敏感指數(shù)時，所得解不穩(wěn)定，且經(jīng)常出現(xiàn)負(fù)值[17]，為避免出現(xiàn)負(fù)值，影響計算的精確性，本文采取水分敏感指數(shù)累積函數(shù)進(jìn)行求解[18]：

(9)

λ(Δti)=z(ti)-z(ti-1)

(10)

式中：ti為作物從某一指定日期算起的生長天數(shù)；λ(Δti)為時段Δt=ti-ti-1的水分敏感指數(shù)值；z(t)為水分敏感指數(shù)累積曲線，a、b、c為待求系數(shù)。根據(jù)試驗(yàn)資料[19]，擬合得阿克蘇棉花a=5.0072、b=0.0497、c=1.23。

1.3.3 水分敏感指數(shù)隨時間的變化過程

根據(jù)上述所求得的a、b、c系數(shù)值，可以推算阿克蘇地區(qū)棉花水分敏感指數(shù)隨時間的變化過程，如圖1和圖2所示。根據(jù)當(dāng)?shù)孛藁ǖ纳趧澐?，得出各階段的缺水敏感指數(shù)：λ1=0.186 4、λ2=0.202 9、λ3=0.678、λ4=0.135 8。

圖1 棉花水分敏感指數(shù)累計函數(shù)變化過程Fig.1 The change process of the cumulative function of cotton moisture sensitivity index

圖2 棉花時段水分敏感指數(shù)變化過程Fig.2 Change process of water sensitive index in cotton

2 優(yōu)化模型與優(yōu)化算法

2.1 優(yōu)化模型

在非充分灌溉條件下，以作物相對產(chǎn)量最大和全生育期總灌水量最小為優(yōu)化目標(biāo)，建立以下優(yōu)化模型。

(1)目標(biāo)函數(shù)。

(11)

(12)

當(dāng)為作物根系提供吸水作用的土層含水量達(dá)到凋萎系數(shù)含水量時，作物相對產(chǎn)量為零。

(2)約束條件。

灌水日期約束：

dmin≤dj≤dmax

(13)

土壤計劃濕潤層深度的約束：

Hmini≤Hi≤Hmaxi

(14)

土壤含水率約束：

θwp≤θi≤θf

(15)

式中：f1為作物生育期總灌水量，mm；f2為作物相對產(chǎn)量；Ei、Ti為作物在第i時段的蒸發(fā)量和蒸騰量，mm；EMi、TMi為作物在第i時段潛在蒸發(fā)量和蒸騰量，mm；xj為作物在第j次灌水量，mm；dmin、dmax為作物生育期的最小和最大天數(shù)，d；dj為作物從生育期開始到第j次灌水時的天數(shù)，mm；Hmini、Hmaxi為土壤在作物第i生長階段內(nèi)最小和最大的濕潤層深度，mm；θi為土壤第i階段含水率；θwp和θf為凋萎系數(shù)和田間持水率。

2.2 優(yōu)化算法

1994年Deb和Srinivas提出非支配排序遺傳算法(Non-dominated Sorting Genetic Algorithm，NSGA)[20]。該算法采用非支配分層的思想，利用適應(yīng)度共享策略使得準(zhǔn)Pareto面上的個體均勻分布，解決了算法早熟收斂的問題。2002年，Deb等[21]提出了快速非支配排序算法(Non-dominated Sorting Genetic Algorithm-II，NSGA-Ⅱ)。該算法主要在原NSGA算法的基礎(chǔ)上引入精英策略，使用擁擠度概念對非支配層級中同一級的個體進(jìn)行排序，提高了算法速度。

NSGA-Ⅱ算法主要用于求解無約束多目標(biāo)優(yōu)化問題，而本研究是以總灌水量最小及相對產(chǎn)量大為目標(biāo)函數(shù)的多約束多目標(biāo)優(yōu)化模型，所以將原算法進(jìn)行改進(jìn)，將Deb約束準(zhǔn)則引入到NSGA-Ⅱ算法中，記作CNSGA-II算法。

Deb約束支配準(zhǔn)則中每個個體根據(jù)Pareto支配、約束違反度兩種信息來選擇較優(yōu)個體。CNSGA-Ⅱ算法中約束處理的方法：如對于一個解x，若其滿足約束條件，則稱該解為可行解，若不滿足，則稱之為不可行解。對于不可行解，一般使用約束違反值(Constraint Violation value，CV)，該值用來定量描述一個解違反約束條件的程度。對于一個解x，其值可如下表達(dá)：

(16)

其中當(dāng)g(x)<0時，g(x)= -g(x),否則g(x)=0，當(dāng)x滿足約束條件中任何條件，即x在可行域內(nèi)，則CV(x)=0；當(dāng)x不完全滿足約束條件時，即x不在可行域內(nèi)，則CV(x)≠0，當(dāng)CV值越小，x越靠近可行域。

3 模型求解及結(jié)果分析

3.1 模型參數(shù)設(shè)定

本研究以豐收灌區(qū)棉花為例，該灌區(qū)棉花的生育期從4月15日到8月20日，共138 d。已收集阿瓦提縣氣象站(東經(jīng)80°24′，北緯40°39′，海拔1045.8m)1960年到2001年共42年的降水資料。模型參數(shù)，如：土壤初始含水量W、田間持水量θFC、作物缺水敏感指數(shù)λi，生育期潛在騰發(fā)量ETM、生育期實(shí)際騰發(fā)量ET。其中ET、ETM分別由公式(2)、(3)計算。

本研究假定棉花在生育期總共灌5次水，將每次灌水量xi與灌水間隔ti作為優(yōu)化參數(shù)，灌水時間間隔設(shè)置為整數(shù)。算法中參數(shù)選?。悍N群規(guī)模P=100，最大進(jìn)化代數(shù)g=500，運(yùn)行次數(shù)20，變異概率Pm=0.1，交叉概率Pc=0.9，交叉分布指數(shù)etac=20，突變分布指數(shù)etam=80。

3.2 優(yōu)化結(jié)果與分析

采用CNSGA-II算法對灌溉制度多目標(biāo)優(yōu)化模型進(jìn)行求解。為提高模型的適用性，本文選取不同典型年的降水優(yōu)化相應(yīng)灌溉制度。典型年的選取方法：收集當(dāng)?shù)?960-2001年共42 a的降水資料，按照總降水量在作物全生育期的大小進(jìn)行頻率計算，選取1996年(生育期降水量109.5 mm)作為豐水年(頻率P=10%)；選取1992年(生育期降水量52.7 mm)作為平水年(頻率P=50%)；選取1973年(生育期降水量13.2 mm)作為枯水年(頻率P=90%)。典型年灌水總量與相對產(chǎn)量的優(yōu)化結(jié)果如圖3～圖5所示。

圖3 枯水年相對產(chǎn)量關(guān)系優(yōu)化結(jié)果Fig.3 Optimal results of the relative output relationship in dry years

圖4 平水年相對產(chǎn)量關(guān)系優(yōu)化結(jié)果Fig.4 Optimal results of the relative output relationship in the average year

圖5 豐水年相對產(chǎn)量關(guān)系優(yōu)化結(jié)果Fig.5 Optimal results of the relationship between the relative output in high water years

由圖3～圖5可知，豐水年、平水年和枯水年圖形中優(yōu)化解的變化趨勢大致相同，在一定范圍內(nèi)，棉花的產(chǎn)量都隨著灌水量的增加而顯著增長。但在平水年，當(dāng)棉花的總灌水量達(dá)到385 mm，其相對產(chǎn)量增加的幅度明顯減小且隨著灌水量的增加而趨近于直線，即灌水所帶來的效益微乎其微，而在豐水年需要總灌水量352 mm、枯水年則需總灌水量416 mm才能達(dá)到相同的效果。

在3個典型年的相對產(chǎn)量穩(wěn)定不變時，他們的相對產(chǎn)量分別為：豐水年0.954、平水年0.929、枯水年0.926，說明降水對平水年與枯水年的相對產(chǎn)量影響不大。在選取典型年為平水年時，該模型相較于當(dāng)?shù)匾酝贫ǖ墓喔戎贫龋涸诠喔榷~一致時，作物相對產(chǎn)量提高了2.5%～3.1%；在相對產(chǎn)量一致時，每公頃土地的灌溉用水節(jié)約了495～705 m3，表明對灌水日期和灌溉水量進(jìn)行合理的安排能在保證作物產(chǎn)量穩(wěn)定的情況下節(jié)約灌溉用水，對于作物生長依靠灌溉的缺水地區(qū)，該模型的合理利用能使農(nóng)戶在作物產(chǎn)量和總灌水量不變的情況下，增加作物種植面積，提高當(dāng)?shù)剞r(nóng)戶的經(jīng)濟(jì)收入。從平水年的優(yōu)化結(jié)果中選取部分Pareto解如表2所示。

4 結(jié) 論

本研究利用農(nóng)田水量平衡模型模擬作物根系層水量之間的關(guān)系，采用作物水分生產(chǎn)函數(shù)模型計算其相對產(chǎn)量，以日為步長，以作物的灌水時間和灌水定額作為自變量，將作物的相對產(chǎn)量和全生育期總灌水量作為目標(biāo)函數(shù)，構(gòu)建內(nèi)陸干旱區(qū)作物灌溉制度多目標(biāo)優(yōu)化模型；為處理含約束的多目標(biāo)優(yōu)化問題，將Deb約束準(zhǔn)則引入到NSGA-Ⅱ算法中，記作CNSGA-Ⅱ算法，并利用CNSGA-Ⅱ算法對優(yōu)化模型進(jìn)行求解，得到了實(shí)驗(yàn)地區(qū)最優(yōu)的灌溉制度。

表2 棉花在灌溉制度優(yōu)化下的部分Pareto解Tab.2 Part of Pareto solution of cotton under the optimization of irrigation system

通過對棉花在3種典型年的灌溉制度優(yōu)化表明：在內(nèi)陸干旱區(qū)作物的生長主要依靠農(nóng)業(yè)灌溉，降水對其相對產(chǎn)量的影響不大。在典型年為平水年時，相較于當(dāng)?shù)匾酝墓喔戎贫龋涸诠喔榷~一致時，作物產(chǎn)量提高了2.5%～3.1%；在相對產(chǎn)量一致時，每公頃土地的灌溉用水節(jié)約了495～705 m3。

棉花灌溉制度多目標(biāo)優(yōu)化模型在優(yōu)化灌水量的同時也將灌水日期一并優(yōu)化，相較于以往只考慮單目標(biāo)的優(yōu)化模型，其優(yōu)化結(jié)果更加合理，而灌水時間間隔設(shè)置為整數(shù)，便于指導(dǎo)農(nóng)業(yè)灌溉實(shí)踐。目前模型輸入的降水?dāng)?shù)據(jù)為歷史數(shù)據(jù)，具有一定的局限性，為制定適用性更強(qiáng)的灌溉制度，在進(jìn)一步研究中將預(yù)測未來氣候變化的模型與灌溉制度優(yōu)化模型相結(jié)合，既能預(yù)測未來極端天氣，又能動態(tài)調(diào)整灌水量與灌水日期。