基于自适应权重粒子群优化算法的地下水污染溯源辨识

2022-12-26高琬玉卢文喜潘紫东白玉堃

中国农村水利水电 2022年12期

高琬玉，卢文喜，潘紫东，白玉堃

（1.吉林大学地下水与资源环境教育部重点实验室，吉林长春 130012；2.吉林大学新能源与环境学院，吉林长春 130012）

0 引言

地下水资源具有分布广泛、不易被污染的特点，是人类生活生产可利用水资源的重要组成部分，其质量好坏是影响环境和生态的重要因素。由于地下水污染［1］具有发生的隐蔽性、发现的滞后性、修复难度大、修复费用高等特点［2］，获取有效的污染源相关信息是地下水污染修复方案设计的前提条件，因此进行地下水污染溯源辨识具有重要的实际意义。

地下水污染溯源辨识是指利用有限的地下水污染现场实际监测数据（水位和污染物浓度），以及现场调查和专业知识等辅助信息，对地下水污染数值模拟模型进行反演求解，识别确定污染源的个数、位置及释放历史（污染质在各时段的释放强度）等相关信息［3］。在数学上，地下水污染溯源辨识问题属于数理方程反问题，具有非线性和不适定性的特点，求解难度较大［4］。

有关地下水污染溯源辨识的理论和方法，至目前已有随机理论与地质统计法、数理方程反演法、同位素法、地球化学指纹法和地球物理探测方法等多种理论方法。总体来讲，目前应用的理论和方法中，数理方程反演法中的模拟-优化方法被较为广泛地应用于辨析确定地下水污染源的分布范围和释放历史等多方面信息、Singh and Datta［5］采用遗传算法对地下水污染源的位置及污染物质释放强度进行同步地溯源辨析。Prakash［6］提出了一种基于克里格方法的监测井网优化设计方法，并在此基础上采用遗传算法对地下水污染源释放污染物质的强度进行溯源辨析。Jha and Datta［7］运用模拟-优化方法溯源辨析地下水污染源的释放历史，并采用模拟退火算法求解地下水污染溯源辨识的优化模型。Guneshwor et al［8］运用基于粒子群优化算法的无网格流动优化模型溯源辨析地下水污染源的释放历史。范小平和李功胜［9］运用一种改进的遗传算法，溯源辨析山东省某区域地下水中硫酸盐的年入渗强度。黄林显等［10］应用模拟-优化方法溯源辨析污染源的位置和污染物质的释放强度，并利用复合进化算法求解优化模型。随着智能优化算法的出现，国内外研究学者将其应用于解决地下水污染溯源辨识问题。然而，传统的智能优化算法存在的早熟收敛问题会导致对反演问题的求解精度不高，针对该问题，国内外研究学者多关注于对算法跳出局部最优能力的提升。Koupaei［11］将混沌映射和黄金分割搜索算法相结合，提高了黄金分割搜索算法的局部搜索能力和快速全局收敛能力。吕石磊［12］提出的基于自适应步长的改进蝙蝠算法有效地避免了过早陷入局部最优和求解精度低的问题。

目前地下水污染溯源识别研究的现状和发展趋势具有如下特点：①研究内容多为对单一污染源释放历史识别［13］，然而在实际场地中，可能存在需要同时对多个污染源释放历史识别的情况。②现阶段研究者多运用传统智能优化算法求解优化模型进行地下水污染源反演识别研究［14］，传统智能优化算法的求解精度对初始点的依赖较大，易陷入局部最优。

针对上述问题，考虑了多个污染源的情况，对多个污染源的释放历史及场地渗透系数联合识别，应用自适应权重粒子群算法求解优化模型，对算法的迭代终止条件进行改进，对自适应权重优化算法收敛结果影响较大的参数进行修改，大幅度提高了算法的收敛速度和计算精度。

1 研究方法

1.1 模拟-优化方法

在模拟-优化方法中，模拟是指地下水溶质运移模拟模型，用来描述地下水系统输入（地下水污染源特征、场地水文地质参数等）与输出（水位及污染物浓度）的激励-响应关系；优化指的是优化模型，其目的是寻找待识别变量的真实值，即优化模型的最优解［15］。

