基于動態(tài)規(guī)劃模型的水資源優(yōu)化配置
——以泰來灌區(qū)為例

徐淑琴，徐恩典，馬圣潔

（東北農(nóng)業(yè)大學水利與土木工程學院，哈爾濱 150030）

近年來，我國針對大型灌區(qū)普遍存在的設計標準低、配套程度差、老化損壞嚴重、用水效率低、管理體制不順以及經(jīng)營體制不活等問題進行系統(tǒng)整改，以促進我國農(nóng)業(yè)發(fā)展、農(nóng)村經(jīng)濟繁榮和農(nóng)民收入提高。特別是我國北方灌區(qū)，存在溝渠沖淤、骨干渠道襯砌不足等問題，水資源嚴重浪費，影響灌區(qū)總體效益發(fā)揮。為保證灌區(qū)正常運行，發(fā)揮自然優(yōu)勢，提高糧食產(chǎn)量和質量，近年來對大型灌區(qū)骨干工程實施配套和節(jié)水改造政策，單位面積凈效益及單位水量經(jīng)濟效益問題仍亟待解決。

針對此問題國內(nèi)外學者開展一系列研究，張帆等以甘肅省黑河中游17 個灌區(qū)間水資源優(yōu)化配置為例，以經(jīng)濟效益、社會效益、生態(tài)效益為目標函數(shù)構建多目標優(yōu)化模型，分別使用傳統(tǒng)方法與復合多目標方法進行求解，所獲優(yōu)化配置方案可將灌溉水利用系數(shù)提高5.42%～7.57%，驗證復合多目標方法所獲得的優(yōu)化方案更能體現(xiàn)決策者對研究區(qū)域種植業(yè)發(fā)展與灌區(qū)水資源配置的多元要求[1]，但單目標模型在其對應目標上表現(xiàn)更優(yōu)。杜麗娟等采用GWAS 模型，通過劃分水資源配置子單元和設置調蓄節(jié)點，采用公平性最優(yōu)和供水缺水率最小為目標函數(shù)，總量控制、供水能力、分質供水等為約束條件，將水資源優(yōu)化配置問題模擬為生物進化問題，采用基于精英策略的非支配遺傳改進算法求解，建立適用于淠史杭灌區(qū)水資源優(yōu)化配置模型[2]，但其研究分析頻率對象較少，難以真實反映情況，出現(xiàn)偏差幾率較大。Cheng等考慮目標函數(shù)各變量在同一階段獨立且約束條件變量逐級遞增的特殊性，提出一種一維動態(tài)規(guī)劃遞推優(yōu)化方法，得到初始解離散差分動態(tài)規(guī)劃（DDDP）算法，解決多維動態(tài)規(guī)劃模型中每個階段包含3 個狀態(tài)變量和3 個決策變量的問題[3]。李茉等針對灌溉水資源優(yōu)化配置中存在的非線性和不確定性等特點，建立區(qū)間線性分式規(guī)劃（ILFP）模型以獲得最大灌溉水生產(chǎn)力和考慮下層農(nóng)民利益的區(qū)間二次規(guī)劃（IQP）模型以獲得最大產(chǎn)量，并在此基礎上構建LFQBP 模型，促進灌區(qū)可持續(xù)發(fā)展[4]。譚倩等針對農(nóng)業(yè)水資源多目標規(guī)劃中存在的權重不確定性難題，建立基于魯棒優(yōu)化方法的農(nóng)業(yè)水資源多目標優(yōu)化配置模型方法（MRPWU），達到節(jié)省灌溉用水，提高生態(tài)效益的目的[5]。綜上表明，在已完成配套和節(jié)水改造的大型灌區(qū)內(nèi)建立合理的優(yōu)化灌溉制度方案是提升灌區(qū)整體經(jīng)濟效益的有效方法。本文針對此問題建立動態(tài)規(guī)劃模型研究灌區(qū)單位面積及單位水量凈效益，以獲取最大經(jīng)濟效益。

1 研究區(qū)域概況

1.1 灌區(qū)概況

