基于模擬優(yōu)化模型的渠井結(jié)合灌區(qū)多目標(biāo)水資源優(yōu)化配置

王芊予，胡天林，芮松楠，褚夢喬，降亞楠，2

（1.西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院，陜西 楊凌 712100；2.西北農(nóng)林科技大學(xué)旱區(qū)農(nóng)業(yè)水土工程教育部重點(diǎn)實(shí)驗(yàn)室，陜西 楊凌 712100）

0 引 言

水是糧食生產(chǎn)和社會經(jīng)濟(jì)發(fā)展不可替代的寶貴資源，伴隨著生活水平的提高和飲食結(jié)構(gòu)的變化，我國對優(yōu)質(zhì)農(nóng)產(chǎn)品的需求不斷增加，導(dǎo)致農(nóng)業(yè)用水量不斷增加。受我國水資源宏觀稀缺的約束，經(jīng)濟(jì)社會的高速發(fā)展和人口數(shù)量的不斷攀升導(dǎo)致工業(yè)用水不斷擠占農(nóng)業(yè)和生態(tài)用水，部分地區(qū)農(nóng)業(yè)用水矛盾突出[1]。尤其在氣候干燥、降水稀少的西北地區(qū)，有限的地表水資源難以滿足灌溉用水的需求，因此常采用井灌來補(bǔ)充灌溉，形成了典型的渠井結(jié)合灌溉模式，實(shí)現(xiàn)了地表水和地下水的聯(lián)合利用，提高了灌溉用水保證率，涵養(yǎng)了寶貴的地下水資源。但近年來由于地區(qū)性降雨較少，上游河流開發(fā)力度加大以及引水工程的老化失修，導(dǎo)致灌溉引水量不斷減少，地表水用水保證率無法得到保障，再加之井灌的靈活性、便捷性和經(jīng)濟(jì)性，部分農(nóng)戶轉(zhuǎn)向采用地下水進(jìn)行灌溉，導(dǎo)致地下水過量開采，渠井用水比例失調(diào)[2]，灌區(qū)渠灌用水下降，地下水補(bǔ)給量減少而無法得到有效涵養(yǎng)，地下水水位持續(xù)下降，形成了地下水降落漏斗[3]，機(jī)井出水量下降、灌溉時(shí)間和抽水能耗激增，甚至機(jī)井報(bào)廢。傳統(tǒng)的灌區(qū)水資源配置沒有與地下水?dāng)?shù)值模擬模型進(jìn)行緊密耦合，因此對灌區(qū)水文地質(zhì)條件、種植結(jié)構(gòu)的空間分布考慮不足，容易導(dǎo)致灌區(qū)地下水開采總量不變，但地下水局部超采的現(xiàn)象，不滿足地下水管理?xiàng)l例[4]提出的取水總量和地下水位雙控指標(biāo)。因此亟需構(gòu)建基于模擬優(yōu)化模型的水資源優(yōu)化配置模型，將地下水?dāng)?shù)值模擬模型與優(yōu)化算法進(jìn)行耦合，充分考慮研究區(qū)的水文地質(zhì)情況、邊界條件，補(bǔ)給、徑流和排泄的時(shí)空分布特性，采用數(shù)值模擬模型對灌溉條件下的地下水的時(shí)空動態(tài)進(jìn)行模擬，最終得到空間分布的模擬結(jié)果，進(jìn)一步提高渠井結(jié)合灌區(qū)水資源配置的合理性，促進(jìn)灌區(qū)灌溉用水方式由粗放低效向節(jié)約集約的精細(xì)化管理轉(zhuǎn)變[5]。

