基于水量水质耦合模拟优化的渠井结合灌区多目标水资源优化配置模型与方法

苏振辉，降亚楠，2，吕婧妤，徐 超，陈 威，李 彬

（1.西北农林科技大学水利与建筑工程学院，陕西 杨凌 712100；2.西北农林科技大学 旱区农业水土工程教育部重点实验室，陕西 杨凌 712100；3.内蒙古自治区农牧业科学院 资源环境与可持续发展研究所，呼和浩特 010031）

0 引 言

我国作为农业大国，由于特殊的地理和气候条件，同时面临着水资源时空分布不均、供需矛盾突出、水质恶化和水资源利用率低的问题，这使得农业发展受到水资源的强烈制约[1,2]。地下水作为水资源的重要组成部分，具有重要的生态功能[3]。农业地下水使用量占全国地下水用水总量的67%，个别地区高达80%～90%[4]。受气候变化和人类开发利用的影响，地下水的补排条件随之改变，一些灌区的地下水动态发生明显变化，水质不可避免受到一定程度的影响，部分地区出现了严重超采导致的地下沉降漏斗和农业用肥用药导致的污染扩散[5,6]。结合我国地下污染情况[7,8]以及国家相关政策，综合考虑灌溉对地下水水量水质的影响，进行地下水的综合管理尤其是地下水污染防治迫在眉睫。而灌区的灌溉方式和用水量、作物产量、水质演变和用水效率等相互影响和制约，形成了复杂的互馈关系：作物的高产需要充足的水肥药供应，大量抽取地下水和化肥农药的使用会导致地下水位的下降、水质的加速恶化和作物水分生产率的下降。因此，灌区水资源优化配置和可持续高效安全利用需要同时考虑作物产量、水量、水质、用水效率等方面的目标和约束。

灌区水资源优化配置一直是水资源领域研究的热点之一[9,10]，早期已经有基于物理模型和优化模型的研究和探索[11,12]，但并没有将地下水数值模拟模型与优化算法进行紧密耦合。近年来，已有部分研究耦合二者构建了模拟优化模型来解决灌区的水资源优化配置问题[13-15]，采用了多种多目标优化算法[16-18]，如非支配排序遗传算法NSGA-Ⅱ和NSGA-Ⅲ等。NSGA-Ⅱ已从理想算例的应用[19,20]发展至实际灌区案例[21-23]，NSGA-Ⅲ也已经有应用于灌区的案例[24,25]，但鲜有对比二者在灌区多目标水资源优化配置中应用效果的研究，此外，已有研究表明灌溉与地下水污染物扩散之间的关系密切[26-28]，但目前的模拟优化模型中均没有考虑灌溉对地下水溶质运移的影响，有加剧污染扩散的风险，亟需构建水量水质耦合模拟优化模型来进行灌区水资源优化配置。

因此本文考虑灌溉收益、作物水分生产率、地下水位变幅和污染扩散四个目标，构建了紧密耦合的综合考虑水量水质的多目标模拟优化模型，并将NSGA-Ⅱ与NSGA-Ⅲ的实际应用效果进行定量对比分析，最后采用多标准分析方法TOPSIS确定推荐的灌溉方案。

1 研究理论框架

1.1 基于FloPy的地下水数值模拟模型

地下水数值模拟模型为精确定量刻画地下水水流和溶质运移提供了强有力的技术手段[29-31]，本文采用二维潜水含水层非稳定地下水流方程来描述研究区的地下水运动过程：

式中：K为含水层的渗流系数；h为区内各点潜水水位；b为含水层底板高程；μ为给水度；t为计算时间；ω为含水层的源汇项；Ω 为研究区域；Γ1、Γ2为第一、第二类边界；q(x,y,t)为流量边界；x,y为横纵坐标。

考虑地下水的对流、弥散、源汇项等，采用单一化学组分的三维迁移偏转方程来描述污染源的扩散情况：

式中：R为延迟因子；C为溶解浓度；为吸附浓度；xi为坐标，xi=x,y,z；qi为x,y,z三个方向上的达西速度分量，qi=qx,qy,qz；Dij为二阶弥散系数张量；qs为源汇处单位体积含水层的流量；Cs为源汇的浓度；λ1、λ2为溶解相、吸附相的反应速率常数；θ为孔隙度；ρb为孔隙介质的体积密度。