泰來灌區(qū)位于泰來農(nóng)場境內(nèi)，灌區(qū)北以大排干末端為界，南至農(nóng)場場界，西部以第四作業(yè)區(qū)田間道為界，東部以嫩江堤防為界。灌區(qū)南北長10.5 km，東西寬8.0 km，水田總控制面積3.75×107m2，設計灌溉面積3.33×107m2。灌區(qū)范圍內(nèi)地面比降為1∶300～1∶3 000，地面高程海拔136～133 m。多年平均氣溫4.2 ℃，最冷1月份平均氣溫-17.6 ℃；極端最低氣溫-35.2 ℃。最熱7月份，平均氣溫23.4 ℃，極端最高氣溫41.6 ℃。無霜期約130 d，多年平均降水量369 mm，6～8月份降水量占年降水量73.6%，多年平均蒸發(fā)量1 798.2 mm，年平均≥10 ℃有效積溫2 500～2 700 ℃，年平均日照2 917.8 h，最大凍土深度1.95 m，本區(qū)屬高溫少雨地區(qū)，一般春季降雨少、風大，最高可達8～12級。

1.2 水資源概況

1.2.1 地表水資源概況

泰來農(nóng)場位于嫩江右岸低平原區(qū)，灌區(qū)屬嫩江流域。本區(qū)多年平均降雨量369 mm，6～8月份降水量占年降水量73.6%，多年平均徑流深118.2 mm。

本區(qū)處于松嫩平原最為干旱的平原地區(qū)，當?shù)氐乇韽搅鳂O少，且年際變化較大，枯水年難以利用。季節(jié)性河流一棵樹溝（大排干）在場區(qū)中部穿過，為主要排水溝，承泄區(qū)為嫩江。流經(jīng)本地區(qū)的嫩江干流是灌區(qū)主要水源地，泰來灌區(qū)提水泵站位于嫩江江橋水文站下游49 km處，區(qū)間無較大支流匯入。

1.2.2 地下水資源概況

泰來灌區(qū)地處大興安嶺東南余坡和松嫩平原西部邊緣過渡地帶，地層巖性均為第四系松散堆積物。

根據(jù)《黑龍江省泰來縣地下水資源調查評價報告》[6]和《黑龍江省齊齊哈爾地區(qū)環(huán)境水文地質調查研究報告》[7]得知，該區(qū)為地下水可廣泛開采的富水區(qū)。表層孔隙潛水埋深為15～28 m，含水層厚度為15～90 m，含水層巖性為砂礫卵石、粉砂、細砂，單井出水量50～100 m3·h-1，單井影響半徑為200～300 m。潛水化學類型以重碳酸鈣、重碳酸鈉為主，鐵離子含量為5 mg·L-1，錳離子含量也較高，地下水礦化度0.16 g·L-1，pH 7.1～7.6，一般超過飲用水標準，但水質符合農(nóng)業(yè)灌溉用水要求。

承壓水層厚1.4～5.0 m，埋深36～48 m，日涌水量543 m3，水量豐富，水質優(yōu)，pH 7.1～7.6，礦化度小于1.0 g·L-1，可用于生活飲用水。

該區(qū)域地下水補給模數(shù)為91 500 m3·a-1·km-2，開采模數(shù)為58 564 m3·a-1·km-2。

本區(qū)地下水資源為當?shù)鼐用裼盟秃堤飮姽嗵峁┏渥闼础?/p>

2 灌區(qū)水資源及水量供需平衡分析

2.1 灌區(qū)水資源分析

2.1.1 一棵樹溝（大排干）水資源分析

一棵樹溝匯水面積1.85×108m2，河道蛇曲發(fā)育，河長61 km，場區(qū)內(nèi)長度為22.2 km，平均寬度30 m。河道雖常年有水，但流量較小，年平均流量0.09 m3·s-1，設計標準下不能作為灌溉水源，主要作用是區(qū)內(nèi)排水。

2.1.2 嫩江水資源分析

嫩江為松花江北源，發(fā)源于大興安嶺伊勒呼里山中段南側。嫩江干流由北向南流經(jīng)嫩江、莫力達瓦、訥河、富裕、甘南、齊齊哈爾、龍江、泰來、杜爾伯特、大安、肇源等市、縣、旗，在吉林省扶余縣三岔河從左岸匯入松花江。嫩江干流全長為1 370 km，流域總面積為2.97×1011m2，泰來灌區(qū)渠首引水口處匯水面積為1.77×1011m2。

