康燕楠，降亞楠,2，蘇振輝

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

我國是農(nóng)業(yè)大國，水資源與農(nóng)業(yè)生產(chǎn)和糧食安全密切相關(guān)。在氣候變化和人類活動的雙重影響下，我國部分灌區(qū)出現(xiàn)了水資源短缺和生態(tài)環(huán)境退化的現(xiàn)象。近年來，位于北方干旱半干旱地區(qū)的大型灌區(qū)，地表水源難以滿足農(nóng)業(yè)發(fā)展以及生產(chǎn)生活的需求，地下水成為主要的灌溉水源，據(jù)統(tǒng)計大約有70%的地下水用于農(nóng)業(yè)灌溉[1]。由此導(dǎo)致部分灌區(qū)的地下水處于超采狀態(tài)，尤其在干旱缺水年份，來源于降水的天然補(bǔ)給量減少，河流和渠道用水不足時，灌區(qū)會大量開采地下水。不合理的地下水開發(fā)利用容易導(dǎo)致采補(bǔ)失衡，致使部分區(qū)域地下水位持續(xù)下降，形成大面積的水位降落漏斗，進(jìn)而引發(fā)了一系列生態(tài)環(huán)境問題，制約了灌區(qū)的可持續(xù)發(fā)展。因此亟需考慮灌溉效益和合理的地下水位上下限，對灌區(qū)內(nèi)有限的水資源進(jìn)行合理配置，使其發(fā)揮最大效益[2-3]。水資源作為灌區(qū)資源和環(huán)境的雙重要素，對農(nóng)業(yè)生產(chǎn)和生態(tài)環(huán)境保護(hù)起到核心作用，對以地表水和地下水為主的農(nóng)業(yè)水資源進(jìn)行合理配置是解決灌區(qū)缺水問題的重要手段。

灌區(qū)水資源合理配置需要綜合考慮經(jīng)濟(jì)、環(huán)境及社會等多方面的目標(biāo)，是一個典型的多目標(biāo)問題。針對地表水和地下水聯(lián)合利用的水資源合理配置，目前國內(nèi)外的研究趨勢是將優(yōu)化算法和地下水?dāng)?shù)值模擬模型進(jìn)行耦合建立模擬優(yōu)化模型來解決該問題。求解水資源優(yōu)化配置的算法主要包括粒子群算法[4]、人工魚群算法[5]、模糊多目標(biāo)線性規(guī)劃[6]、遺傳算法[7]及魯棒規(guī)劃方法[8]等，由于它們能夠處理較為復(fù)雜的多目標(biāo)問題，因此在水資源優(yōu)化配置中得到廣泛應(yīng)用。Davijani等[9]采用遺傳算法和粒子群算法兩種啟發(fā)式算法，將水資源合理分配給工業(yè)、農(nóng)業(yè)和市政部門。譚倩等[8]采用魯棒規(guī)劃方法，對多個目標(biāo)進(jìn)行權(quán)重分析以便于決策者確定最優(yōu)方案。

