基于水量水質(zhì)耦合模擬優(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；3.內(nèi)蒙古自治區(qū)農(nóng)牧業(yè)科學(xué)院 資源環(huán)境與可持續(xù)發(fā)展研究所，呼和浩特 010031）

0 引 言

我國(guó)作為農(nóng)業(yè)大國(guó)，由于特殊的地理和氣候條件，同時(shí)面臨著水資源時(shí)空分布不均、供需矛盾突出、水質(zhì)惡化和水資源利用率低的問(wèn)題，這使得農(nóng)業(yè)發(fā)展受到水資源的強(qiáng)烈制約[1,2]。地下水作為水資源的重要組成部分，具有重要的生態(tài)功能[3]。農(nóng)業(yè)地下水使用量占全國(guó)地下水用水總量的67%，個(gè)別地區(qū)高達(dá)80%～90%[4]。受氣候變化和人類開(kāi)發(fā)利用的影響，地下水的補(bǔ)排條件隨之改變，一些灌區(qū)的地下水動(dòng)態(tài)發(fā)生明顯變化，水質(zhì)不可避免受到一定程度的影響，部分地區(qū)出現(xiàn)了嚴(yán)重超采導(dǎo)致的地下沉降漏斗和農(nóng)業(yè)用肥用藥導(dǎo)致的污染擴(kuò)散[5,6]。結(jié)合我國(guó)地下污染情況[7,8]以及國(guó)家相關(guān)政策，綜合考慮灌溉對(duì)地下水水量水質(zhì)的影響，進(jìn)行地下水的綜合管理尤其是地下水污染防治迫在眉睫。而灌區(qū)的灌溉方式和用水量、作物產(chǎn)量、水質(zhì)演變和用水效率等相互影響和制約，形成了復(fù)雜的互饋關(guān)系：作物的高產(chǎn)需要充足的水肥藥供應(yīng)，大量抽取地下水和化肥農(nóng)藥的使用會(huì)導(dǎo)致地下水位的下降、水質(zhì)的加速惡化和作物水分生產(chǎn)率的下降。因此，灌區(qū)水資源優(yōu)化配置和可持續(xù)高效安全利用需要同時(shí)考慮作物產(chǎn)量、水量、水質(zhì)、用水效率等方面的目標(biāo)和約束。

灌區(qū)水資源優(yōu)化配置一直是水資源領(lǐng)域研究的熱點(diǎn)之一[9,10]，早期已經(jīng)有基于物理模型和優(yōu)化模型的研究和探索[11,12]，但并沒(méi)有將地下水?dāng)?shù)值模擬模型與優(yōu)化算法進(jìn)行緊密耦合。近年來(lái)，已有部分研究耦合二者構(gòu)建了模擬優(yōu)化模型來(lái)解決灌區(qū)的水資源優(yōu)化配置問(wèn)題[13-15]，采用了多種多目標(biāo)優(yōu)化算法[16-18]，如非支配排序遺傳算法NSGA-Ⅱ和NSGA-Ⅲ等。NSGA-Ⅱ已從理想算例的應(yīng)用[19,20]發(fā)展至實(shí)際灌區(qū)案例[21-23]，NSGA-Ⅲ也已經(jīng)有應(yīng)用于灌區(qū)的案例[24,25]，但鮮有對(duì)比二者在灌區(qū)多目標(biāo)水資源優(yōu)化配置中應(yīng)用效果的研究，此外，已有研究表明灌溉與地下水污染物擴(kuò)散之間的關(guān)系密切[26-28]，但目前的模擬優(yōu)化模型中均沒(méi)有考慮灌溉對(duì)地下水溶質(zhì)運(yùn)移的影響，有加劇污染擴(kuò)散的風(fēng)險(xiǎn)，亟需構(gòu)建水量水質(zhì)耦合模擬優(yōu)化模型來(lái)進(jìn)行灌區(qū)水資源優(yōu)化配置。