嫩江在泰來農(nóng)場東部經(jīng)過，流經(jīng)農(nóng)場段長12.5 km，河道自然比降0.35‰，河寬250～300 m，行洪最小寬度5.5 km，堤防長度10 km，防洪標準為20年一遇。嫩江流域徑流主要集中在6～9月，約占全年徑流量70%，春季4～5月份枯水期徑流量約占全年15%，封凍期徑流占全年10%。

根據(jù)《泰來農(nóng)場水利工程初步設計》，嫩江中下游泰來農(nóng)場段多年平均年徑流量為2.10×1010m3，設計保證率情況下年徑流量4.82×109m3?？菟髁?5.8 m3·s-1。

2004年全國第二次水資源規(guī)劃提出嫩江流域多年平均地表水資源量為2.94×1010m3，其中江橋以上為2.47×1010m3，比原規(guī)劃增加6.66×109m3，江橋以上增加3.69×109m3，江橋以下增加2.97×109m3。

2.1.3 地表水資源開發(fā)利用現(xiàn)狀

嫩江流域地表水資源開發(fā)利用包括黑龍江省、吉林省、內(nèi)蒙古自治區(qū)用水。根據(jù)尼爾基水庫初設統(tǒng)計，嫩江流域干流現(xiàn)狀地表水總用水量約2.25×109m3。其中城市居民生活及工業(yè)用水量6.4×108m3，農(nóng)牧業(yè)灌溉用水量1.23×109m3，濕地生態(tài)用水量3.75×108m3；嫩江流域七大支流現(xiàn)狀總用水量約1.72×109m3，其中居民生活及工業(yè)用水量3×107m3；農(nóng)牧業(yè)灌溉用水量1.58×109m3，濕地生態(tài)用水量1.12×108m3。干支流總用水量約3.97×109m3，占嫩江流域地表水資源量17.5%。江橋以上現(xiàn)狀總用水量約2.41×109m3，占天然徑流11.5%。

綜上所述，本區(qū)水資源量豐富，為當?shù)厮锓N植提供充足資源，灌區(qū)水源有保證。

2.1.4 灌區(qū)現(xiàn)狀水量供需分析

2.1.4.1 灌區(qū)水田現(xiàn)狀水量供需分析

灌區(qū)水源地嫩江水量充足，根據(jù)渠首泵站設計成果，現(xiàn)有嫩江渠首灌溉站設計流量4.86 m3·s-1。

灌區(qū)水田實灌面積3.33×107m2，現(xiàn)狀灌溉水利用系數(shù)為0.65，保證率P=50%情況下水稻凈灌溉定額8 220 m3·hm-1，灌溉毛定額12 645 m3·hm-1，年需水總量為4.215×107m3。如按灌水率1.158×105m3·s-1·hm-1，灌區(qū)最大需水流量5.94 m3·s-1，泵站設計流量4.86 m3·s-1，供水流量無法滿足3.33×107m2水田同時引水灌溉要求。但因2號水庫蓄水容積存蓄泵站引水，按照水泵每年供水時間122 d、每日工作時間20 h 計，泵站年供水能力4.268×107m3，滿足灌區(qū)現(xiàn)狀水田灌溉用水要求。

2.1.4.2 區(qū)內(nèi)其他各業(yè)現(xiàn)狀水量供需分析

根據(jù)對泰來縣水資源開發(fā)利用狀況的分析及建議得知[8]，該區(qū)域地下水補給模數(shù)為9.15×104m3·a-1·km-2，開采模數(shù)為5.8 564×104m3·a-1·km-2。區(qū)域1.02×108m2內(nèi)年可開采量5.92×106m3。