为便于与多目标遗传算法进行耦合，本文采用Flopy 来构建MODFLOW 和MT3D 模型。FloPy 是由基于python 语言开发的第三方库[32]，可以实现MODFLOW、MT3D、SEAWAT 等模型的读写、构建、运行以及结果显示，并支持最新的MODFLOW6 版本，在FloPy 构建、运行地下水数值模拟模型时可以调用支持并行计算的mf2k_h5_dbl_parallel. exe、mf2k_h5_parallel.exe等程序来实现地下水数值模拟模型的并行运算。

1.2 基于Pymoo的多目标模拟优化框架

Pymoo 是由Blank 等[33]人基于Python 语言开发的一个能实现多目标优化功能的第三方库，它可以采用不同的算法解决多目标优化问题，并提供丰富的可视化结果。其中遗传算法是用来分析和解决多目标优化问题的一种进化算法，其核心就是协调各个目标函数之间的关系，找出使得各个目标函数都尽可能达到比较大的（或者比较小的）函数值的最优解集。在多目标优化的遗传算法中，由Deb 等[34]提出的NSGA-Ⅱ是影响最大和应用范围最广的多目标遗传算法之一。2014年Deb等[35]又改进后提出NSGA-Ⅲ，也就是根据NSGA-Ⅱ框架提出了一种基于参考点的多目标进化算法，增加了非优势群体成员的权重。NSGA 适用于许多具有3到15个目标的复杂优化问题。Pymoo 中采用超立方体积指标值（Hypervolume，简称Hv）作为评价指标，Hv值收敛时可认为求解所得非劣解集已经趋于合理。

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

NSGA 计算得到的非劣解集包括多个优化方案，生产实践中需要从中选出最终的推荐方案。本文采用一种多标准决策方法TOPSIS(Technique for Order Preference by Similarity to Ideal Solution)，该方法已在水资源管理领域被广泛采用[36,37]。利用TOPSIS，每一方案都被赋予一个相对分数。然后，将各备选方案按其相对分数进行排序，并选择排名最高的方案作为推荐方案。

TOPSIS 方法的主要步骤包括[38]：①构建判断矩阵；②对矩阵进行归一化处理，得到归一化矩阵；③根据熵的定义确定评价指标的熵；④评价指标的熵权，计算各指标权重集；⑤确定理想解与负理想解；⑥计算各方案与理想解的相对接近程度：

式中：S+i为各方案与理想解的接近程度；S-i为各方案与负理想解的接近程度；Ci∈[0,1]，且该方案Ci值越大越好。

1.4 总体框架

本文通过构建多目标模拟优化模型的方法来进行灌区水资源优化配置，其中的决策变量有4 个，分别为3 个灌季的灌溉用水量和灌溉水源的渠井用水比例。采用经验公式计算不同灌水条件下的作物产量，调用FloPy 运行MODFLOW 和MT3D 模型来计算抽水灌溉对地下水位和污染物运移的影响，采用NSGA 算法来求解所构建的多目标模拟优化模型，在得到非劣解集后采用TOPSIS 方法来确定最终的推荐方案。本文的总体框架如图1所示。

图1 多目标模拟优化模型框架Fig.1 The framework of multi-objective simulation optimization model

2 典型灌区案例概化及模拟优化模型的构建

2.1 典型灌区案例概化

本文根据干旱半干旱地区牧草灌溉的野外实地调查概化出典型案例。灌区地势平坦，渠井等灌溉设施完善，同一地块可根据渠道来水情况选择采用渠水或井水进行灌溉。如图2所示，该区域为由一口机井控制的灌溉范围，其面积约为36 hm2，将其概化为一个正方形区域，所开采的地下水为一厚度为35 m 的均质各向同性潜水含水层。含水层顶板高程为635 m，底板高程为600 m，潜水水位631 m。该含水层上层为弱透水层，可接受降水和灌溉水的入渗补给，区域内的作物为苜蓿，每年按照灌溉制度和灌溉习惯分为3 个灌季进行灌溉。灌溉水源来自衬砌完好的渠水（地表水）和机井抽取的地下水。机井位于区域的中心位置。在典型灌区内设置一面状污染源，另外结合垃圾填埋场实际污染物类型特征，选取硫化物作为本次模拟分析的污染物，按照地下水质量标准（GB/T 14848--2017）中Ⅲ类水（主要用于生活和农业用水）水质要求为标准（其中硫化物浓度不大于0.02 mg/L）判断灌区是否受到污染。并结合灌溉抽水的不同情况来探究污染扩散随之发生的变化。本案例将灌溉用水中的渠井用水比例和3个灌季的灌溉用水量作为决策变量，分别根据平水年和枯水年的实测降水资料和灌溉情况考虑作物产量、污染物扩散规模等目标构建水量水质耦合的多目标模拟优化模型来进行多目标水资源优化配置。

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