部分研究通過構(gòu)建耦合地下水?dāng)?shù)值模擬模型和灌區(qū)水資源模擬優(yōu)化模型來進(jìn)行灌區(qū)水資源優(yōu)化配置。目前常用的地下水流數(shù)值模擬軟件主要有MODFLOW、FloPy[10]、GMS[11]、FEFLOW[12]等。Sedki等[13]將地下水非穩(wěn)定流模擬模型與基于遺傳算法構(gòu)建的優(yōu)化模型相結(jié)合，構(gòu)建了多目標(biāo)地下水資源可持續(xù)管理模型，在尋找最佳抽水方案的同時最大程度降低了對環(huán)境的不利影響。高玉芳等[14]為解決沿海地下水超采所引發(fā)的海水入侵問題，將MODFLOW的變密度地下水流及溶質(zhì)運(yùn)移模型與人工魚群算法相耦合，構(gòu)建了地下水模擬優(yōu)化管理模型。宋健[15]基于多目標(biāo)進(jìn)化算法構(gòu)建了復(fù)雜地下水系統(tǒng)的模擬優(yōu)化模型，用來解決濱海地下水管理和地下水污染監(jiān)測網(wǎng)優(yōu)化問題。為了解決實(shí)際應(yīng)用中多個目標(biāo)函數(shù)(大于3)的高維多目標(biāo)優(yōu)化問題，Deb等[16]提出了基于參考點(diǎn)的第3代非支配排序遺傳算法(NSGA-Ⅲ)，Blank等[17]基于Python語言開發(fā)了一個多目標(biāo)遺傳算法框架Pymoo。金江躍等[18]通過假想案例構(gòu)建了一個多目標(biāo)地下水模擬優(yōu)化模型。當(dāng)前灌區(qū)的水資源優(yōu)化配置傾向于將優(yōu)化配置結(jié)果代入地下水模擬模型進(jìn)而評估優(yōu)化配置方案對地下水的影響，屬于單向傳輸過程，而能夠動態(tài)雙向傳輸?shù)倪z傳算法與地下水耦合模型目前多應(yīng)用于地下水超采造成的海水入侵問題，針對灌區(qū)地表水地下水聯(lián)合優(yōu)化配置進(jìn)行緊密耦合的研究較少。因此本文構(gòu)建了基于NSGA-Ⅲ和FloPy的灌區(qū)水資源多目標(biāo)模擬優(yōu)化模型，實(shí)現(xiàn)了優(yōu)化算法與模擬模型之間的緊密耦合和動態(tài)的數(shù)據(jù)交換，并通過一個假想算例驗(yàn)證了該模擬優(yōu)化模型在渠井結(jié)合灌區(qū)的有效性和可靠性，為灌區(qū)水資源優(yōu)化配置和可持續(xù)高效安全利用提供了強(qiáng)有力的技術(shù)手段。

1 研究方法

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

MODFLOW是一個采用有限差分法進(jìn)行地下水流數(shù)值模擬的軟件，它由一個主程序和若干個相對獨(dú)立的子程序包構(gòu)成，每個子程序包又包含數(shù)個子模塊來分別完成模擬的一部分[19]。例如用井Well包模擬抽水井和注水井對地下水的影響。

FloPy是Bakker等[10]基于Python編程語言開發(fā)的可以運(yùn)行MODFLOW、MT3D、SEAWAT等模型的地下水建模工具，它可以通過Python語言構(gòu)建和運(yùn)行地下水?dāng)?shù)值模擬模型，且便于和優(yōu)化算法耦合構(gòu)建地下水模擬優(yōu)化模型。Python是目前功能強(qiáng)大的軟件開發(fā)語言之一，可以借助擴(kuò)展包完成繪圖、數(shù)據(jù)分析和優(yōu)化等過程，具有腳本語言中最強(qiáng)大的類庫和更完整的數(shù)據(jù)結(jié)構(gòu)集。基于Python語言開發(fā)的FloPy具有創(chuàng)建新模型、導(dǎo)入運(yùn)行已有模型、讀取和繪制模擬結(jié)果等強(qiáng)大的功能。在Flopy中建立MODFLOW模型通常包含以下八個步驟：(1) 導(dǎo)入FloPy包；(2) 創(chuàng)建MODFLOW模型對象；(3) 設(shè)置模型的基本參數(shù)離散化、邊界條件、水文地質(zhì)參數(shù)和圖層類型等；(4) 添加對應(yīng)的子程序包(例如井Well、補(bǔ)給Recharge、河流River)；(5) 定義MODFLOW水頭求解方案；(6) 定義MODFLOW需要保存文件的輸出類型；(7) 生成MODFLOW輸入文件并調(diào)用MODFLOW獲得求解方案；(8) 讀取輸出文件并可視化進(jìn)行分析。

import flopy

dis=flopy.modflow.ModflowDis(model, nlay, nrow, ncol, top=ztop, botm=botm)

pcg=flopy.modflow.ModflowPcg(mf)

hds=bf.HeadFile(‘Example’+“.hds”)