因此本文考慮灌溉收益、作物水分生產(chǎn)率、地下水位變幅和污染擴(kuò)散四個(gè)目標(biāo)，構(gòu)建了緊密耦合的綜合考慮水量水質(zhì)的多目標(biāo)模擬優(yōu)化模型，并將NSGA-Ⅱ與NSGA-Ⅲ的實(shí)際應(yīng)用效果進(jìn)行定量對(duì)比分析，最后采用多標(biāo)準(zhǔn)分析方法TOPSIS確定推薦的灌溉方案。

1 研究理論框架

1.1 基于FloPy的地下水?dāng)?shù)值模擬模型

地下水?dāng)?shù)值模擬模型為精確定量刻畫(huà)地下水水流和溶質(zhì)運(yùn)移提供了強(qiáng)有力的技術(shù)手段[29-31]，本文采用二維潛水含水層非穩(wěn)定地下水流方程來(lái)描述研究區(qū)的地下水運(yùn)動(dòng)過(guò)程：

式中：K為含水層的滲流系數(shù)；h為區(qū)內(nèi)各點(diǎn)潛水水位；b為含水層底板高程；μ為給水度；t為計(jì)算時(shí)間；ω為含水層的源匯項(xiàng)；Ω 為研究區(qū)域；Γ1、Γ2為第一、第二類邊界；q(x,y,t)為流量邊界；x,y為橫縱坐標(biāo)。

考慮地下水的對(duì)流、彌散、源匯項(xiàng)等，采用單一化學(xué)組分的三維遷移偏轉(zhuǎn)方程來(lái)描述污染源的擴(kuò)散情況：

式中：R為延遲因子；C為溶解濃度；為吸附濃度；xi為坐標(biāo)，xi=x,y,z；qi為x,y,z三個(gè)方向上的達(dá)西速度分量，qi=qx,qy,qz；Dij為二階彌散系數(shù)張量；qs為源匯處單位體積含水層的流量；Cs為源匯的濃度；λ1、λ2為溶解相、吸附相的反應(yīng)速率常數(shù)；θ為孔隙度；ρb為孔隙介質(zhì)的體積密度。

為便于與多目標(biāo)遺傳算法進(jìn)行耦合，本文采用Flopy 來(lái)構(gòu)建MODFLOW 和MT3D 模型。FloPy 是由基于python 語(yǔ)言開(kāi)發(fā)的第三方庫(kù)[32]，可以實(shí)現(xiàn)MODFLOW、MT3D、SEAWAT 等模型的讀寫(xiě)、構(gòu)建、運(yùn)行以及結(jié)果顯示，并支持最新的MODFLOW6 版本，在FloPy 構(gòu)建、運(yùn)行地下水?dāng)?shù)值模擬模型時(shí)可以調(diào)用支持并行計(jì)算的mf2k_h5_dbl_parallel. exe、mf2k_h5_parallel.exe等程序來(lái)實(shí)現(xiàn)地下水?dāng)?shù)值模擬模型的并行運(yùn)算。

1.2 基于Pymoo的多目標(biāo)模擬優(yōu)化框架

Pymoo 是由Blank 等[33]人基于Python 語(yǔ)言開(kāi)發(fā)的一個(gè)能實(shí)現(xiàn)多目標(biāo)優(yōu)化功能的第三方庫(kù)，它可以采用不同的算法解決多目標(biāo)優(yōu)化問(wèn)題，并提供豐富的可視化結(jié)果。其中遺傳算法是用來(lái)分析和解決多目標(biāo)優(yōu)化問(wèn)題的一種進(jìn)化算法，其核心就是協(xié)調(diào)各個(gè)目標(biāo)函數(shù)之間的關(guān)系，找出使得各個(gè)目標(biāo)函數(shù)都盡可能達(dá)到比較大的（或者比較小的）函數(shù)值的最優(yōu)解集。在多目標(biāo)優(yōu)化的遺傳算法中，由Deb 等[34]提出的NSGA-Ⅱ是影響最大和應(yīng)用范圍最廣的多目標(biāo)遺傳算法之一。2014年Deb等[35]又改進(jìn)后提出NSGA-Ⅲ，也就是根據(jù)NSGA-Ⅱ框架提出了一種基于參考點(diǎn)的多目標(biāo)進(jìn)化算法，增加了非優(yōu)勢(shì)群體成員的權(quán)重。NSGA 適用于許多具有3到15個(gè)目標(biāo)的復(fù)雜優(yōu)化問(wèn)題。Pymoo 中采用超立方體積指標(biāo)值（Hypervolume，簡(jiǎn)稱Hv）作為評(píng)價(jià)指標(biāo)，Hv值收斂時(shí)可認(rèn)為求解所得非劣解集已經(jīng)趨于合理。