區(qū)內(nèi)其他用水包括生活用水、旱田噴滴灌用水和生態(tài)用水等，均取自地下水。生活用水主要為人畜用水，項目區(qū)內(nèi)現(xiàn)狀總人口0.28 萬人，牲畜1 245 頭，根據(jù)《黑龍江省地方標準〈用水定額〉（DB23/T727-2003）》[9]，農(nóng)村生活人均用水定額0.080 m3·d-1，牲畜每頭平均用水定額0.1 m3·d-1。每年人畜共消耗水量1.3×105m3。目前旱田灌溉面積6.67×106m2，按用水毛定額1 050 m3·hm-1計，年用水7×105m3?，F(xiàn)狀生態(tài)用水量為場區(qū)綠化地用水，體量較小，忽略不計。上述各業(yè)年用水總量為8.3×105m3。區(qū)內(nèi)年可開采量遠大于生活用水和旱田灌溉用水。

綜上，灌區(qū)水利工程現(xiàn)狀供水能力4.351×107m3，各業(yè)需水量4.298×107m3，水量富余5.3×105m3。

3 基于單位面積經(jīng)濟效益最大的動態(tài)規(guī)劃模型及ASGA算法

3.1 模型設定

此動態(tài)規(guī)劃模型實質是在保證灌水量前提下，在單位面積上獲得最大經(jīng)濟效益。通過優(yōu)化灌溉制度，實現(xiàn)灌區(qū)合理灌水。在穩(wěn)步提升經(jīng)濟效益前提下，提高灌區(qū)水資源利用效率。以泰來灌區(qū)水稻灌溉制度為研究對象，以單位面積實際產(chǎn)量Ya與最高產(chǎn)量Ym比值最大為優(yōu)化目標[10]。

3.2 單作物優(yōu)化灌溉制度模型

3.2.1 目標函數(shù)

目標函數(shù)采用在供水不足條件下，水量和農(nóng)作物實際產(chǎn)量連乘模型，目標函數(shù)為單位面積實際產(chǎn)量Ya與最高產(chǎn)量Ym的比值最大[11]，即

式中，ETa為實際蒸發(fā)量（mm）；ETm為最大蒸發(fā)蒸騰量（mm）；n為模型階段總數(shù)；λi為缺水敏感指數(shù),；YD為泰來灌區(qū)水稻多年平均最高產(chǎn)量，YD＝7 800 kg·hm-2；T1為泰來水稻市場平均價格，T1＝3.0 CNY·kg-1；Q為灌水總量（mm）；T2為水價，其中泰來灌區(qū)現(xiàn)狀水價為0.1 CNY·m-3；T3為其他農(nóng)業(yè)費用，包括種子、化肥、農(nóng)藥、人工、電費、物資等，平均8 000 CNY·hm-2。

3.2.2 階段變量

根據(jù)作物生育過程，將作物全生育期劃分為插秧返青期、分蘗期、拔節(jié)期、抽穗期、乳熟期。

3.2.3 狀態(tài)變量

狀態(tài)變量為各階段初可用于分配的灌溉水量qi及初始田面蓄水深度hi。

3.2.4 系統(tǒng)方程

水量分配方程：

式中，qi、qi+1為第i及第i+1階段初系統(tǒng)可用于分配水量（mm）；mi為第i階段灌水量（mm）。

田間水量平衡方程：

式中，hi、hi+1為第i及第i+1 階段初田面水層深度（mm）；Ci為第i階段排水量（m3）；其余符號意義同前。

3.2.5 約束條件

3.2.5.1 決策約束

式中，Q為全生育期單位面積上可供分配水量（mm）；（ETm）i、（ETmin）i為第i階段最大及最小蒸發(fā)蒸騰量（mm）。

3.2.5.2 田面水層hi約束

3.2.6 初始條件

3.2.6.1 初始田面水深

3.2.6.2 第一階段可用于分配水量

第一階段可用于分配水量為水稻全生育期可用于分配水量，即q1=Q1。

3.3 ASGA算法介紹