但對于邊界條件和水文地質(zhì)情況復(fù)雜的區(qū)域，難以直接采用FloPy構(gòu)建模型，F(xiàn)loPy提供了功能強(qiáng)大的API接口來導(dǎo)入其他軟件構(gòu)建的MODFLOW模型。如基于GMS構(gòu)建的模型在選擇保存為通用的MODFLOW格式后，即可通過FloPy的fmp.Modflow.load語句導(dǎo)入，model = fpm.Modflow.load(‘Example.mfn’, version=‘mf2005’,verbose=False,check=False)。FloPy可以正確地識別定義模型的離散化、活動單元格、頂?shù)装甯叱獭⑺μ匦院推渌磪R項(xiàng)，進(jìn)而修改和運(yùn)行模型。FloPy兼容MODFLOW-2000、2005、NWT、USG、LGR和最新版本的MODFLOW6。

FloPy可以通過編程的方式修改所導(dǎo)入模型的滲透系數(shù)、井模塊、補(bǔ)給模塊以及河流模塊等，也可以便捷地讀取地下水位模擬結(jié)果來計算遺傳算法中適應(yīng)度函數(shù)的值，進(jìn)而實(shí)現(xiàn)與多目標(biāo)遺傳算法的緊密耦合。

1.2 多目標(biāo)優(yōu)化框架Pymoo

Pymoo是為解決多目標(biāo)優(yōu)化的問題，基于Python語言開發(fā)的遺傳算法框架，可以通過模塊化和基于對象的編程來實(shí)現(xiàn)擴(kuò)展。主要包括優(yōu)化問題構(gòu)建、優(yōu)化算法實(shí)現(xiàn)以及數(shù)據(jù)分析在內(nèi)的三個主要的功能模塊。Pymoo提供了相關(guān)的性能指標(biāo)來衡量多目標(biāo)優(yōu)化算法獲得結(jié)果的質(zhì)量并將其可視化輸出。

本研究中選用Deb在2014年提出的多目標(biāo)遺傳算法NSGA-Ⅲ，它引入?yún)⒖键c(diǎn)機(jī)制，對于那些非支配并且接近參考點(diǎn)的種群個體進(jìn)行保留，來維持種群的多樣性。假設(shè)M為目標(biāo)數(shù)目，P為每個方向上分割的份數(shù)，其計算公式為H=(M+P-1P)。

目前研究中，NSGA-Ⅲ對求解3-15個目標(biāo)函數(shù)問題有良好的適用性，該算法的基本框架如圖1所示。

圖1 NSGA-Ⅲ的基本框架

1.3 耦合模型的構(gòu)建

在FloPy與Pymoo耦合時，需要在Pymoo中調(diào)用FloPy模擬結(jié)果來計算與地下水位和地下水開采費(fèi)用有關(guān)的目標(biāo)函數(shù)值，并讀取地下水位來判斷開采方案是否可行(地下水位是否在合理范圍內(nèi))，將FloPy作為能被Pymoo調(diào)用的子程序，實(shí)現(xiàn)二者之間數(shù)據(jù)的無縫傳輸。本研究采用Pymoo框架中的遺傳算法NSGA-Ⅲ進(jìn)行多目標(biāo)優(yōu)化模型的求解，其尋優(yōu)是一個循環(huán)過程。當(dāng)循環(huán)開始時，需要根據(jù)所設(shè)置的種群大小，對決策變量中各參數(shù)進(jìn)行初始化，形成第一代種群。將該種群作為可行解集，其中的每一個體都需要調(diào)用一次子程序FloPy進(jìn)行對應(yīng)方案下的地下水模擬計算，得到相應(yīng)的模擬結(jié)果。然后Pymoo將根據(jù)模擬結(jié)果判斷該方案是否滿足設(shè)置的約束條件，并計算目標(biāo)函數(shù)值，對各可行解進(jìn)行評價，得到新一代的種群。重復(fù)上述步驟，直至達(dá)到收斂條件，耦合模型的構(gòu)建流程如圖2所示。

圖2 耦合模型構(gòu)建流程圖

2 算例研究

2.1 算例背景