當(dāng)前的水資源優(yōu)化配置?？紤]經(jīng)濟(jì)、生態(tài)、社會等多方面目標(biāo)。邵東國[6]等通過考慮生態(tài)環(huán)境保護(hù)、水權(quán)轉(zhuǎn)讓、利益補(bǔ)償?shù)纫蛩兀跐M足人們生存、工程供水能力，用水公平性的基礎(chǔ)上，構(gòu)建了水資源凈效益最大的目標(biāo)函數(shù)，并采用水資源系統(tǒng)的熵變關(guān)系對水資源配置進(jìn)行合理性評價(jià)。付銀環(huán)[7]等基于區(qū)間兩階段隨機(jī)規(guī)劃的方法，以灌區(qū)用水成本最小為目標(biāo)結(jié)合作物水分生產(chǎn)函數(shù)，建立了不確定性的水資源優(yōu)化配置模型。粟曉玲[8]等以生態(tài)、經(jīng)濟(jì)綜合效益最大為目標(biāo)建立單元種植結(jié)構(gòu)優(yōu)化模型，再利用水資源轉(zhuǎn)化模擬模型對結(jié)果進(jìn)行求解，實(shí)現(xiàn)了先優(yōu)化后模擬的模擬優(yōu)化松散耦合模型。Farhadi[9]等以灌溉缺水量最少，水量分配公平度最大和地下水下降最小為目標(biāo)，將MODFLOW 求解出的分配方案作為數(shù)據(jù)訓(xùn)練集，進(jìn)而構(gòu)建基于神經(jīng)網(wǎng)絡(luò)的多目標(biāo)優(yōu)化模型并結(jié)合納什議價(jià)模型進(jìn)行方案選取。綜上，在目標(biāo)的量化方面，目前研究大多通過貨幣價(jià)值衡量水資源價(jià)值，而忽略其社會屬性、環(huán)境屬性和公有性[10]，當(dāng)前盡管有部分模型在選取配置方案的過程中考慮了行政區(qū)間配水的公平性問題，但針對渠井結(jié)合灌區(qū)考慮用水公平度的模型較少，在實(shí)際應(yīng)用中不利于灌區(qū)的統(tǒng)觀統(tǒng)管，在我國種植結(jié)構(gòu)受政府政策、市場和政府引導(dǎo)的大背景下，在水資源配置中需要充分考慮特色名優(yōu)經(jīng)濟(jì)作物的用水特點(diǎn)，也應(yīng)保證糧油等作物的用水需求，進(jìn)一步提升水利對地方特色經(jīng)濟(jì)的推動作用，因此有必要在水資源配置中進(jìn)一步考慮各區(qū)域的配水公平性。

隨著計(jì)算機(jī)領(lǐng)域的優(yōu)化算法和水資源模型理論的高速發(fā)展，水資源模擬優(yōu)化模型的研究日益深入。譚倩[11]等人將魯棒規(guī)劃的多目標(biāo)方法引入水資源優(yōu)化研究中，建立基于魯棒規(guī)劃方法的農(nóng)業(yè)水資源多目標(biāo)優(yōu)化配置模型（MRPWU 模型），通過引入保護(hù)函數(shù)和非線性保護(hù)函數(shù)線性化的方式有效處理雙目標(biāo)規(guī)劃中權(quán)重的不確定性。邵東國[12]等通過構(gòu)建AquaCrop 作物模型，對不同作物的產(chǎn)量與灌排關(guān)系進(jìn)行模擬計(jì)算，以水稻產(chǎn)量、灌水量、排水量為優(yōu)化目標(biāo)，提出基于AquaCrop 模型與熵值法耦合的多目標(biāo)多情景優(yōu)化方法。孫月峰[13]等人采用混合遺傳模擬退火算法建立了以經(jīng)濟(jì)、社會和環(huán)境的綜合效益最優(yōu)為目標(biāo)的優(yōu)化配置模型，有效提高了求解速度和解的精度。楊蘊(yùn)[14]等人將基于過渡帶理論的海水入侵變密度流數(shù)值模擬程序SEAWAT 同遺傳混合算法NPTSGA 相耦合，開發(fā)了海水入侵條件下地下水多目標(biāo)模擬優(yōu)化管理模型，以實(shí)現(xiàn)多重管理的目標(biāo)。傳統(tǒng)水資源多目標(biāo)模擬模型的優(yōu)化方法采用目標(biāo)加權(quán)法、距離函數(shù)法或最小最大準(zhǔn)則將多目標(biāo)問題轉(zhuǎn)化為單目標(biāo)問題進(jìn)行求解[15]，但由于權(quán)重和需求水平的高敏感性使得解集的優(yōu)劣性無法得到保證。而現(xiàn)階段的優(yōu)化方法也可以耦合數(shù)值模擬模型進(jìn)行求解，但多采用松散耦合的方式，雖然提高了結(jié)果的合理性和可行性，但采用的耦合方式使得計(jì)算機(jī)內(nèi)存在優(yōu)化算法和模擬模型之間無法快速的共享數(shù)據(jù)，極大地降低了計(jì)算速度和求解效率[16]。緊密耦合的模擬優(yōu)化模型可采用統(tǒng)一的編程語言把優(yōu)化算法與數(shù)值模擬模型進(jìn)行緊密耦合，將數(shù)值模擬模型計(jì)算的分布式結(jié)果作為目標(biāo)函數(shù)值，極大地增強(qiáng)了目標(biāo)函數(shù)對空間信息的提取能力，也可將空間分布的約束條件施加到相應(yīng)的位置（單元格）上，在空間上對地下水的開采進(jìn)行調(diào)控。目前研究中采用的多目標(biāo)優(yōu)化算法如蟻群算法[17]、粒子群算法[18]、NSGA-II[19]等啟發(fā)式算法，在處理多個(gè)目標(biāo)（5 個(gè)及以上，Many-Objective Optimization）函數(shù)時(shí)搜尋能力較低，采用專門針對多目標(biāo)改進(jìn)的NSGA-III算法的研究較少。