文章以待识别变量为决策变量、以模拟模型的污染物浓度输出值与污染物浓度实测差值最小为目标函数建立优化模型，并将地下水溶质运移模拟模型作为待识别变量的等式约束条件嵌入优化模型中，然后运用群智能优化算法求解优化模型［14］。

1.2 拉丁超立方方法

拉丁超立方方法是一种基于分层采样的方法［16］，能够反映随机变量的整体分布。与其它抽样方法相比［17］，它具有效率高、样本覆盖程度好、样本更具有代表性等优点［18］。对于一维输入变量，拉丁超立方抽样原理为：根据需要，确定抽样数目n，将变量服从的概率密度函数等间隔的分为n个子区间，然后在每个子区间内随抽取一个值，最后将抽取的n个值进行随机排列即完成一次抽样［19］。

根据前人研究，渗透系数对模拟模型的结果影响较大［20］，因此，文章在识别污染源释放历史的同时识别渗透系数。采用拉丁超立方方法，在建立模拟模型的替代模型时，对待识别变量进行抽样，获得训练和检验替代模型的输入样本。

1.3 BP神经网络

在模拟-优化算法中，优化模型的迭代求解过程需反复多次调用模拟模型，这会带来庞大的计算负荷，严重制约了模拟-优化方法在反演识别实际应用中的可行性［21］，替代模型在功能上逼近模拟模型，能够以很小的计算负荷逼近模拟模型的输入-输出响应关系，因此，研究建立模拟模型的替代模型成为近年来研究进展中的前沿问题之一。采用BP 神经网络建立替代模型拟合模拟模型的输出结果。

BP（back propagation）神经网络是1986 年由Rumelhart 和McClelland 为首的科学家提出的概念，是一种根据误差逆向传播使误差极小化进行训练的多层前馈神经网络。它的学习过程由信号的正向传播与误差的反向传播两个过程组成［22］。它的基本思想为：

（1）先计算每一层的状态和激活值，直到最后一层（即信号是前向传播的）；

（2）计算每一层的误差，误差的计算过程是从最后一层向前推进的（即信号是反向传播的）；

（3）更新参数，不断迭代前两个步骤，直到满足停止准则。

其求解过程如下所示：

（1）正向传播。第l(2 ≤l≤L) 层神经元的状态及激活值为：

对于L层感知器，网络的最终输出为a(l)。前馈神经网络中信息的前向传递过程如下：

为了调整权重和偏置使总体误差最小，采用列文伯格-马夸尔特法求总体误差最小值，并求此时所对应的各个神经元的参数（即权重和偏置）。

1.4 自适应权重粒子群优化算法

粒子群优化算法（Particle Swarm Optimization，PSO）是模拟鸟类觅食行为的群智能优化算法［23］。鸟类寻找栖息地的过程与寻找特定问题的解的过程类似。设每个优化问题的解是搜索空间中的一只鸟，把鸟视为空间中一个微粒，每个粒子都有一个由被优化函数所决定的适应度值，还有一个速度决定它们的飞行方向和距离。粒子的运动通过追随当前的最优例子在解空间中搜索最优解。

设n维搜索空间中，粒子i的当前位置为X(i)、当前飞行速度V(i)及所经历的最好位置P(i)（即具有最好适应度值的位置）分别表示为：

对于最小化问题，若f(X)为最小化的目标函数，则微粒i的当前最好位置由下式确定：

设群体中的粒子数为S，群体中所有粒子所经历过的最好位置为Pg(t)，称为全局最好位置，即：

基本粒子群算法粒子i的进化方程可描述为：

式中：vij(t)为粒子第j维第t代的运动速度；C1、C2为加速度常数（学习因子）；r1j、r2j分别为两个相互独立的随机数；Pg(t)为全局最好粒子的位置。式（8）描述了粒子i在搜索空间中以一定的速度飞行，这个速度要根据自身的飞行经历［式（8）中右侧第2项］和同伴的飞行经历［式（8）中右侧第3项］进行动态调整。

带有惯性因子的粒子群优化算法是对于式（8）中的vij(t)项加以惯性权重ω，即：