由于實(shí)際灌區(qū)情況非常復(fù)雜，缺乏構(gòu)建模擬優(yōu)化模型所需的基礎(chǔ)數(shù)據(jù)(如具體的井位)，因此本文在對西北干旱半干旱地區(qū)的典型渠井結(jié)合灌區(qū)——陜西省寶雞峽灌區(qū)進(jìn)行調(diào)研的基礎(chǔ)上，綜合考慮灌區(qū)實(shí)際條件設(shè)計了一個包含若干水文地質(zhì)分區(qū)的典型灌域來進(jìn)行模擬優(yōu)化模型的構(gòu)建驗(yàn)證和結(jié)果的分析論述。該區(qū)域長5 km，寬3.5 km，在MODFLOW模型中將該區(qū)域在空間上剖分為50行50列的長方形網(wǎng)格，根據(jù)水文地質(zhì)條件和種植作物的不同分為5個區(qū)域。其中分區(qū)1、5和2、3、4的地下水埋深不同，來模擬寶雞峽灌區(qū)塬上和塬下兩個灌溉系統(tǒng)。假定研究區(qū)為非均質(zhì)各向同性的潛水含水層，地下水流向?yàn)閺奈鞯綎|，兩邊水頭差設(shè)為5 m。東西邊界概化為定水頭邊界，南北邊界概化為零流量邊界。灌域內(nèi)共布置50眼抽水井，其井網(wǎng)密度為3眼/km2，構(gòu)建的水文地質(zhì)概念模型如圖3所示。

圖3 水文地質(zhì)概念模型

采用2009年灌區(qū)某氣象站的逐日降水資料進(jìn)行計算，計算時段為2009年1月1日至2009年12月31日，計算時間步長和應(yīng)力期均設(shè)為1 d，基于GMS構(gòu)建了地下水非穩(wěn)定流數(shù)值模擬模型。將其保存為原始MODFLOW格式后導(dǎo)入FloPy用來計算不同方案下的地下水位。本研究將降水入滲補(bǔ)給、田間灌溉入滲補(bǔ)給、井灌回歸補(bǔ)給按照面狀補(bǔ)給形式補(bǔ)給地下水，地下水開采采用MODFLOW中的井模塊進(jìn)行模擬。算例相關(guān)的水文地質(zhì)參數(shù)和種植結(jié)構(gòu)及種植面積分別如表1和表2所示。

表1 水文地質(zhì)參數(shù)

表2 種植結(jié)構(gòu)及種植面積 單位：km2

2.2 決策變量

針對本研究問題，共設(shè)置19個決策變量，分別為抽取地下水量Gij以及各作物生長過程的供水比例anj。其中Gij為各分區(qū)抽取地下水量，m3；anj為各作物第n次灌水量與N次總灌水量之比；i為分區(qū)編號；j為作物類型編號。

(1)

2.3 目標(biāo)函數(shù)的建立

灌區(qū)水資源優(yōu)化模型中目標(biāo)之間往往存在競爭關(guān)系，所以需要考慮不同層次利益主體的目標(biāo)[16]。因此，本算例設(shè)置的優(yōu)化目標(biāo)是滿足最大抽水降深約束的條件下，實(shí)現(xiàn)下層決策主體(農(nóng)民)的目標(biāo)：作物總產(chǎn)量最大以及經(jīng)濟(jì)效益最大，并同時實(shí)現(xiàn)上層決策主體(灌區(qū)管理者)的目標(biāo)：地下水平均累計降深最小。

(1) 產(chǎn)量目標(biāo)：各分區(qū)作物總產(chǎn)量最大

(2)

式中：Z1是每公頃作物總產(chǎn)量，kg；i是分區(qū)編號，i=1，2，3，4，5；j是作物類型編號，j=1，2，3，4，5，分別是小麥，玉米，棉花，果樹，蔬菜；Qij為不同分區(qū)作物的總需水量，m3；Aij為不同分區(qū)作物的種植面積，hm2；f是描述產(chǎn)量與生育期耗水量的二次函數(shù)關(guān)系，通過查閱相關(guān)文獻(xiàn)確定灌溉量與作物產(chǎn)量之間的具體關(guān)系如公式(2)—公式(6)所示[20-24]：

(3)

(4)

(5)

(6)

(7)