基于此本研究從一個(gè)假想算例著手，基于Python 語言將優(yōu)化框架Pymoo 中的NSGA-Ⅲ算法與地下水?dāng)?shù)值模擬軟件包FloPy 進(jìn)行緊密耦合，通過編寫模塊化目標(biāo)函數(shù)建立典型渠井結(jié)合灌區(qū)的多目標(biāo)水資源配置模擬優(yōu)化模型，以探索經(jīng)濟(jì)效益、區(qū)域配水公平、地下水可持續(xù)利用等多個(gè)目標(biāo)對配置方案的影響，同時(shí)給決策者提供基于不同偏好和基于博爾達(dá)計(jì)數(shù)規(guī)則確定的兩類水資源配置方案，為灌區(qū)地表水地下水可持續(xù)安全高效利用，地下水取水總量和地下水水位雙控提供強(qiáng)有力的技術(shù)支撐。

1 研究方法

1.1 地下水?dāng)?shù)值模擬模型

MODFLOW 是由美國地質(zhì)調(diào)查局采用有限差分法開發(fā)的用于孔隙介質(zhì)中三維地下水流模擬的模塊化數(shù)值模擬模型軟件[20]，由一個(gè)主程序與若干相對獨(dú)立的子程序包組成，通過定義子程序（補(bǔ)給、河流、井、排水等）中的相關(guān)模塊進(jìn)行數(shù)值模擬過程。鑒于其顯著的模塊化結(jié)構(gòu)，Bakker[21]等人基于Python編寫了地下水?dāng)?shù)值模擬FloPy庫。

由于灌區(qū)承壓水開發(fā)利用程度較低，資料匱乏，目前灌溉主要開發(fā)利用潛水，且淺層地下水以垂向運(yùn)動為主，故將研究區(qū)潛水水流系統(tǒng)概化為均質(zhì)各向同性二維非穩(wěn)定流單層潛水含水層，地下水流數(shù)學(xué)模型[22]如下：

式中：K為潛水含水層滲透系數(shù)，m/d；H為潛水水位，m；Z為含水層底板高程，m；W為源匯項(xiàng)，m/d；μ為潛水含水層給水度；t為時(shí)間，d；Ω為計(jì)算區(qū)域；x,y為坐標(biāo)；Γ1，Γ2為第一、二類邊界；q(x,y,t)為流量邊界單寬流量，m2/d。

本研究通過調(diào)用FloPy 中相應(yīng)的子程序包基于上述數(shù)學(xué)模型定義并修改含水層網(wǎng)格、水文地質(zhì)參數(shù)、初始條件、邊界條件等方面，通過自定義灌區(qū)補(bǔ)給項(xiàng)（包括降雨、田灌入滲、井回歸補(bǔ)給等）和灌區(qū)排泄項(xiàng)，開發(fā)、運(yùn)行精細(xì)化模型并讀取分析結(jié)果。

1.2 多目標(biāo)遺傳算法NSGA-Ⅲ

多目標(biāo)問題的傳統(tǒng)求解方法是通過加權(quán)或約束折衷等方法將多目標(biāo)優(yōu)化問題轉(zhuǎn)化為單目標(biāo)優(yōu)化問題，當(dāng)需獲得多個(gè)不同解決方案時(shí)只能通過重復(fù)應(yīng)用模擬來實(shí)現(xiàn)，而進(jìn)化算法能夠在一次模擬運(yùn)行的過程中找到一系列帕累托最優(yōu)解集，非支配排序遺傳算法（NSGA）則是最早的進(jìn)化算法之一。由于在排序過程中出現(xiàn)計(jì)算復(fù)雜度較高、易丟失良好解決方案和需提前設(shè)定共享參數(shù)的問題，Deb[23]等人提出改進(jìn)的非支配排序遺傳算法（NSGA-Ⅱ），通過快速非支配排序的方法降低計(jì)算復(fù)雜度，引入精英策略以增大采樣空間，提出擁擠距離代替共享參數(shù)，以保證種群在進(jìn)化過程中的多樣性。隨著優(yōu)化目標(biāo)數(shù)量的進(jìn)一步增多，求解出的非支配個(gè)體在種群中的占比呈指數(shù)型上升，非支配解集的收斂性顯著下降，為保證多樣性而引入的擁擠距離算子計(jì)算代價(jià)增加，為解決以上問題，Deb 和Jain[24,25]進(jìn)一步提出基于參考點(diǎn)的第三代非支配排序遺傳算法（NSGA-Ⅲ），采用邊界交叉構(gòu)造權(quán)重的方法產(chǎn)生參考點(diǎn)，在臨界層的環(huán)境選擇中通過預(yù)定義多個(gè)分布良好的參考點(diǎn)以維持種群的多樣性，修改精英選擇算子和后代種群以求解一般約束的多目標(biāo)優(yōu)化問題，有效解決了高維目標(biāo)優(yōu)化問題中的解集收斂性問題。