粒子群优化算法有记忆，所以所有粒子都能保留关于当前最优解的信息［24］。粒子之间有建设性的合作，群中的粒子之间共享信息［25］。然而，简单粒子群优化算法的性能在很大程度上取决于其参数（包括学习因子C1、C2和惯性权重ω），并且经常遇到陷入局部最优从而过早收敛的问题，适当控制全局搜索和局部搜索对于有效地找到最优解至关重要。应用自适应惯性权重因子（AIWF）控制全局搜索求解地下水污染溯源识别问题的优化模型。基本思想如下：如果适应度开始停滞时，粒子群搜索会从邻域模式向全局模式转换，一旦适应度开始下降，则又恢复到邻域模式，以免陷入局部最优。当适应度的停滞次数足够大时，惯性系数开始逐渐变小，从而利于局部搜索［26］。其中，自适应权重公式如下：

式中：ωmin和ωmax是预先给定的最小惯性系数和最大惯性系数，一般取0.4 和0.9；式（12）为第t次迭代时所有粒子的平均适应度；式（13）为第t次迭代时所有粒子的最小适应度。

当粒子已经找到最佳位置后，再增加迭代次数会浪费计算时间，因此，采用了自动退出迭代循环的方法。

（1）初始化最大迭代次数、计数器以及最大计数值，分别取1 200，0，20；

（2）定义“函数变化量容忍度”，一般取非常小的正数，取10-6；

（3）在迭代的过程中，每次计算出来最佳适应度后，都计算该适应度和上一次迭代时最佳适应度的变化量（取绝对值）；

（4）判断这个变化量和“函数变化量容忍度”的相对大小，如果前者小，则计数器加1；否则计数器清0；

（5）不断重复这个过程，有以下两种可能：①此时还没有超过最大迭代次数，计数器的值超过了最大计数值，那么直接跳出迭代循环，搜索结束。②此时已经达到了最大迭代次数，那么直接跳出循环，搜索结束。

2 案例应用

2.1 问题概述

文章参考潘紫东等作者的假想算例，考虑了3 个污染源同时污染的情况，根据9口监测井5次监测数据（45维），同时识别3个污染源五段释放历史以及含水层参数（19维）［27］。如图1所示，研究区概化为二维非均质各向同性的不规则承压含水层（2 000 m×2 500 m），用边长为20 m 的正方形将研究区剖分为100×125个有限差分网格。研究区的渗透系数按介质颗粒分为K1、K2、K3、K44 个分区，南北边界为已知流量边界，东西边界为已知水头边界，含水层和污染源的相关值和范围见表1。模拟时间为5 年，每1 年记为1 个应力期。研究区内存在3 个潜在的污染源，在应力期内向含水层排放污染物。建立9 口监测井（Obs1～Obs9），在应力期内每年对含水层中的污染物浓度进行监测。

表1 含水层与污染源相关参数Tab.1 Parameters of pollution sources and aquifer

图1 研究区概况Fig.1 Overview of the study area

2.2 模拟模型建立

根据研究区水文地质概念模型，建立研究区地下水水流的数学模型：

式中：t为时间变量；K为渗透系数；H为水位高程；B为底板高程；w为源汇项；μ为给水度；Γ1，Γ3为已知流量边界；Γ2，Γ4为已知水头边界；q(x，y，t)和φ(x，y，t)为已知函数；n⇀为边界上某点(x，y)处外法线方向上的单位向量。

在地下水流数学模拟模型的基础上，建立地下水溶质运移数学模型：

式中：n为孔隙度；M为承压含水层的厚度；fM为源汇项，表示单位时间含水层单位面积溶质的增减量；c为地下水溶质浓度；Dxx、Dyy是x、y轴方向上的水动力弥散系数；ux、uy分别为实际平均流速向量u⇀在x、y轴向上的分量；S为研究区；Γ1为已知浓度边界；Γ2和Γ4为已知水动力弥散通量边界；Γ3为已知对流-弥散通量边界；c0(x，y)，c1(x，y，t)，c2(x，y，t)，c3(x，y，t)为已知函数。

利用GMS 软件中的MODFLOW 和MT3DMS 对地下水流动和污染物运移过程进行了计算，模拟模型的计算结果如图2所示。