2.2 地下水数值模拟模型构建

根据典型案例概化的水文地质概念模型如图3所示。潜水含水层水头西高东低。则模拟区南北可概化为由流线组成的零流量二类边界，东西为由等水位线组成的一类水位边界。在FloPy 中编制代码来设定模型的初始条件、边界条件和源汇项[24]（降水入渗补给、地表水灌溉入渗补给、井灌回归补给、机井抽水），并输入渗透系数、弹性给水度、污染源位置和污染物浓度等信息，分别构建非稳定流MODFLOW 水流模拟模型和MT3D溶质运移模拟模型。为保证模拟精度，反映降水和抽水灌溉对地下水的影响，模型的应力期设置为1 d，模拟时段为一个日历年，模拟区域横向和纵向均剖分为60 等份，共计3 600个单元格（图3）。构建好的模型可以定量模拟降水入渗、渠水灌溉、抽水灌溉等对地下水位和和污染物扩散运移的影响。

图3 典型模拟区水文地质概念模型概化图Fig.3 Schematic diagram of hydrogeological conceptual model in typical simulation area

2.3 水量水质耦合的多目标模拟优化模型的构建

模拟优化模型中考虑的4 个决策变量设定为（Q1,Q2,Q3,x）,其中Qi（i=1,2,3）为各个灌季的灌溉水量（m³/hm2），x为灌溉用水中地下水所占的比例。灌溉用水量根据灌溉制度在灌溉时间内进行逐日平均分配并考虑入渗补给后作为源汇项输入地下水数值模拟模型中。

模拟优化共考虑4个目标，分别为灌溉效益、作物水分生产率、灌溉对地下水动态的影响和污染物扩散范围。约束条件为灌溉抽水量的上下限。其中4个目标函数以及约束条件表述如下。

（1）效益最大化目标：

其中产量计算公式采用Li等[39]于2017年实验所得产量公式：

式中：I1、I2、I3为第一次、第二次、第三次的实际灌溉量，mm；A为灌区面积，hm2；Ys为苜蓿总产量，kg；P、Ps、Pg分别为苜蓿单价，元/kg、地表水单价，元/m³、地下水用电单价，元/kWh；Y1、Y2、Y3分别为第一次、第二次、第三次刈割的干苜蓿产量，kg/hm2。

（2）作物水分生产率目标[40]：

式中：θ为作物水分生产率，kg/m³；λ为复种指数；G为粮食平均单产，kg/hm2；IR为平均灌溉用水量，m³/hm2；Pe为粮食作物生长期内实际补充到作物根分布层可被其利用的有效降水量，mm。

（3）地下水位变化目标，统计日尺度的地下水水位降深之和。

调用MODFLOW 地下水数值模拟组件计算得到模拟区单元格水位矩阵。再由水位矩阵中的数据计算求得每个单元格的水位降深。

式中：k为各网格编号，k=1，2，…，3 600；t为决策时段长度，t=1，2，…，365，d。hkt为地下水模型计算所得逐日的各单元格的水位数据，m。

（4）污染的扩散面积目标：

式中：S为总扩散面积，hm2；A为灌区总面积，hm2；s为MT3D 模型根据地下水运动过程计算所得的污染物扩散的单元格数量，会随着不同的地下水抽水过程（即Q1,Q2,Q3的大小以及地下水占比x）而变化；n为总单元格数。

（5）约束条件。设置不同水平年每个灌季灌水量上下限以及总灌水量约束：

平水年：

枯水年：

式中：Q1、Q2、Q3分别为3个灌季的灌水量，m³/hm2

3 模型求解

根据已有的4个目标，结合地下水数值模型，利用多目标优化工具Pymoo建立以遗传算法框架为基础的多目标管理优化模型对模型进行求解。经过多次反复试验，在既保证优化结果的收敛性，又兼顾模型运行一次所需时间的前提下，设定种群数为50，遗传代数为50 代，交叉率0.8，变异率0.05。本次运行求解所用电脑CPU 为Intel(R) Core(TM) i5-7200 U，频率为2.50 GHz，运行内存为8 GB，进行一次完整的求解过程（50×50）大约需要8 h。

4 优化结果分析

4.1 NSGA-Ⅱ和NSGA-Ⅲ应用效果对比分析