2 模擬優(yōu)化模型的構(gòu)建

2.1 數(shù)值模擬模型的構(gòu)建

康燕楠[26]等在對西北典型渠井結(jié)合灌區(qū)寶雞峽灌區(qū)調(diào)研的基礎(chǔ)上，綜合灌區(qū)的實(shí)際情況設(shè)計(jì)了一個(gè)灌溉設(shè)施完善包含不同水文地質(zhì)分區(qū)的典型灌域，該灌域僅采用渠灌或井灌均可滿足灌溉要求，長5 km，寬3.5 km，在空間上剖分成50 行50 列的長方形網(wǎng)格，東西邊界概化為定水頭邊界，南北邊界概化為零流量邊界，灌域內(nèi)共設(shè)置50 眼抽水井，具體見文獻(xiàn)[26]。本文在此基礎(chǔ)上，根據(jù)不同作物種類將灌區(qū)在空間上劃分為5 個(gè)種植區(qū)域，同時(shí)基于行政管理的考慮將灌域劃分為5個(gè)行政分區(qū)，通過行政分區(qū)1、5 和2、3、4 不同的地下水埋深來模擬灌區(qū)塬上和塬下兩個(gè)灌溉系統(tǒng)，在各行政分區(qū)設(shè)置一口觀測井以便監(jiān)測地下水位，渠井結(jié)合灌域的模型示意圖如圖1所示，相關(guān)水文地質(zhì)參數(shù)參考鄧康婕[27]根據(jù)涇惠渠灌區(qū)成果報(bào)告等資料確定，種植結(jié)構(gòu)與灌溉周期[28-30]如表1和表2所示。

表1 水文地質(zhì)參數(shù)表Tab.1 Hydrogeological parameters of the case study

表2 種植結(jié)構(gòu)與灌溉制度表Tab.2 Crop patterns and irrigation schedule

圖1 典型灌域示意圖（單位：km）Fig.1 A diagram of geological model of typical irrigation area

本文使用FloPy 構(gòu)建MODFLOW 模型，首先導(dǎo)入FloPy 并創(chuàng)建一個(gè)MODFLOW 的模型對象；利用FloPy 中的Dis 程序包創(chuàng)建含水層，設(shè)定應(yīng)力期和計(jì)算時(shí)間步長為默認(rèn)值1d；采用FloPy 中的Bas 程序包模擬活動單元格的屬性并設(shè)定初始水頭值、Lpf 程序包對含水層的水力特性與類型進(jìn)行定義、Ghb 程序包對含水層的邊界條件進(jìn)行設(shè)定、Rch 程序包將降雨入滲、田灌入滲和井回歸補(bǔ)給給地下水、Wel程序包通過輸入的抽水井位置、抽水速度和應(yīng)力期對含水層的井開采情況進(jìn)行模擬；最后調(diào)用Sip 求解器進(jìn)行求解并使用Oc 程序包對MODFLOW模型的輸出頻率和類型進(jìn)行定義。

2.2 優(yōu)化模型

2.2.1 決策變量

本文針對五種作物分別考慮其灌溉用水量與渠井用水比例，因此共設(shè)置10 個(gè)決策變量，分別為PCWi、Fi，其中：PCWi為各作物用于灌溉的渠井用水比例；Fi為分配給各作物的灌溉用水量，m3；i為作物種類數(shù)量。以各作物在不同灌季不同灌水定額的比例為依據(jù)，將各作物的渠灌用水量與井灌用水量在不同灌季進(jìn)行分配，根據(jù)作物種植結(jié)構(gòu)與行政分區(qū)進(jìn)行空間劃分，以模擬灌溉水入滲補(bǔ)給，包括渠灌回歸、井回歸補(bǔ)給和地下水灌溉開采情況。

2.2.2 目標(biāo)函數(shù)

灌區(qū)水資源優(yōu)化配置是一個(gè)涉及多方利益主體的復(fù)雜問題，因此本文在選擇優(yōu)化目標(biāo)時(shí)，通過引入灌區(qū)凈經(jīng)濟(jì)效益最大目標(biāo)以滿足農(nóng)戶需求，設(shè)置地下水平均累計(jì)下降變化最小目標(biāo)以滿足生態(tài)和地下水管理的可持續(xù)發(fā)展要求，考慮缺水量最小目標(biāo)以滿足作物生長發(fā)育要求，設(shè)置各行政區(qū)配水公平度和用水效率最大目標(biāo)以滿足管理者實(shí)行配水方案的需求。

（1）經(jīng)濟(jì)效益目標(biāo)。灌區(qū)內(nèi)用水產(chǎn)生凈效益最大，即：