1.3 基于理想解的多層次分析方法

NSGA 計(jì)算得到的非劣解集包括多個(gè)優(yōu)化方案，生產(chǎn)實(shí)踐中需要從中選出最終的推薦方案。本文采用一種多標(biāo)準(zhǔn)決策方法TOPSIS(Technique for Order Preference by Similarity to Ideal Solution)，該方法已在水資源管理領(lǐng)域被廣泛采用[36,37]。利用TOPSIS，每一方案都被賦予一個(gè)相對(duì)分?jǐn)?shù)。然后，將各備選方案按其相對(duì)分?jǐn)?shù)進(jìn)行排序，并選擇排名最高的方案作為推薦方案。

TOPSIS 方法的主要步驟包括[38]：①構(gòu)建判斷矩陣；②對(duì)矩陣進(jìn)行歸一化處理，得到歸一化矩陣；③根據(jù)熵的定義確定評(píng)價(jià)指標(biāo)的熵；④評(píng)價(jià)指標(biāo)的熵權(quán)，計(jì)算各指標(biāo)權(quán)重集；⑤確定理想解與負(fù)理想解；⑥計(jì)算各方案與理想解的相對(duì)接近程度：

式中：S+i為各方案與理想解的接近程度；S-i為各方案與負(fù)理想解的接近程度；Ci∈[0,1]，且該方案Ci值越大越好。

1.4 總體框架

本文通過(guò)構(gòu)建多目標(biāo)模擬優(yōu)化模型的方法來(lái)進(jìn)行灌區(qū)水資源優(yōu)化配置，其中的決策變量有4 個(gè)，分別為3 個(gè)灌季的灌溉用水量和灌溉水源的渠井用水比例。采用經(jīng)驗(yàn)公式計(jì)算不同灌水條件下的作物產(chǎn)量，調(diào)用FloPy 運(yùn)行MODFLOW 和MT3D 模型來(lái)計(jì)算抽水灌溉對(duì)地下水位和污染物運(yùn)移的影響，采用NSGA 算法來(lái)求解所構(gòu)建的多目標(biāo)模擬優(yōu)化模型，在得到非劣解集后采用TOPSIS 方法來(lái)確定最終的推薦方案。本文的總體框架如圖1所示。

圖1 多目標(biāo)模擬優(yōu)化模型框架Fig.1 The framework of multi-objective simulation optimization model

2 典型灌區(qū)案例概化及模擬優(yōu)化模型的構(gòu)建

2.1 典型灌區(qū)案例概化

本文根據(jù)干旱半干旱地區(qū)牧草灌溉的野外實(shí)地調(diào)查概化出典型案例。灌區(qū)地勢(shì)平坦，渠井等灌溉設(shè)施完善，同一地塊可根據(jù)渠道來(lái)水情況選擇采用渠水或井水進(jìn)行灌溉。如圖2所示，該區(qū)域?yàn)橛梢豢跈C(jī)井控制的灌溉范圍，其面積約為36 hm2，將其概化為一個(gè)正方形區(qū)域，所開(kāi)采的地下水為一厚度為35 m 的均質(zhì)各向同性潛水含水層。含水層頂板高程為635 m，底板高程為600 m，潛水水位631 m。該含水層上層為弱透水層，可接受降水和灌溉水的入滲補(bǔ)給，區(qū)域內(nèi)的作物為苜蓿，每年按照灌溉制度和灌溉習(xí)慣分為3 個(gè)灌季進(jìn)行灌溉。灌溉水源來(lái)自襯砌完好的渠水（地表水）和機(jī)井抽取的地下水。機(jī)井位于區(qū)域的中心位置。在典型灌區(qū)內(nèi)設(shè)置一面狀污染源，另外結(jié)合垃圾填埋場(chǎng)實(shí)際污染物類型特征，選取硫化物作為本次模擬分析的污染物，按照地下水質(zhì)量標(biāo)準(zhǔn)（GB/T 14848--2017）中Ⅲ類水（主要用于生活和農(nóng)業(yè)用水）水質(zhì)要求為標(biāo)準(zhǔn)（其中硫化物濃度不大于0.02 mg/L）判斷灌區(qū)是否受到污染。并結(jié)合灌溉抽水的不同情況來(lái)探究污染擴(kuò)散隨之發(fā)生的變化。本案例將灌溉用水中的渠井用水比例和3個(gè)灌季的灌溉用水量作為決策變量，分別根據(jù)平水年和枯水年的實(shí)測(cè)降水資料和灌溉情況考慮作物產(chǎn)量、污染物擴(kuò)散規(guī)模等目標(biāo)構(gòu)建水量水質(zhì)耦合的多目標(biāo)模擬優(yōu)化模型來(lái)進(jìn)行多目標(biāo)水資源優(yōu)化配置。