遺傳算法作為一種有效全局搜索方法，其應用領域已滲透到多個學科中，如函數(shù)優(yōu)化、組合優(yōu)化、生產(chǎn)調度、自動控制、機器學習、圖像處理、人工生命、遺傳編程、機器學習、數(shù)據(jù)挖掘等[12-13]。傳統(tǒng)GA（見圖1）主要通過初始化種群，選擇、交叉和變異算子尋求最優(yōu)解以解決目標函數(shù)優(yōu)化問題[14]，但由于其進化鏈條單一，有其自身缺點，如易產(chǎn)生早熟收斂、收斂速度慢、局部尋優(yōu)能力較差等，且二進制編碼遺傳算法進行數(shù)值優(yōu)化時，存在精度低缺點[15-16]。

圖1 傳統(tǒng)GA流程Fig.1 Flow of traditional GA

ASGA（見圖2）是在傳統(tǒng)GA 中添加鏡像算子（Mirror）、元組（Cell）操作和對比種群淘汰庫：

圖2 ASGA流程Fig.2 ASGA flow

①鏡像算子（Mirror）是一種對二進制編碼染色體改變較強的算法，轉碼變化如下：

式中，X為十進制編碼；L為二進制編碼長度；x為二進制編碼；i為左起二進制位數(shù)，以L=20，隨機生成512組二進制編碼種群進行鏡像算子變化效果試驗，結果如下：

由概率密度圖（見圖3）可知，染色體鏡像前與鏡像后各區(qū)域十進制數(shù)值取值概率分布相差較大，且二者呈現(xiàn)出反向極值現(xiàn)象，即在鏡像算子優(yōu)化前較多出現(xiàn)的十進制數(shù)值在鏡像算子優(yōu)化后出現(xiàn)次數(shù)驟減，因此可看出鏡像算子是將種群內(nèi)個體值向其相反性質方向改變，可使劣質個體在一定程度上向優(yōu)質種群改變，以提升種群個體的最優(yōu)概率。

圖3 鏡像前后概率密度Fig.3 Probability density before and after mirroring

②元組操作是將二進制數(shù)據(jù)切塊，種群分為n節(jié)（n為變量數(shù)），形成n個相對獨立染色體模塊并鏈接，形成染色體鏈（Chromosome-link），然后元組化地進行變異、交叉、鏡像等操作，可節(jié)省代碼量，條理清晰，提高程序運行速度，且各變量間變化維持相對獨立性，在種群中可生成多種混合變量，以避免早熟問題[11]。

初始化種群染色體長度如下：

式中，Lgro為種群染色體總長度；Li為第i染色體模塊長度，n為求解變量數(shù)。

③對比種群淘汰庫。在優(yōu)化過程中對比種群個體，淘汰較差個體，列入黑名單。淘汰庫隨優(yōu)化過程不斷更新，且淘汰庫可用系數(shù)矩陣建立，占用空間小。

3.4 求解參數(shù)

參考充分灌溉作物生育時段劃分方法，結合灌區(qū)多年耐冷優(yōu)質水稻松粳系列品種葉齡模式判斷生育進程，并參考2002年、1998年、2017年水稻試驗結果，P=50%情況下水稻生育期時段參數(shù)見表1，P=75%情況下水稻生育期時段參數(shù)見表2，P=90%情況下水稻生育期時段參數(shù)見表3。

表1 泰來灌區(qū)水稻非充分灌溉下優(yōu)化灌溉制度求解數(shù)據(jù)（P=50%）Table 1 Solution data of optimized irrigation schedule under insufficient rice irrigation in Tailai Irrigation District

表2 泰來灌區(qū)水稻非充分灌溉下優(yōu)化灌溉制度求解數(shù)據(jù)（P=75%）Table 2 Solution data of optimized irrigation schedule under insufficient rice irrigation in Tailai Irrigation District

表3 泰來灌區(qū)水稻非充分灌溉下優(yōu)化灌溉制度求解數(shù)據(jù)（P=90%）Table 3 Solution data of optimized irrigation schedule under insufficient rice irrigation in Tailai Irrigation District

4 求解結果及分析

對不同用水保證率（50%、75%、90%）水資源進行優(yōu)化配置。水稻單作物優(yōu)化應用ASGA優(yōu)化原理，選定父代初始種群規(guī)模為N=400，交叉概率pc=0.80，變異概率pm=0.20，鏡像概率pj=0.70，優(yōu)秀個體數(shù)目選定為20 個，將原始數(shù)據(jù)輸入，加速運行，得出不同供水量下目標函數(shù)最優(yōu)值及優(yōu)化結果見表4～7。保證率P=50%情況下單位面積凈效益隨灌溉定額變化見圖4。