式中：f1、f2、f3、f4、f5、Q1、Q2、Q3、Q4、Q5分別為冬小麥、夏玉米、果樹、棉花以及蔬菜的產(chǎn)量和全生育期需水量。

(2) 經(jīng)濟(jì)目標(biāo)：農(nóng)業(yè)用水產(chǎn)生的凈效益最大

(8)

式中：Pij為各作物的單價，元；Csi為渠道引水的價格，0.27元/m3；Cg1和Cg2分別塬上和塬下抽水灌溉的價格，0.16元/m3和0.11元/m3；Gij為各分區(qū)抽取的地下水量，m3。

(3) 地下水平均累計降深目標(biāo)

該指標(biāo)表示每日某網(wǎng)格的地下水位相較于前一天該網(wǎng)格地下水位下降值(不考慮前一天地下水位值低于當(dāng)天水位值的情況)的累計年平均值。

(9)

式中：k是各網(wǎng)格編號，k=1，2，3…2500；t是決策時段長度，t=1，2，3…365，d。

2.4 約束條件

(1) 地下水變幅約束：

hmin

(10)

式中：hmax為該區(qū)域平均合理地下水埋深的上限閾值，m；hmin為該區(qū)域平均合理地下水埋深的下限閾值，m。參考地下水配置的相關(guān)研究，將塬下地下水埋深的上限閾值設(shè)置為2.43 m，下限閾值設(shè)置為14.79 m；塬上地下水埋深的上限閾值設(shè)置為3.05，下限閾值設(shè)置為58.05 m。

(2) 需水約束：

ETmin,ij

(11)

式中：ETmin,ij為各分區(qū)不同作物的最小需水量，m3/hm2，本文取灌溉定額的60%；ETmax,ij為各分區(qū)不同作物的最大需水量，m3/hm2，本文取灌溉定額。

3 優(yōu)化結(jié)果及分析

3.1 參數(shù)設(shè)置和求解過程

NSGA-Ⅲ的參數(shù)設(shè)置：種群大小N=200；優(yōu)化搜索代數(shù)為100；使用二進(jìn)制交叉和多項(xiàng)式變異，其中交叉率Pc=0.7，變異率Pm=0.05。

多目標(biāo)優(yōu)化算法求解結(jié)果可用超體積指數(shù)(Hypervolume)來評價，其表示由非劣解集中的個體與參考點(diǎn)在目標(biāo)空間中所圍成的超立方體的體積，以評價解集的收斂性與分布性，該指標(biāo)評價方法與Pareto保持一致，指標(biāo)越大，表明Pareto前沿的占優(yōu)性越強(qiáng)。算例的Hypervolume指標(biāo)見圖4，初始種群在第40代左右可達(dá)收斂水平，并逐步尋找最優(yōu)解集。

圖4 HV評價指標(biāo)

3.2 優(yōu)化方案分析

水資源多目標(biāo)優(yōu)化不能得到單一的全局最優(yōu)解，而是由多個非劣解構(gòu)成的非劣解集。其中包含傾向于不同目標(biāo)的多個方案，本文從非劣解集中選擇三個典型方案進(jìn)行對比分析(見表3)，具體結(jié)果如圖5所示。由表3可以看出，三個典型方案分別側(cè)重于產(chǎn)量最大、地下水平均降深最小和經(jīng)濟(jì)效益最大三個目標(biāo)。

表3 三種方案的目標(biāo)值

圖5 三個水資源配置方案對比分析