圖2 典型案例概化圖Fig.2 Schematic diagram of typical cases

2.2 地下水?dāng)?shù)值模擬模型構(gòu)建

根據(jù)典型案例概化的水文地質(zhì)概念模型如圖3所示。潛水含水層水頭西高東低。則模擬區(qū)南北可概化為由流線組成的零流量二類邊界，東西為由等水位線組成的一類水位邊界。在FloPy 中編制代碼來(lái)設(shè)定模型的初始條件、邊界條件和源匯項(xiàng)[24]（降水入滲補(bǔ)給、地表水灌溉入滲補(bǔ)給、井灌回歸補(bǔ)給、機(jī)井抽水），并輸入滲透系數(shù)、彈性給水度、污染源位置和污染物濃度等信息，分別構(gòu)建非穩(wěn)定流MODFLOW 水流模擬模型和MT3D溶質(zhì)運(yùn)移模擬模型。為保證模擬精度，反映降水和抽水灌溉對(duì)地下水的影響，模型的應(yīng)力期設(shè)置為1 d，模擬時(shí)段為一個(gè)日歷年，模擬區(qū)域橫向和縱向均剖分為60 等份，共計(jì)3 600個(gè)單元格（圖3）。構(gòu)建好的模型可以定量模擬降水入滲、渠水灌溉、抽水灌溉等對(duì)地下水位和和污染物擴(kuò)散運(yùn)移的影響。

圖3 典型模擬區(qū)水文地質(zhì)概念模型概化圖Fig.3 Schematic diagram of hydrogeological conceptual model in typical simulation area

2.3 水量水質(zhì)耦合的多目標(biāo)模擬優(yōu)化模型的構(gòu)建

模擬優(yōu)化模型中考慮的4 個(gè)決策變量設(shè)定為（Q1,Q2,Q3,x）,其中Qi（i=1,2,3）為各個(gè)灌季的灌溉水量（m3/hm2），x為灌溉用水中地下水所占的比例。灌溉用水量根據(jù)灌溉制度在灌溉時(shí)間內(nèi)進(jìn)行逐日平均分配并考慮入滲補(bǔ)給后作為源匯項(xiàng)輸入地下水?dāng)?shù)值模擬模型中。

模擬優(yōu)化共考慮4個(gè)目標(biāo)，分別為灌溉效益、作物水分生產(chǎn)率、灌溉對(duì)地下水動(dòng)態(tài)的影響和污染物擴(kuò)散范圍。約束條件為灌溉抽水量的上下限。其中4個(gè)目標(biāo)函數(shù)以及約束條件表述如下。

（1）效益最大化目標(biāo)：

其中產(chǎn)量計(jì)算公式采用Li等[39]于2017年實(shí)驗(yàn)所得產(chǎn)量公式：

式中：I1、I2、I3為第一次、第二次、第三次的實(shí)際灌溉量，mm；A為灌區(qū)面積，hm2；Ys為苜?？偖a(chǎn)量，kg；P、Ps、Pg分別為苜蓿單價(jià)，元/kg、地表水單價(jià)，元/m3、地下水用電單價(jià)，元/kWh；Y1、Y2、Y3分別為第一次、第二次、第三次刈割的干苜蓿產(chǎn)量，kg/hm2。