分别将目标函数灌区效益、作物水分生产率、地下水位变化3 个模型计算所得的目标函数值设置为x、y、z轴，污染扩散面积用颜色标记，绘制枯水年和平水年的非劣解集四维空间分布图（图4和图5）。从模拟优化的结果图中可以初步看出，NSGA-Ⅱ在模拟优化过程中所能保存的非劣解集数量明显多于NSGA-Ⅲ。反复验证分析后发现，NSGA-Ⅱ是按照所设定的种群数量来保存非劣解集的，即本次所设定的每代50个种群，模拟优化所留下的非劣解集的数量仍然是50；但由于NSGA-Ⅲ增加了非优势群体成员的权重，种群之间的竞争淘汰会更加激烈，这使得NSGA-Ⅲ算法下所能保存的非劣解集的数量只有20个左右，明显少于NSGA-Ⅱ。

图4 枯水年优化结果空间分布图Fig.4 Spatial distribution of optimum results of dry year

图5 平水年优化结果空间分布图Fig.5 Spatial distribution of optimum results of normal year

按照效益最大，作物水分生产率最高，地下水位变化最小分析，最优的非劣解集应该出现在空间的左下角，其中颜色越偏向红色说明污染扩散越严重，越偏向绿色说明污染控制的越好。总体来说非劣解集的空间分布符合4个目标的权衡关系。再分别对比枯水年与平水年条件下两种遗传算法保留的非劣解集可以看出，其空间分布基本一致，各目标函数变换范围也趋于相同，但是如之前所分析，NSGA-Ⅲ会保留更少的精英方案，非劣解集中地下水位变化和污染扩散面积这两个目标函数值明显偏大的方案显著减少。

两种年份下的NSGA-Ⅱ和NSGA-Ⅲ的Hv值分布图如图6和图7 所示，其中在进行50 代的求解过程中NSGA-Ⅱ对于结果的收敛效果明显不如NSGA-Ⅲ，其中枯水年优化过程中NSGA-Ⅱ前30代繁衍hv值波动较大，而NSGA-Ⅲ在繁衍不到10 代的情况下即已趋于稳定；平水年优化过程中NSGA-Ⅱ甚至出现未收敛的情况，而NSGA-Ⅲ依旧可以在繁衍30～40 代后基本趋于稳定收敛。由此可以看出在本次模拟优化模型的求解效率上NSGA-Ⅱ低于NSGA-Ⅲ，不能实现快速收敛。

图6 枯水年优化结果Hv值Fig.6 Hypervolume value of optimization results in dry year (NSGA -Ⅱ(left) NSGA -Ⅲ(right))

图7 平水年优化结果Hv值Fig.7 Hypervolume value of optimization results in normal year

综合以上来看，NSGA-Ⅲ保存的非劣解少于NSGA-Ⅱ，但是NSGA-Ⅲ可以更快的使模拟结果收敛而寻找最优解，在本次模拟种群数目为50 进行50 代进化计算所得的优化方案中，其效率是明显优于NSGA-Ⅱ的。所以后续的方案筛选将主要在枯水年和平水年的NSGA-Ⅲ的优化结果中进行。

4.2 多标准分析法选择推荐方案

由NSGA-Ⅲ计算得到的非劣解集中不同水平年的灌溉方案分别对应4个目标函数值（表1和表2），其中平水年保留的非劣解集包含20 个方案，枯水年包含18 个方案。各方案目标函数值对比如图8和图9所示。

表2 NSGA-Ⅲ优化枯水年方案结果Tab.2 Dry year results of NSGA-Ⅲ optimized solution

图8 平水年非劣解集目标函数值Fig.8 Objective function value of Pareto front in normal year

图9 枯水年非劣解集目标函数值Fig.9 Objective function value of Pareto front in dry year

为更好的看出非劣解集中每个方案权衡4个不同目标的表现，绘制平水年各方案的百分比堆积图如图8和图9所示。

从各方案对比的图表中可以看出，每个方案对不同目标的倾向不同，难以直接评价某一方案的优劣程度。因此本文采用TOPSIS 分析评价方法，利用各方案中目标函数值的最值自动生成默认的理想解与负理想解，然后由公式（4）依据各方案与理想解和负理想解的距离进行综合打分，依据打分结果分别选出了枯水年和平水年的推荐方案如表3所示。

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

结合前文所提到的问题，重点分析两个推荐方案灌溉模式下的地下水水位变化以及污染扩散情况，利用地下水数值模拟模型模拟其一整年的地下水位变化与污染扩散的最终结果并作图比较分析。