式中：Z1為灌區(qū)用水總凈效益，元；i,j分別代表作物類型編號與行政分區(qū)編號；Fij為各作物在各行政區(qū)的灌溉用水量，m3；PCWij為各作物在各行政區(qū)用于灌溉的渠井用水比例；Pij為各作物在各行政區(qū)的單價(jià)，即蔬菜、夏玉米、冬小麥、棉花、大豆的單價(jià)分別為3.9、3.4、3.2、7.0、2.4 元/kg；Aij為各作物在各行政區(qū)的面積，hm2；Pc為渠道引水灌溉的價(jià)格，0.27 元/m3；Pws和Pwx分別為塬上塬下抽水灌溉的價(jià)格0.16 和0.11 元/m3；Si為不同作物產(chǎn)量與灌溉量的具體關(guān)系。

通過查閱作物耗水量與產(chǎn)量相關(guān)關(guān)系文獻(xiàn)[31-34]，得出如下公式：

式中：S1、S2、S3、S4、S5分別為蔬菜、夏玉米、冬小麥、棉花以及大豆的產(chǎn)量函數(shù)，kg/hm2；Q1、Q2、Q3、Q4、Q5分別為蔬菜、夏玉米、冬小麥、棉花以及大豆的全生育內(nèi)的灌溉水量，m3/hm2。

（2）單元格地下水位平均累計(jì)下降變化目標(biāo)。灌區(qū)單個(gè)網(wǎng)格每日地下水位較前一日的下降年累計(jì)值最?。ú豢紤]地下水位上升的情況），即：

式中：Z2為單個(gè)網(wǎng)格的年累計(jì)下降值，m；N為灌區(qū)在空間內(nèi)劃分的總網(wǎng)格數(shù)；hd為灌區(qū)內(nèi)各網(wǎng)格每日地下水位之和，m。

（3）缺水目標(biāo)。各作物需水量與灌溉用水量之差最小，即：

式中：Z3為各作物的缺水量之和，m3；Qi為灌區(qū)各作物在概化模型中的最大灌溉需水量，m3。

（4）公平度目標(biāo)。塬上與塬下的各行政分區(qū)灌溉用水量與需水量比值差異最小，即：

式中：Z4為灌區(qū)塬上塬下公平度之和，公平度目標(biāo)值越小表明越公平；m,n為塬上、塬下行政區(qū)數(shù)；factorj為各行政區(qū)的公平度因子；factorum為塬上各行政區(qū)公平度因子的均值；factordm為塬下各行政區(qū)公平度因子的均值。

（5）用水效率目標(biāo)。各行政區(qū)內(nèi)不同作物總產(chǎn)量與總分配水的比值最大，即：

式中：Z5是灌區(qū)總用水效率之和；Si為各作物產(chǎn)量與需水量的函數(shù)關(guān)系。

2.2.3 約束條件

為保持地下水采補(bǔ)平衡以及灌區(qū)用水特點(diǎn)，考慮約束如下：

（1）供水能力約束：

式中：Wi為用于灌溉各作物的地下水開采量，m3；Ci為用于灌溉各作物的渠灌量，m3；Wmax,i為種植各作物行政區(qū)內(nèi)的所有井的開采能力之和，m3；Cmax,i為種植各作物行政區(qū)內(nèi)的渠道輸水能力之和，m3。

令布爾變量μ(pi,kα)表示協(xié)同成員pi與知識點(diǎn)kα間是否存在關(guān)聯(lián)關(guān)系。協(xié)同成員pi與知識點(diǎn)kα間存在關(guān)聯(lián)關(guān)系，即協(xié)同成員pi掌握知識點(diǎn)kα，則μ(pi,kα)=1；反之，則有μ(pi,kα)=0。因此，K-K子網(wǎng)絡(luò)與P-P子網(wǎng)絡(luò)之間的映射關(guān)系可以表示為

（2）需水量約束：

式中：Fmin,i為各作物所分配用于灌溉的最小水量值，m3，本文取枯水年灌溉定額的60%，F(xiàn)max,i為各作物所分配用于灌溉的最大水量值，m3，本文取枯水年的灌溉定額。

（3）渠井用水比例約束：

式中：PCWmax為使寶雞峽灌區(qū)地下水位抬升兩米的最大渠井用水比例；PCWmin為使寶雞峽灌區(qū)地下水位降低兩米的最小渠井用水比例，參考寶雞峽灌區(qū)渠井用水比例對地下水位影響的相關(guān)研究[35]，取PCWmax為1.97，PCWmin為0（即全井灌）。

（4）非負(fù)約束：

式中：Wi為用于灌溉各作物的井灌量，m3；Ci為用于灌溉各作物的渠灌量，m3。