（2）作物水分生產(chǎn)率目標(biāo)[40]：

式中：θ為作物水分生產(chǎn)率，kg/m3；λ為復(fù)種指數(shù)；G為糧食平均單產(chǎn)，kg/hm2；IR為平均灌溉用水量，m3/hm2；Pe為糧食作物生長(zhǎng)期內(nèi)實(shí)際補(bǔ)充到作物根分布層可被其利用的有效降水量，mm。

（3）地下水位變化目標(biāo)，統(tǒng)計(jì)日尺度的地下水水位降深之和。

調(diào)用MODFLOW 地下水?dāng)?shù)值模擬組件計(jì)算得到模擬區(qū)單元格水位矩陣。再由水位矩陣中的數(shù)據(jù)計(jì)算求得每個(gè)單元格的水位降深。

式中：k為各網(wǎng)格編號(hào)，k=1，2，…，3 600；t為決策時(shí)段長(zhǎng)度，t=1，2，…，365，d。hkt為地下水模型計(jì)算所得逐日的各單元格的水位數(shù)據(jù)，m。

（4）污染的擴(kuò)散面積目標(biāo)：

式中：S為總擴(kuò)散面積，hm2；A為灌區(qū)總面積，hm2；s為MT3D 模型根據(jù)地下水運(yùn)動(dòng)過(guò)程計(jì)算所得的污染物擴(kuò)散的單元格數(shù)量，會(huì)隨著不同的地下水抽水過(guò)程（即Q1,Q2,Q3的大小以及地下水占比x）而變化；n為總單元格數(shù)。

（5）約束條件。設(shè)置不同水平年每個(gè)灌季灌水量上下限以及總灌水量約束：

平水年：

枯水年：

式中：Q1、Q2、Q3分別為3個(gè)灌季的灌水量，m3/hm2

3 模型求解

根據(jù)已有的4個(gè)目標(biāo)，結(jié)合地下水?dāng)?shù)值模型，利用多目標(biāo)優(yōu)化工具Pymoo建立以遺傳算法框架為基礎(chǔ)的多目標(biāo)管理優(yōu)化模型對(duì)模型進(jìn)行求解。經(jīng)過(guò)多次反復(fù)試驗(yàn)，在既保證優(yōu)化結(jié)果的收斂性，又兼顧模型運(yùn)行一次所需時(shí)間的前提下，設(shè)定種群數(shù)為50，遺傳代數(shù)為50 代，交叉率0.8，變異率0.05。本次運(yùn)行求解所用電腦CPU 為Intel(R) Core(TM) i5-7200 U，頻率為2.50 GHz，運(yùn)行內(nèi)存為8 GB，進(jìn)行一次完整的求解過(guò)程（50×50）大約需要8 h。

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

4.1 NSGA-Ⅱ和NSGA-Ⅲ應(yīng)用效果對(duì)比分析

分別將目標(biāo)函數(shù)灌區(qū)效益、作物水分生產(chǎn)率、地下水位變化3 個(gè)模型計(jì)算所得的目標(biāo)函數(shù)值設(shè)置為x、y、z軸，污染擴(kuò)散面積用顏色標(biāo)記，繪制枯水年和平水年的非劣解集四維空間分布圖（圖4和圖5）。從模擬優(yōu)化的結(jié)果圖中可以初步看出，NSGA-Ⅱ在模擬優(yōu)化過(guò)程中所能保存的非劣解集數(shù)量明顯多于NSGA-Ⅲ。反復(fù)驗(yàn)證分析后發(fā)現(xiàn)，NSGA-Ⅱ是按照所設(shè)定的種群數(shù)量來(lái)保存非劣解集的，即本次所設(shè)定的每代50個(gè)種群，模擬優(yōu)化所留下的非劣解集的數(shù)量仍然是50；但由于NSGA-Ⅲ增加了非優(yōu)勢(shì)群體成員的權(quán)重，種群之間的競(jìng)爭(zhēng)淘汰會(huì)更加激烈，這使得NSGA-Ⅲ算法下所能保存的非劣解集的數(shù)量只有20個(gè)左右，明顯少于NSGA-Ⅱ。