灌溉方案A、B 的抽水井和污染源处地下水位变化如图10和图11 所示，抽水井处水位随灌溉进行变化明显，污染源处水位也会随灌溉的进行发生一定的变化，由此可以看出灌溉导致的地下水位变化抽水井处向四周逐渐减弱，但其波动也已经影响到污染源处的地下水。其中不同水平年条件下的灌溉方案A、B 灌溉时典型灌区抽水井和污染源处的地下水位逐日变化趋势基本一致，但是枯水年的灌溉方案B地下水位变化幅度相对更大。为更进一步的体现灌溉对灌区地下水水位变化的影响，再分别选取其中地下水位波动最大的一日，作出整个灌区的地下水水位图（图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最大单日灌溉量时地下水位Fig.12 The groundwater level at maximum daily irrigation volume in solution A and solution B

而在实际模拟一整年的灌溉抽水后污染物的总体扩散情况如图13 所示，由于本案例的典型灌区采用的节水灌溉模式中的喷灌模式，总体灌溉水量不大，又引进了地表水水源极大的减轻了地下水灌溉压力，所以最后的污染物扩散情况控制比较好。

图13 方案A和方案B污染扩散情况Fig.13 Pollution diffusion of solution A and solution B

而若不引进地表水水源且一直按照各水平年最高灌溉定额灌溉，其对应的最大单日灌溉抽水时的水位变化和最终的污染物扩散如图14 和图15 所示。其污染物的扩散面积要明显大于前面的两个推荐方案。根据模型进一步计算结果，平水年的推荐灌溉方案A 要比全地下水最大灌溉量灌溉方案的地下水位变化幅度减少了0.55 m（56%），污染扩散面积减小了1.24 hm2（54%），枯水年的推荐灌溉方案B 要分别比全地下水灌溉方案的地下水位变化幅度减少了0.51 m（43%），污染扩散面积减小了1.13 hm2（40%），模型对重点关注的灌区地下水水量水质的优化效果显著。

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

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

5 结论与讨论

基于FloPy 和Pymoo 构建了典型渠井结合灌区水量水质耦合的多目标模拟优化模型。模型采用地下水数值模拟模型MODFLOW 和溶质输移模型MT3D 模拟地下水位变化和污染源扩散过程。分别采用多目标遗传算法NSGA-Ⅱ和NSGA-Ⅲ考虑4个目标（灌区收益最大化、作物水分生产率最大、地下水位变幅最小、污染扩散面积最小）获得了该复杂问题的非劣解集。最终采用TOPSIS 方法获得不同水平年的推荐方案。主要结论如下：①运用基于Python 版本的FloPy 和Pymoo 构建了耦合水量水质的灌区水资源多目标模拟优化模型，优化结果表明该模型可以较好的权衡多个目标之间的竞争关系；②构建的模型可以使用Python 语言命令直接调用MODFLOW 和MT3D 的可执行程序来运行模型，且模型设置目标函数以及修改参数简洁方便，可以直接使用代码循环多次运行模型，并能够基于不同版本的遗传算法进行求解，显著的提高了模型运行调试和问题求解的效率；③对比了不同版本遗传算法（NSGA-Ⅱ、NSGA-Ⅲ）在灌区水资源多目标优化方面的实际应用效果。结果表明，相比NSGA-Ⅱ，NSGA-Ⅲ会因为对非优势种群的兼顾，使得竞争更加激烈，以至于会丢失一部分的种群。但是NSGA-Ⅲ收敛效果更好、效率更高、提供的非劣解集更加稳定；④以NSGA-Ⅲ优化所得非劣解集为基础，利用多标准分析方法（TOPSIS）选择出了不同水平年的推荐灌溉方案。在保证收益和作物水分生产率的前提下，推荐灌溉方案明显的减少了对地下水的影响（水位变化以及污染的扩散），平枯水年分别可以最高减少地下水位变幅56%、43%，最高减少平枯水年的污染扩散面积54%、40%。验证了多目标模拟优化模型在灌区水资源优化配置方面的有效性与可靠性。

本研究从工具以及技术角度探究了多目标遗传算法在灌区进行考虑多重因素的灌区水资源优化配置，后续研究可以与实际情况相结合，扩大研究范围，耦合作物模型并考虑农业农村污染和能源消耗情况，进一步构建符合我国灌区实际情况的耦合作物、水量水质和能源等多目标的模拟优化模型。而针对模拟优化过程中数值模拟模型运行时间较长的情况，可以尝试使用基于机器学习的替代模型，以提高模拟优化模型运行效率。