2.3 模擬優(yōu)化模型

優(yōu)化問題中灌區(qū)總經(jīng)濟(jì)效益目標(biāo)通過調(diào)用灌溉水量—產(chǎn)量函數(shù)對個(gè)體進(jìn)行遍歷計(jì)算得出對應(yīng)目標(biāo)結(jié)果；灌區(qū)地下水位平均累計(jì)下降變化目標(biāo)需調(diào)用2.1 中基于FloPy 構(gòu)建的MODFLOW 地下水?dāng)?shù)值模擬模型對其輸出的水頭數(shù)據(jù)進(jìn)行處理得到所需目標(biāo)結(jié)果；缺水量目標(biāo)則可通過決策變量間的簡單計(jì)算得出相應(yīng)目標(biāo)結(jié)果；公平度目標(biāo)可根據(jù)塬上塬下不同行政分區(qū)的基本數(shù)組運(yùn)算得出；用水效率目標(biāo)可基于各行政區(qū)的經(jīng)濟(jì)效益與用水量進(jìn)行數(shù)組運(yùn)算求解，采用外部耦合的方式，將各優(yōu)化目標(biāo)以目標(biāo)函數(shù)的形式內(nèi)嵌至優(yōu)化模型中，通過決策變量在模擬模型和優(yōu)化模型間進(jìn)行參數(shù)傳遞，不斷重復(fù)調(diào)用優(yōu)化問題中的目標(biāo)函數(shù)和約束條件，使用NSGA-Ⅲ算法進(jìn)行優(yōu)化求解。優(yōu)化問題部分將模擬模型內(nèi)嵌進(jìn)目標(biāo)函數(shù)，優(yōu)化求解則通過不斷調(diào)用優(yōu)化問題計(jì)算目標(biāo)函數(shù)值和約束條件實(shí)現(xiàn)模擬模型和優(yōu)化模型的緊密耦合，具體耦合過程如圖2所示。

圖2 模擬優(yōu)化緊密耦合流程圖Fig.2 Tightly coupled simulation optimization flowchart

3 結(jié)果與分析

3.1 參數(shù)設(shè)置與模型評價(jià)

NSGA-Ⅲ參數(shù)設(shè)置：參考方向采用Das-Dennis 方法，分區(qū)數(shù)partitions=5；選擇空間隨機(jī)采樣、隨機(jī)選擇方式；取用模擬二進(jìn)制交叉和多項(xiàng)式突變，其中交叉率Pc= 0.73，變異率Pm= 0.05；種群大小popsize=130，優(yōu)化搜索代數(shù)generation=100。

將采用模擬優(yōu)化模型求解出的Pareto 解集使用Hypervolume 指標(biāo)[36]（表示解集中的個(gè)體與參考點(diǎn)在目標(biāo)空間中所圍成的超立方體體積）進(jìn)行評價(jià)，該指標(biāo)可以有效評價(jià)解集的收斂性、均勻性以及廣泛性，同時(shí)當(dāng)Hypervolume 指標(biāo)達(dá)到最大時(shí)表明該解集收斂至Pareto 最優(yōu)[37]。經(jīng)過反復(fù)調(diào)試，最終選取參考點(diǎn)（40.34，-14 045 971.15，5 777 949.98，1.019，-0.05）并計(jì)算得出相應(yīng)的Hypervolume 指 標(biāo)。Hypervolume指標(biāo)見圖3，算法在50代左右達(dá)到收斂。

圖3 Hypervolume指標(biāo)Fig.3 Hypervolume indicator

3.2 基于不同偏好的配置方案選擇

由于模擬優(yōu)化模型各優(yōu)化目標(biāo)間的排斥性，無法得出唯一的全局最優(yōu)解，而是得到一系列的非劣解集。求解得出的非劣解共57 組，五目標(biāo)對應(yīng)Pareto 非劣解集見圖4（其中散點(diǎn)尺寸越大代表該目標(biāo)方案缺水量越大，顏色接近紅色代表該目標(biāo)方案經(jīng)濟(jì)效益越高，顏色越接近藍(lán)色代表該目標(biāo)方案經(jīng)濟(jì)效益越低）。

圖4 配置方案Pareto解集空間分布圖Fig.4 Spatial distribution of Pareto front

為方便決策者進(jìn)行方案選擇，本研究分別選取灌區(qū)經(jīng)濟(jì)效益最高、地下水位平均累計(jì)降深最小、缺水量最小、公平度最大和用水效率最高5個(gè)不同偏好的典型方案進(jìn)行分析，各不同偏好方案見圖5（圖中各目標(biāo)函數(shù)值經(jīng)極差標(biāo)準(zhǔn)化處理，為避免標(biāo)準(zhǔn)化后目標(biāo)函數(shù)值出現(xiàn)0值，故對標(biāo)準(zhǔn)化后各目標(biāo)值加0.5），具體目標(biāo)值如表3所示。