圖4 枯水年優(yōu)化結(jié)果空間分布圖Fig.4 Spatial distribution of optimum results of dry year

圖5 平水年優(yōu)化結(jié)果空間分布圖Fig.5 Spatial distribution of optimum results of normal year

按照效益最大，作物水分生產(chǎn)率最高，地下水位變化最小分析，最優(yōu)的非劣解集應(yīng)該出現(xiàn)在空間的左下角，其中顏色越偏向紅色說(shuō)明污染擴(kuò)散越嚴(yán)重，越偏向綠色說(shuō)明污染控制的越好?？傮w來(lái)說(shuō)非劣解集的空間分布符合4個(gè)目標(biāo)的權(quán)衡關(guān)系。再分別對(duì)比枯水年與平水年條件下兩種遺傳算法保留的非劣解集可以看出，其空間分布基本一致，各目標(biāo)函數(shù)變換范圍也趨于相同，但是如之前所分析，NSGA-Ⅲ會(huì)保留更少的精英方案，非劣解集中地下水位變化和污染擴(kuò)散面積這兩個(gè)目標(biāo)函數(shù)值明顯偏大的方案顯著減少。

兩種年份下的NSGA-Ⅱ和NSGA-Ⅲ的Hv值分布圖如圖6和圖7 所示，其中在進(jìn)行50 代的求解過(guò)程中NSGA-Ⅱ?qū)τ诮Y(jié)果的收斂效果明顯不如NSGA-Ⅲ，其中枯水年優(yōu)化過(guò)程中NSGA-Ⅱ前30代繁衍hv值波動(dòng)較大，而NSGA-Ⅲ在繁衍不到10 代的情況下即已趨于穩(wěn)定；平水年優(yōu)化過(guò)程中NSGA-Ⅱ甚至出現(xiàn)未收斂的情況，而NSGA-Ⅲ依舊可以在繁衍30～40 代后基本趨于穩(wěn)定收斂。由此可以看出在本次模擬優(yōu)化模型的求解效率上NSGA-Ⅱ低于NSGA-Ⅲ，不能實(shí)現(xiàn)快速收斂。

圖6 枯水年優(yōu)化結(jié)果Hv值Fig.6 Hypervolume value of optimization results in dry year (NSGA -Ⅱ(left) NSGA -Ⅲ(right))

圖7 平水年優(yōu)化結(jié)果Hv值Fig.7 Hypervolume value of optimization results in normal year

綜合以上來(lái)看，NSGA-Ⅲ保存的非劣解少于NSGA-Ⅱ，但是NSGA-Ⅲ可以更快的使模擬結(jié)果收斂而尋找最優(yōu)解，在本次模擬種群數(shù)目為50 進(jìn)行50 代進(jìn)化計(jì)算所得的優(yōu)化方案中，其效率是明顯優(yōu)于NSGA-Ⅱ的。所以后續(xù)的方案篩選將主要在枯水年和平水年的NSGA-Ⅲ的優(yōu)化結(jié)果中進(jìn)行。

4.2 多標(biāo)準(zhǔn)分析法選擇推薦方案

由NSGA-Ⅲ計(jì)算得到的非劣解集中不同水平年的灌溉方案分別對(duì)應(yīng)4個(gè)目標(biāo)函數(shù)值（表1和表2），其中平水年保留的非劣解集包含20 個(gè)方案，枯水年包含18 個(gè)方案。各方案目標(biāo)函數(shù)值對(duì)比如圖8和圖9所示。

表2 NSGA-Ⅲ優(yōu)化枯水年方案結(jié)果Tab.2 Dry year results of NSGA-Ⅲ optimized solution

圖8 平水年非劣解集目標(biāo)函數(shù)值Fig.8 Objective function value of Pareto front in normal year

圖9 枯水年非劣解集目標(biāo)函數(shù)值Fig.9 Objective function value of Pareto front in dry year