優(yōu)化結(jié)果中，方案1作物產(chǎn)量最大的同時地下水平均累計降深最大，經(jīng)濟(jì)效益最小，不利于灌區(qū)的經(jīng)濟(jì)發(fā)展；方案2灌區(qū)地下水平均累計降深最小，作物產(chǎn)量和經(jīng)濟(jì)效益介于方案1與3之間，對地下水的可持續(xù)高效安全利用最有利；方案3灌區(qū)收益最大，作物產(chǎn)量最小，地下水平均累計降深介于方案1與2之間。大多數(shù)研究表明作物產(chǎn)量與經(jīng)濟(jì)效益呈正相關(guān)關(guān)系，但是在本研究中，由于灌區(qū)分為兩類灌溉系統(tǒng)-塬上和塬下，二者是競爭關(guān)系。塬上地下水埋深大，抽水費(fèi)用高于塬下，導(dǎo)致保證作物產(chǎn)量的同時經(jīng)濟(jì)效益較低。因此當(dāng)考慮經(jīng)濟(jì)效益最大時，應(yīng)優(yōu)先滿足價格最高的果樹和抽水費(fèi)用較低的塬下小麥、玉米的灌水需求；當(dāng)考慮產(chǎn)量最大時，應(yīng)優(yōu)先滿足每公頃產(chǎn)量最大的小麥、玉米以及蔬菜的灌水需求；當(dāng)考慮地下水平均累計降深最小時，應(yīng)優(yōu)先滿足塬下種植作物的灌水需求，塬上作物僅滿足最低需水量。

觀測時段內(nèi)，方案1地下水位最大變幅為3.6 m，最低水位為373.9 m；方案2地下水位最大變幅為1.9 m，最低水位為375.1 m；方案3地下水位最大變幅為2.5 m，最低水位為374.2 m(見圖6、圖7)。三種方案最低水位均出現(xiàn)在7月28日，原因是該時刻灌溉已持續(xù)90 d且需要同時灌溉玉米、棉花及蔬菜，地下水開采量大于補(bǔ)給量，造成地下水位嚴(yán)重下降。

圖6 不同方案地下水位變化趨勢

圖7 不同方案灌區(qū)地下水位分布圖

綜上所述，三個目標(biāo)之間是相互競爭的關(guān)系，若要改善一方面是以犧牲其它兩個目標(biāo)為代價。決策者可根據(jù)自己偏重方向判斷，若決策者考慮增大作物產(chǎn)量，則選擇方案1；若決策者側(cè)重于地下水的保護(hù)，則選擇方案2；若決策者傾向于增加大經(jīng)濟(jì)效益，則建議選擇方案3。

4 結(jié)論與討論

為實(shí)現(xiàn)灌區(qū)水資源多目標(biāo)優(yōu)化配置，本文基于NSGA-Ⅲ和FloPy構(gòu)建了一個緊密耦合的灌區(qū)水資源多目標(biāo)模擬優(yōu)化模型，并基于假想算例進(jìn)行計算分析，驗(yàn)證了模擬優(yōu)化模型的有效性和可靠性。所構(gòu)建的灌區(qū)水資源模擬優(yōu)化模型以作物產(chǎn)量最大、經(jīng)濟(jì)效益最大和地下水平均累計降深最小為目標(biāo)，在考慮地下水位變幅約束的同時利用遺傳算法NSGA-Ⅲ獲得該問題的非劣解集。算例研究表明，基于FloPy和Pymoo可以在灌區(qū)構(gòu)建緊密耦合的水資源模擬優(yōu)化模型，能夠同時處理經(jīng)濟(jì)、產(chǎn)量和地下水位3個目標(biāo)函數(shù)之間相互競爭的復(fù)雜關(guān)系，并獲得傾向不同目標(biāo)的多個水資源配置方案，便于決策者根據(jù)當(dāng)?shù)貙?shí)際情況進(jìn)行決策。

本研究僅構(gòu)建了3個目標(biāo)函數(shù)，沒有充分發(fā)揮NSGA-Ⅲ算法在求解多個目標(biāo)優(yōu)化問題方面的優(yōu)勢，在后續(xù)研究中將進(jìn)一步考慮生態(tài)、社會等方面的目標(biāo)，針對真實(shí)灌區(qū)構(gòu)建考慮多個目標(biāo)的水資源模擬優(yōu)化模型。此外，本文僅針對一個概化的典型灌域進(jìn)行地下水模擬優(yōu)化模型的構(gòu)建，運(yùn)算相對簡單，因此采用mf2005.exe執(zhí)行文件進(jìn)行地下水?dāng)?shù)值模擬模型的計算，在進(jìn)行大型灌區(qū)模型計算時，可采用并行運(yùn)算mf2k5_h5_parallel.exe執(zhí)行文件以節(jié)省計算時間。