表3 5個(gè)不同偏好典型方案的目標(biāo)值Tab.3 Typical plan values of five different allocation plans

圖5 各不同偏好方案柱狀圖Fig.5 Objectives values of different allocation plans

在5 個(gè)不同偏好的典型方案中，方案3 經(jīng)濟(jì)效益最高，但缺水量相較均值高11.27%，不利于灌區(qū)作物的生長發(fā)育；方案25 單網(wǎng)格地下水位平均累計(jì)下降變化最小，對灌區(qū)地下水的安全利用最有利；方案15 缺水量最少，但經(jīng)濟(jì)效益相較于均值而言低0.61%，公平度相較于均值而言低16.42%，用水效率相較于均值低26.76%，不利于灌區(qū)配水方案的實(shí)施和經(jīng)濟(jì)的可持續(xù)發(fā)展；方案49 公平度最高，但單網(wǎng)格地下水位平均累計(jì)下降相較于均值而言高10.76%，用水效率相較于均值低19.72%，雖對灌區(qū)配水方案的實(shí)施最有利，但不利于灌區(qū)整體的高效發(fā)展；方案30 用水效率最高，但缺水量相較均值而言高11.06%，影響灌區(qū)作物的產(chǎn)量。

綜上所述，以上5個(gè)目標(biāo)間存在相互排斥關(guān)系，某一目標(biāo)最優(yōu)勢必要犧牲其余目標(biāo)值，決策者可以根據(jù)自己的偏好目標(biāo)進(jìn)行方案選擇，若需要提高經(jīng)濟(jì)效益則選擇方案3，若偏好保護(hù)地下水則選擇方案25（雖渠井用水比例中井灌用水占比較多，但實(shí)際灌溉用水量較少，即對地下水抽取總量較少），而若側(cè)重灌區(qū)作物的生長發(fā)育情況則選擇方案15，若想要更好推進(jìn)配水方案的實(shí)施則選擇方案49，若提倡用水的高效性則選擇方案30，各方案對應(yīng)的灌溉用水量、渠井用水比例見表4。

表4 5個(gè)不同偏好典型方案的灌溉用水量和渠井用水比例Tab.4 Water distribution and water used proportion of five different allocation plans

3.3 基于博爾達(dá)計(jì)數(shù)規(guī)則的配置方案選擇

灌區(qū)水資源配置問題是一個(gè)涉及多因素影響的復(fù)雜過程，單純追求某一因素最優(yōu)而忽略其他方面是不夠合理的，應(yīng)充分考慮各目標(biāo)的實(shí)際情況，綜合選取決策方案。故本文基于博爾達(dá)計(jì)數(shù)規(guī)則[38]，針對不同目標(biāo)對配水方案進(jìn)行打分賦值，分?jǐn)?shù)從1 到m，對應(yīng)目標(biāo)函數(shù)值最優(yōu)的方案得到m分，目標(biāo)函數(shù)值最劣的方案得到1 分，最終將各方案的目標(biāo)分?jǐn)?shù)值累加，取得分最高的方案作為最優(yōu)方案。該計(jì)數(shù)規(guī)則更多地考慮各目標(biāo)的綜合影響，避免對不同目標(biāo)不同偏好引起的投票悖論，充分考慮各目標(biāo)的偏好信息，具有一致性和共識性，各方案總得分情況如圖6所示。

圖6 各方案得分Fig.6 Final score of each pareto front

最終選擇得分較高的前三名作為最優(yōu)方案（見表5和圖7，其中圖7中各數(shù)據(jù)經(jīng)極差標(biāo)準(zhǔn)化處理后加0.5），其中方案31的經(jīng)濟(jì)效益增加0.14%，單元格地下水平均累計(jì)下降變化低4.49%，用水效率高16.9%，缺水量低0.15%，但公平度低52.2%；方案42 的單元格地下水平均累計(jì)下降變化低3.35%，公平度高8.96%，用水效率高15.49%，但經(jīng)濟(jì)效益低0.54%，缺水量高5.83%；方案35 的單元格地下水平均累計(jì)下降變化低4.32%，公平度高1.49%，作物缺水量低5.21%，但用水效率低14.08%，經(jīng)濟(jì)效益低0.08%（以上變幅均基于均值進(jìn)行比較）。

表5 基于博爾達(dá)計(jì)數(shù)規(guī)則的最優(yōu)方案Tab.5 Best solutions chosen acquired based on Borda count

圖7 3個(gè)最優(yōu)方案對比Fig.7 Comparison of three optimal solutions