為更好的看出非劣解集中每個(gè)方案權(quán)衡4個(gè)不同目標(biāo)的表現(xiàn)，繪制平水年各方案的百分比堆積圖如圖8和圖9所示。

從各方案對(duì)比的圖表中可以看出，每個(gè)方案對(duì)不同目標(biāo)的傾向不同，難以直接評(píng)價(jià)某一方案的優(yōu)劣程度。因此本文采用TOPSIS 分析評(píng)價(jià)方法，利用各方案中目標(biāo)函數(shù)值的最值自動(dòng)生成默認(rèn)的理想解與負(fù)理想解，然后由公式（4）依據(jù)各方案與理想解和負(fù)理想解的距離進(jìn)行綜合打分，依據(jù)打分結(jié)果分別選出了枯水年和平水年的推薦方案如表3所示。

表3 不同水平年的推薦方案Tab.3 Optimal solution in different standard years

結(jié)合前文所提到的問(wèn)題，重點(diǎn)分析兩個(gè)推薦方案灌溉模式下的地下水水位變化以及污染擴(kuò)散情況，利用地下水?dāng)?shù)值模擬模型模擬其一整年的地下水位變化與污染擴(kuò)散的最終結(jié)果并作圖比較分析。

灌溉方案A、B 的抽水井和污染源處地下水位變化如圖10和圖11 所示，抽水井處水位隨灌溉進(jìn)行變化明顯，污染源處水位也會(huì)隨灌溉的進(jìn)行發(fā)生一定的變化，由此可以看出灌溉導(dǎo)致的地下水位變化抽水井處向四周逐漸減弱，但其波動(dòng)也已經(jīng)影響到污染源處的地下水。其中不同水平年條件下的灌溉方案A、B 灌溉時(shí)典型灌區(qū)抽水井和污染源處的地下水位逐日變化趨勢(shì)基本一致，但是枯水年的灌溉方案B地下水位變化幅度相對(duì)更大。為更進(jìn)一步的體現(xiàn)灌溉對(duì)灌區(qū)地下水水位變化的影響，再分別選取其中地下水位波動(dòng)最大的一日，作出整個(gè)灌區(qū)的地下水水位圖（圖12）。

圖10 平水年地下水位變化Fig.10 The groundwater level fluctuation in normal year

圖11 枯水年地下水位變化Fig.11 The groundwater level fluctuation level in dry year

圖12 方案A和方案B最大單日灌溉量時(shí)地下水位Fig.12 The groundwater level at maximum daily irrigation volume in solution A and solution B

而在實(shí)際模擬一整年的灌溉抽水后污染物的總體擴(kuò)散情況如圖13 所示，由于本案例的典型灌區(qū)采用的節(jié)水灌溉模式中的噴灌模式，總體灌溉水量不大，又引進(jìn)了地表水水源極大的減輕了地下水灌溉壓力，所以最后的污染物擴(kuò)散情況控制比較好。

圖13 方案A和方案B污染擴(kuò)散情況Fig.13 Pollution diffusion of solution A and solution B

而若不引進(jìn)地表水水源且一直按照各水平年最高灌溉定額灌溉，其對(duì)應(yīng)的最大單日灌溉抽水時(shí)的水位變化和最終的污染物擴(kuò)散如圖14 和圖15 所示。其污染物的擴(kuò)散面積要明顯大于前面的兩個(gè)推薦方案。根據(jù)模型進(jìn)一步計(jì)算結(jié)果，平水年的推薦灌溉方案A 要比全地下水最大灌溉量灌溉方案的地下水位變化幅度減少了0.55 m（56%），污染擴(kuò)散面積減小了1.24 hm2（54%），枯水年的推薦灌溉方案B 要分別比全地下水灌溉方案的地下水位變化幅度減少了0.51 m（43%），污染擴(kuò)散面積減小了1.13 hm2（40%），模型對(duì)重點(diǎn)關(guān)注的灌區(qū)地下水水量水質(zhì)的優(yōu)化效果顯著。

圖14 全地下水灌溉最大單日灌溉量時(shí)地下水位（平水年、枯水年）Fig.14 The groundwater level at maximum daily irrigation volume under the condition of full groundwater irrigation (normal year and dry year)