图2 污染质分布情况示意Fig.2 Distribution of contaminants

由于是一个假想案例，因此，需要人为设定一组待识别变量的值作为真实值，输入模拟模型，运行获得每个时期观测井的污染物浓度，将模拟获得的观测井污染物浓度视为实际监测井浓度。待识别的渗透系数真实值见表2；待识别的真实污染源信息见表3。

表2 渗透系数Tab.2 Hydraulic conductivity

表3 污染源信息Tab.3 Information of pollution sources

2.3 替代模型的建立

利用拉丁超立方方法对19 个输入变量（4 个渗透系数和3个污染源各5 个时段的释放强度）进行抽样，抽取400 组分布较均匀的样本作为输入值。将400 组样本输入值分别代入GMS软件进行计算，计算得到400组9个观测点5年的污染物浓度观测数据作为样本输出值。

将400 组样本随机选取280 组作为样本训练，60 组作为验证集，60 组作为测试集。运用列文伯格-马夸尔特方法训练BP神经网络。

训练集的拟合情况如图3（a）所示，验证集拟合情况如图3（b）所示，测试集的拟合情况如图3（c）所示，整体的拟合情况如图3（d）所示。

由图3的数据可以看出，由BP神经网络建立的替代模型所求得的输出与模拟模型计算所求得的输出，二者之间拟合精度较高，将真实的渗透系数和污染源释放强度等输入代入到替代模型中，将求得的输出与观测数据进行对比，其拟合程度如图4所示。

图3 训练集、验证集、测试集及整体的拟合情况Fig.3 Regression of training，validation，test and all

图4 替代模型与模拟模型输出对比Fig.4 Comparison of output between alternative model and simulation model

因此可以用BP 神经网络建立的替代模型代替地下水溶质运移模拟模型进行运算。

2.4 优化模型建立及求解

优化模型由3 个部分组成：①待求的地下水污染源释放历史及场地的渗透系数作为决策变量；②各个监测井污染质浓度实际监测值与模拟计算值之差的绝对值极小化作为目标函数；③决策变量满足地下水溶质运移规律和合理地范围作为约束条件（替代模型作为等式约束条件）。

分别运用传统粒子群算法和自适应权重粒子群优化算法对建立的优化模型进行求解，运用自适应权重粒子群算法时采用了自动退出迭代循环的方法，使其能够在寻找到最优解时自动提前跳出迭代。研究最终跳出迭代循环时，迭代次数为707代，运行时间为2 分45 秒。传统粒子群算法不能提前跳出迭代，且运行一次迭代的时间为1.5 秒，为了对比两种算法的运行结果，同样运行了707 代，耗时5 分54.5 秒。目标函数值收敛曲线如图5 所示，同等条件下，传统粒子群算法易陷入局部最优，耗时更长，而自适应权重粒子群算法的求解精度更高，耗时较短。

图5 优化识别过程中目标函数值收敛曲线Fig.5 Convergence curve of objective function value

函数中的参数初始值对算法的寻优能力和计算时间有着很大的影响，因此本文运用控制变量法对函数中的参数逐一进行测试，研究发现，当粒子的数量在200 以下时，目标函数收敛曲线下降慢，计算结果精度低，容易陷入局部最优解，当粒子数量在500以上时，目标函数收敛曲线基本不变，但运算速度随着粒子的数增加而明显减慢；邻域内粒子的比例越大，目标函数收敛曲线下降速度越慢，邻域内粒子的比例越小，目标函数收敛曲线下降速度越快；惯性权重范围、个体学习因子和社会学习因子之间相互影响。

经过反复试验，最终确定粒子数量为500，邻域内粒子的比例为0.5，设置个体学习因子和社会学习因子为0.8，惯性权重范围为0.4～1.6；最大停滞迭代数为20。运用自适应权重粒子群算法求解优化模型最终得到渗透系数和污染源释放历史的反演结果。

由图5、表4、5 和表6 可知，应用自适应权重优化算法求解优化模型，能够以较快的速度搜索到全局最优，渗透系数反演结果的相对误差均小于5%，污染源信息反演结果的相对误差绝大多数小于5%，说明自适应权重粒子群优化算法的求解精度较高。