取博爾達(dá)最優(yōu)水資源配置方案31 的地下水位變化趨勢（圖8）進(jìn)行分析，各個(gè)行政分區(qū)的地下水位隨降水和灌溉開采變化明顯，分區(qū)1年地下水位變幅不超過0.5 m，分區(qū)2年地下水位變幅不超過2 m，分區(qū)3年地下水位變幅不超過0.5 m，分區(qū)4年地下水位變幅不超過1 m，分區(qū)5年地下水位變幅不超過0.5 m。

圖8 博爾達(dá)最優(yōu)配置方案各分區(qū)觀測井地下水位變化趨勢Fig.8 The changing trend of groundwater level for each observed wells based on Borda count

綜上所述，若決策者對經(jīng)濟(jì)效益、地下水可持續(xù)發(fā)展、用水效率以及缺水量四方面較為重視，而對配水公平度要求程度較低，則選取方案31；若決策者傾向配水公平度較高而對作物缺水量要求程度較低，則選取方案42；若決策者對用水效率要求程度較低，可選取方案35，基于博爾達(dá)計(jì)數(shù)規(guī)則選取的各方案對應(yīng)的渠灌、井灌用水量見表6。

表6 基于博爾達(dá)計(jì)數(shù)規(guī)則選取方案的渠灌和井灌用水量m3/hm2Tab.6 The amount of canal and well irrigation water based on Borda count selected solutions

4 結(jié)論與討論

本文基于Python 語言構(gòu)建了一個(gè)緊密耦合的渠井結(jié)合灌區(qū)多目標(biāo)水資源配置模擬優(yōu)化模型，充分考慮降雨入滲補(bǔ)給，地表水灌溉下滲補(bǔ)給地下水，井灌開采抽取地下水以及井灌回歸補(bǔ)給地下水等水循環(huán)環(huán)節(jié)，同時(shí)結(jié)合經(jīng)濟(jì)效益、地下水水位平均下降變化、農(nóng)作物缺水程度、各行政區(qū)域配水公平度和用水效率5 個(gè)目標(biāo)，并采用NSGA-Ⅲ算法進(jìn)行了優(yōu)化求解。

根據(jù)優(yōu)化求解的結(jié)果，本文提供傾向于不同目標(biāo)的配水方案，若偏好全年的經(jīng)濟(jì)效益最高，灌區(qū)蔬菜、玉米、小麥、棉花、大豆渠井用水比例分別為1∶99、2∶98、0、0、2∶98；若傾向維持地下水的可持續(xù)發(fā)展，灌區(qū)五種作物的渠井用水比例分別為2∶98、0、34∶66、1∶99、4∶96；若要保證作物的生長發(fā)育，灌區(qū)五種作物的渠井用水比例分別為1∶98、1∶99、24∶76、37∶63、1∶99；若要強(qiáng)調(diào)各行政區(qū)配水量的公平性，灌區(qū)五種作物的渠井用水比例分別為2∶98、1∶99、0、1∶99、1∶99；若要關(guān)注行政區(qū)用水的高效性，灌區(qū)5 種作物的渠井用水比例分別為0、2∶98、27∶73、1∶99、5∶95。

同時(shí)本文也提供基于博爾達(dá)計(jì)數(shù)規(guī)則綜合考量的配水方案，最終選擇3 個(gè)最優(yōu)方案，其中博爾達(dá)最優(yōu)方案對應(yīng)的蔬菜、玉米、小麥、棉花、大豆的渠灌、井灌用水量分別為19.5、1 930.35；68.34、1298.41 m3/hm2；641.49、723.39；83.91、2 013.87 m3/hm2；533.08、554.84 m3/hm2,決策者可以根據(jù)當(dāng)?shù)氐膶?shí)際情況進(jìn)行決策。

本文通過考慮行政區(qū)劃的配水公平度和用水效率模擬決策者的配水過程，使得模型的模擬結(jié)果更易于實(shí)際應(yīng)用，通過調(diào)整渠井用水比例和灌溉水量來進(jìn)行水資源優(yōu)化配置，為灌區(qū)水資源配置提供了新的參考依據(jù)。但本模型僅基于一個(gè)假想案例，未采用連續(xù)多年資料進(jìn)行多年調(diào)節(jié)，在考慮作物產(chǎn)量—耗水量關(guān)系時(shí)也并未采用目前通用的AquaCrop 作物模型[39]進(jìn)行緊密耦合，同時(shí)并未分區(qū)約束地下水水位以及考慮地下水補(bǔ)償價(jià)格的影響，因此在后續(xù)的研究中可進(jìn)一步將上述不足納入考量并基于真實(shí)灌區(qū)構(gòu)建模型。此外也可以采用不同優(yōu)化算法進(jìn)行求解，橫向?qū)Ρ仍诟髂繕?biāo)情況下不同算法的性能情況和求解結(jié)果的優(yōu)劣性，以實(shí)現(xiàn)配水方案更為精準(zhǔn)的管理。