地下水是重要的饮用水源. 由于人类发展忽略了环境保护，导致地下水污染严重，其中硝酸盐污染问题突出[1]. 过量的硝酸盐会对人类健康造成严重危害，例如蓝婴病[2]. 因此，治理地下水硝酸盐刻不容缓[3].
目前，国际上地下水硝酸盐污染治理技术分为原位处理技术和抽出处理技术. 抽出处理技术仍然是最常用的一种. 抽出处理技术的核心在于如何有效地把污染的地下水从含水层中抽出来，通常考虑的因素包括抽水井的位置、 抽水井数、 抽水速率、 安装井和抽水处理成本，需要找到一个合适的方法得到最优解，从而使治理成本最小化[4]. 因此，其关键是多目标优化问题. 目前，国外已有许多模拟优化技术运用到抽出处理技术设计[5]. 传统的优化方法包括线性规划、 非线性规划、 动态规划、 二次规划和整数规划[6]. 较新的优化方法包括遗传算法(GA)[7, 8, 9, 10, 11]和模拟退火法(SA)[12, 13, 14, 15, 16]. 这些优化方法可以耦合到地下水模拟模型中.
抽出处理在优化过程中容易陷入局部优化. GA和SA采用了多种全局优化方法，能够实现全局优化. Dougherty等[17]首次将模拟退火法(SA)运用到地下水修复最优设计，在其设计中包含井的安装费用和抽水费用. McKinney等[18]运用遗传算法解决了最大抽水量和最低抽水成本等地下水优化问题. Yoon等[19]运用实数编码遗传算法来确定最优抽水速率和最佳井的位置，最大限度减小了地下水修复系统成本. Erickson等[20]利用小生境技术遗传算法解决了抽出处理地下水修复问题，实现了修复成本和剩余污染物浓度最小化. 吴剑锋等[21]开发改进了一种小生境Pareto遗传算法(INPGA)，以美国麻省军事保护区为实例，通过建立研究区复杂地下水污染治理的多目标优化管理模型，并添加消息传递接口(MPI)的并行计算和个体适应度，提高了运行库计算速度.
本文以北京某场地生活垃圾填埋场硝酸盐污染地下水为研究对象，采用Visual MODFLOW 4.6中Modular Groundwater Optimizer (MGO)[22]，将水流模拟软件 MODFLOW、 溶质运移MT3DMS与遗传算法和模拟退火法相结合，尝试解决地下水硝酸盐污染抽出处理优化模拟，优化井的位置、 井的抽水速率及最小化治理成本. 其执行程序通过FORTRAN自动编码器来汇编，属于内置模块，在优化过程中采用DOS界面运行.
1 地下水硝酸盐抽出处理优化系统设计
1.1 遗传算法
遗传算法(Genetic Algorithm)是模拟达尔文生物进化论的自然选择和遗传学机理的生物进化过程的计算模型，是一种通过模拟自然进化过程搜索最优解的方法，它最初由美国Michigan大学Holland教授于1975年首先提出. 采用概率化的寻优方法，能自动获取和指导优化的搜索空间，自适应地调整搜索方向，不需要确定的规则. 遗传算法的核心是选择、 交叉、 变异. 更多遗传算法的描述参见文献[23, 24]. 从任一初始种群出发，通过选择、 交叉、 变异操作，产生一群更适应环境的个体，检查其结果是否符合优化的要求，如果不符合，重复其上述过程，直到优化出满意的结果.
1.2 模拟退火法
模拟退火法，源于统计热力物理学，它模拟熔融状态下物体逐渐冷却达到结晶状态的物理过程. 材料中的原子原来会停留在使内能有局部最小值的位置，加热使能量变大，原子会离开原来位置，而随机在其他位置中移动. 退火冷却时速度较慢，使得原子有较多可能可以找到内能比原先更低的位置. 不依赖初始模型的选择，能寻找全局最小点而不陷入局部极小[19]. 模拟退火算法的操作步骤如下.
步骤1：给定初始温度T0解X0算f0=f(x0); 步骤2：随机产生扰动Δx,计算f1=f(x0+Δx)和Δf=f1-f0; 步骤3：若Δf≤0则转步骤5，否则产生随机数R∈[0,1]; 步骤4：若p=exp(-Δf/T)≥R，则用x′=x0+Δx代替x0， f1代替f0; 步骤5：检查M过程是否结束，若未结束则转步骤2; 步骤6：降低温度参数T; 步骤7：判断终止准则是否满足，若满足，则算法终止，否则，转步骤2.
1.3 地下水多目标优化模型
地下水优化模型通常包含两组变量：决策变量和状态变量. 在优化过程中的目的是识别这些决策变量的最佳组合. 状态变量依赖地下水流方程的水头和地下水溶质运移方程的浓度. MODFLOW水流方程如下：
式中，Kx、 Ky、 Kz分别为渗透系数在x、 x、 y方向的分量; h表示水头(m); qs表示单位时间从单位体积含水层流入或流出的水量(d-1); Ss表示贮水率(m-1); t表示时间(d). MT3DMS地下水溶质运移方程：
式中,θ表示含水层的孔隙度(无量纲); cK表示溶质组分K的浓度(mg ·L-1); D ij是水动力弥散系数张量(m2 ·d-1); vi是孔隙中实际水流速度(m ·d-1); qs表示单位时间从单位体积含水层流入或流出的水量(d-1); cKs是源汇项溶质组分K的浓度(mg ·L-1); Rn表示化学反应项总和[mg ·(L ·d)-1].
优化设计过程中，目标函数包含井的安装费用、 抽水费用、 抽水地下水硝酸盐治理费用等. 其数学表达式如下.
最小化：
约束条件：
式中,J表示目标函数，a1表示安装一口井的总费用，W表示研究区备选井的数量，N是总的优选井数量，即从备选井中优化出来井的数量，a2是抽出和处理单位体积的水的费用[4]，Qi单位时间抽水的体积(m3 ·d-1)，Qmin和Qmax分别是最小抽水和最大抽水量，Δti是井i持续抽水时间(d)，hi是井i的水头(m)，hmin和hmax分别是最小水头和最大水头(m); ci是井i的溶质的浓度(mg ·L-1)，c*表示目标治理浓度(mg ·L-1). 上述定义多目标优化方案也可以通过试错法解决. 虽然试错法简单、 易用，但测试和检验烦琐，且不能保证组合的最优方案，容易产生局部优化，往往很难找到最优的抽出处理方案. 所以在满足所有约束条件下，应用遗传算法和模拟退火法比试错法更有效解决井的位置及组合问题[22].
图 2 研究区三维地质剖面
2 材料与方法
2.1 研究区水文地质条件
研究区位于北京某场地，所处区域地势平坦. 场地东侧为养殖场，西南侧为当地农民耕作农田，农村垃圾和畜禽养殖废水均排放到池塘(图 1).
图 1 研究区初始流场示意
依据场地调查资料，利用CTECH对研究区构建三维地质剖面图，第一层主要含水层为非承压含水层，地下水标高16~23 m，厚度约6.92 m，以粉细砂为主，局部夹有砂质粉土、 黏质粉土和粉质黏土(图 2). 研究区第一含水层下面主要为砂质粉土、 黏质粉土，渗透性很差，能够起到较好的隔水作用.
在研究区进行单孔抽水试验. 通过Aquifer Test软件分析，结果见表 1. 研究区含水层的平均渗透系数为3.37 m ·d-1.
表 1 不同求参方法中渗透系数和导水系数结果表
2.2 研究区地下水污染状况
利用Geoprobe钻机对生活垃圾填埋场地下水及池塘污水进行取样调查. 现场测试数据：池塘地表水的NH+4-N高达46.8 mg ·L-1，同时地表 50 cm渗坑NO-3-N高达65 mg ·L-1. 研究区地下水硝酸盐仍处于不断变化的情况，同时场地调查尚不够充分，不足以刻画出实际地下水硝酸盐的污染羽状况. 因此，本研究以农村垃圾和畜禽养殖废水为污染源，硝酸盐入渗浓度为65 mg ·L-1，模拟在天然地下水流动情况下释放5 a形成的硝酸盐污染羽作为抽出处理的初始污染羽. 研究区硝酸盐初始污染羽总量为36.99 kg，见图 3.
2.3 抽出处理备选井方案
为了设计合理的抽出处理治理优化方案，根据模拟生成硝酸盐污染羽形成的范围及通过试错法确定的较为有效的抽水井位置，在研究区布置了18口备选井方案(图 3). 目的是在满足治理成本最小情况下，抽出处理后达到目标治理浓度，同时确定最佳井数及其井的位置和抽水速率.
研究区周围地下水硝酸盐背景浓度约为1 mg ·L-1. 将研究区红色框作为浓度约束区，在100 d抽出处理修复时间内，其硝酸盐目标修复浓度低于10 mg ·L-1.
2.4 模型概化及参数确定
研究区场地范围较小，可以将研究区假设为均质各向异性含水层. 地下水硝酸盐污染浓度主要分布在第一含水层，该层下部主要为砂质粉土、 黏质粉土，渗透系数较差，能够起到较好的隔水作用，故 而将含水层概化成一层. 研究场地为225 m×200m，将其剖分成2.25 m×2 m. 地下水从北东方向向南西方向流(图 1). 研究区的初始水力坡度约为0.01. 西北方向和东南方向为隔水边界. 地下水硝酸盐的阻滞因子设为1. 地下水流模拟设置为稳定流模拟(两种优化方法只适合稳定流[4])，模型输入参数见表 2.
图 3 地下水硝酸盐初始污染羽分布示意和备选井位置
表 2 模型输入参数
表 3列出了安装一口井的总费用、 抽出和处理单位体积水的费用[4]. 最小与最大水头值、 最小与最大抽水量及硝酸盐目标治理浓度等约束条件.
表 3 费用系数和约束变量
3 优化结果分析与讨论
3.1 抽出处理优化
本研究中GA设定的相关参数如下：种群数目为100，最大遗传代数为100，每口井的抽水速率离散化区间数为32，交叉概率为0.90，变异概率为0.05. SA设定的相关参数如下：初始温度为1，抽样次数10000，降温率为0.9，最大迭代次数为100.
两种优化方法都要求浓度约束区治理100 d后的地下水硝酸盐(以氮计)达到10 mg ·L-1. GA优化结果：8号井优化抽水速率为155 m3 ·d-1，14号井优化抽水速率为10 m3 ·d-1; SA优化结果：8号井优化抽水速率为82 m3 ·d-1，14号井优化抽水速率为39 m3 ·d-1. 其他井都为零(表 4). GA在8号井优化的抽水速率比SA在8号井的要大，在14号井优化的抽水速率比SA在14号优化的速率要小. 但SA总抽水速率(121 m3 ·d-1)比GA(165 m3 ·d-1)小.
表 4 优化井位置与抽水速率
基于GA和SA优化的抽水速率，100 d后硝酸盐总去除量分别为28.45 kg和31.42 kg，其总量去除率分别76.89%和84.92%，见图 4. 优化过程中发现，总抽水速率和硝酸盐总去除量存在近似线性增加的趋势.
GA、 SA同时优化了井的位置和抽水量. 图 4可以看出研究区中18口备选井方案，两种优化方法最优的两口井都是8号井和14号井. 说明在该场地抽出处理中，只需要2口井就可以满足约束条件，到达治理目标，比初始方案减少16口井，节省了井的安装成本. 同时说明最优井的位置在污染羽的中轴线上中游和中下游抽水布置最为合理.
图 4 GA和SA抽水治理100 d后地下水硝酸盐分布示意
3.2 GA与SA对比分析
GA、 SA不仅可以优化井的位置、 抽水速率及井数，而且还可以优化系统治理成本(井安装成本与抽出处理成本之和)，见表 5. GA和SA优化井数一样，其安装成本也就一样，但SA优化总的抽水量比GA要低，所以SA优化的系统成本要比GA要低，说明SA在多目标地下水硝酸盐修复优化管理要比GA体现出一定的优越性，体现在节省系统治理成本.
图 5为GA与SA遗传代数与目标函数优化收敛对比. GA在40代后优化目标函数趋于收敛，而SA在20代后趋于稳定. GA在前20代目标函数优化波动性比SA要大. 两者最终都趋于稳定. 两种优化方法识别决策变量的最佳组合. 在最小化治理成本下，全局优化井数、 井的抽水量、 安装成本、 抽水处理成本，使得优化结果可信. SA的系统治理成本要比GA少3338.45元，节省成本约6.8%，见表 5具体参见 污水处理技术资料或污水技术资料更多相关技术文档。
图 5 GA、 SA目标函数优化收敛对比
表 5 GA、 SA优化系统治理成本/元
4 结论
(1)本研究依据北京某场地调查资料，利用遗传算法(GA)和模拟退火法(SA)对地下水硝酸盐污染羽治理区中布置18口备选井方案进行优化. 为达到100 d将研究区硝酸盐浓度降低到10 mg ·L-1的目标，GA和SA优化的井都为8号井和14号井，优化总抽水速率分别为165 m3 ·d-1和121 m3 ·d-1，同时其硝酸盐总量去除率分别为76.89%和84.92%.
(2)优化结果说明，最佳井的位置应位于硝酸盐污染羽中游及中下游的中轴线上，且中游抽水速率比下游要大. 两种优化算法对比表明模拟退火法(SA)优化系统治理成本比遗传算法(GA)节省6.8%，且波动性小，更收敛.
(3)模拟退火法(SA)和遗传算法(GA)耦合到地下水流模拟程序MODFLOW和溶质运移程序MT3DMS中，可以设计合理的抽出处理方案和符合管理目标修复系统，为地下水污染现场修复提供科学的指导.
致谢： 本研究场地基础信息调查由中国环境科学研究院、 北京建工环境修复有限责任公司、 轻工环保研究所等共同完成. 此外，南京大学吴剑锋老师对遗传算法优化方法工作给予指导!在此一并表示感谢.(来源及作者：中国地质大学(北京)水资源与环境学院, 水资源与环境工程北京市重点实验室 姜烈、何江涛、姜永海、刘菲)