圖15 全地下水灌溉污染物擴(kuò)散情況（平水年、枯水年）Fig.15 Pollution diffusion under the condition of full groundwater irrigation (normal year and dry year)

5 結(jié)論與討論

基于FloPy 和Pymoo 構(gòu)建了典型渠井結(jié)合灌區(qū)水量水質(zhì)耦合的多目標(biāo)模擬優(yōu)化模型。模型采用地下水?dāng)?shù)值模擬模型MODFLOW 和溶質(zhì)輸移模型MT3D 模擬地下水位變化和污染源擴(kuò)散過(guò)程。分別采用多目標(biāo)遺傳算法NSGA-Ⅱ和NSGA-Ⅲ考慮4個(gè)目標(biāo)（灌區(qū)收益最大化、作物水分生產(chǎn)率最大、地下水位變幅最小、污染擴(kuò)散面積最小）獲得了該復(fù)雜問(wèn)題的非劣解集。最終采用TOPSIS 方法獲得不同水平年的推薦方案。主要結(jié)論如下：①運(yùn)用基于Python 版本的FloPy 和Pymoo 構(gòu)建了耦合水量水質(zhì)的灌區(qū)水資源多目標(biāo)模擬優(yōu)化模型，優(yōu)化結(jié)果表明該模型可以較好的權(quán)衡多個(gè)目標(biāo)之間的競(jìng)爭(zhēng)關(guān)系；②構(gòu)建的模型可以使用Python 語(yǔ)言命令直接調(diào)用MODFLOW 和MT3D 的可執(zhí)行程序來(lái)運(yùn)行模型，且模型設(shè)置目標(biāo)函數(shù)以及修改參數(shù)簡(jiǎn)潔方便，可以直接使用代碼循環(huán)多次運(yùn)行模型，并能夠基于不同版本的遺傳算法進(jìn)行求解，顯著的提高了模型運(yùn)行調(diào)試和問(wèn)題求解的效率；③對(duì)比了不同版本遺傳算法（NSGA-Ⅱ、NSGA-Ⅲ）在灌區(qū)水資源多目標(biāo)優(yōu)化方面的實(shí)際應(yīng)用效果。結(jié)果表明，相比NSGA-Ⅱ，NSGA-Ⅲ會(huì)因?yàn)閷?duì)非優(yōu)勢(shì)種群的兼顧，使得競(jìng)爭(zhēng)更加激烈，以至于會(huì)丟失一部分的種群。但是NSGA-Ⅲ收斂效果更好、效率更高、提供的非劣解集更加穩(wěn)定；④以NSGA-Ⅲ優(yōu)化所得非劣解集為基礎(chǔ)，利用多標(biāo)準(zhǔn)分析方法（TOPSIS）選擇出了不同水平年的推薦灌溉方案。在保證收益和作物水分生產(chǎn)率的前提下，推薦灌溉方案明顯的減少了對(duì)地下水的影響（水位變化以及污染的擴(kuò)散），平枯水年分別可以最高減少地下水位變幅56%、43%，最高減少平枯水年的污染擴(kuò)散面積54%、40%。驗(yàn)證了多目標(biāo)模擬優(yōu)化模型在灌區(qū)水資源優(yōu)化配置方面的有效性與可靠性。

本研究從工具以及技術(shù)角度探究了多目標(biāo)遺傳算法在灌區(qū)進(jìn)行考慮多重因素的灌區(qū)水資源優(yōu)化配置，后續(xù)研究可以與實(shí)際情況相結(jié)合，擴(kuò)大研究范圍，耦合作物模型并考慮農(nóng)業(yè)農(nóng)村污染和能源消耗情況，進(jìn)一步構(gòu)建符合我國(guó)灌區(qū)實(shí)際情況的耦合作物、水量水質(zhì)和能源等多目標(biāo)的模擬優(yōu)化模型。而針對(duì)模擬優(yōu)化過(guò)程中數(shù)值模擬模型運(yùn)行時(shí)間較長(zhǎng)的情況，可以嘗試使用基于機(jī)器學(xué)習(xí)的替代模型，以提高模擬優(yōu)化模型運(yùn)行效率。