圖4 單位面積（hm2）凈效益隨灌溉定額變化（P=50%）Fig.4 Changes in economic benefits per unit area(hm2)with irrigation volume

表4 泰來灌區(qū)水稻非充分灌溉制度優(yōu)化成果（P=50%）Table 4 Optimization results of insufficient irrigation system for rice in Tailai Irrigation District

表5 泰來灌區(qū)水稻非充分灌溉制度優(yōu)化成果（P=75%）Table 5 Optimization results of insufficient irrigation system for rice in Tailai Irrigation District

表6 泰來灌區(qū)水稻非充分灌溉制度優(yōu)化成果（P=90%）Table 6 Optimization results of insufficient irrigation system for rice in Tailai Irrigation District

表7 基于ASGA的灌區(qū)水資源優(yōu)化配置結果Table 7 Result of optimal allocation of water resources in irrigation area based on ASGA

由圖4可知，當灌區(qū)分配灌溉水量較小時，灌區(qū)水稻凈效益值隨灌水量增加而大幅增加，這是因水稻在極度缺水情況下，每增加單位灌水量即大大提高水稻單位面積凈效益。隨灌水量增加，水稻凈效益曲線坡度不斷變緩，即灌水量對單位面積經(jīng)濟效益影響因子逐漸減小，當灌溉定額為8 220 m3時，單位面積凈效益達到最大值。

保證率P=50%情況下單位灌水量凈效益隨單位灌水量變化見圖5，由此可見，隨灌區(qū)灌水量增加，單位灌水量產(chǎn)生單位面積凈效益值增加梯度逐漸減小，當灌水量達8 000 m3時，單位水量所產(chǎn)生凈效益接近峰值，此時再增加灌水量時，單位水量凈效益值出現(xiàn)下降趨勢。

圖5 單位灌水量經(jīng)濟效益隨單位灌水量變化（P=50%）Fig.5 Economic benefit of unit irrigation volume varies with the unit irrigation volume

5 結 論

a.本文充分考慮灌區(qū)灌溉水量對作物各生育期的影響，建立基于單位面積經(jīng)濟效益最大的動態(tài)規(guī)劃模型，運用ASGA 優(yōu)化算法對目標函數(shù)求解，結果表明，該模型與算法適合該灌區(qū)水資源優(yōu)化配置。對比張帆等以經(jīng)濟效益、社會效益、生態(tài)效益為目標函數(shù)構建的多目標優(yōu)化模型[1]，本模型對經(jīng)濟效益單目標優(yōu)化結果更好，但優(yōu)化目標單一，無法衡量灌區(qū)整體效益。對比杜麗娟等構建的GWAS 模型[2]，本模型優(yōu)化效率更高，但也與其存在類似情況，因分析頻率對象較少，無法真實反映情況，易出現(xiàn)偏差。

b.泰來灌區(qū)水資源量豐富，為當?shù)厮锓N植提供充足資源，灌區(qū)水源有保證。在不同灌溉設計保證率情況下，隨灌水量增加，灌區(qū)單位面積凈效益呈上升趨勢，但提升幅度逐漸減小，說明在灌區(qū)灌水量逐漸增加情況下，供水量對單位面積凈效益的影響因子越來越小。所以通過增加灌水量提高灌區(qū)整體經(jīng)濟效益的方法并不可取，在水資源不足或水價上漲情況下，采用非充分灌溉制度更為合理。

c.本研究采用ASGA 優(yōu)化算法具有一定創(chuàng)新性，解決傳統(tǒng)算子陷入局部最優(yōu)解的問題，且模型計算速度快、優(yōu)化效率高、方案新穎，有較高實用價值。本文也有不足之處，未考慮水價變化因素，未來研究可將水價因素并入優(yōu)化算法中進行灌區(qū)整體水資源優(yōu)化配置，且優(yōu)化目標單一，加入生態(tài)效益、社會效益等目標組合成多目標優(yōu)化算法。