基于仿生原理的粒子群算法求解江安校区水动力弥散系数_梅杰.pdf

第 47 卷 第 6 期 煤田地质与勘探 Vol. 47 No.6 2019 年 12 月 COAL GEOLOGY 2. China Electrical Construction Group Chengdu Engineering Corporation Limited, Chengdu 611130, China Abstract Hydrodynamic dispersion coefficient is an important parameter in the study of solute transport in groundwater. In order to understand the transport law of pollutants in groundwater, the particle swarm optimiza- tionPSO algorithm based on bionic principle is used to solve the hydrodynamic dispersion coefficient of phreatic aquifer under natural flow field in the dispersion test site of Jiang’an Campus of Sichuan University. Compared with the least-square and the standard curve comparison , the results show that the results of the standard curve are subject to subjective influence, and the errors are relatively large. The results of least square fit well with the measured data, but the calculation process is relatively complex. Particle swarm optimizationPSO is a reliable solution , has the highest accuracy, faster calculation, good convergence and stability. Keywords hydrodynamic dispersion coefficient; particle swarm optimization algorithm; groundwater pollution; disper- sion experiment; phreatic aquifer 随着工业发展和人口增长，我国大中城市地下 水存在着不同程度的污染，部分地区地下水污染呈 加重趋势，并有从点状污染向带状和面状污染发展 的趋势[1]。因此，研究污染物在地下水中的运移规 律，从而了解污染物在潜水含水层中的时空分布特 征，对地下水污染防治具有重要研究意义[2-3]。 在多孔介质含水层中，水动力弥散系数为表征 在一定流速下多孔介质对某种污染物质弥散能力的 参数， 它在宏观上反映了多孔介质中地下水流动过程 和空隙结构特征对溶质运移过程的影响[2,4]。水动力 ChaoXing 第 6 期 梅杰等 基于仿生原理的粒子群算法求解江安校区水动力弥散系数 99 弥散系数可通过室内实验或现场弥散试验确定， 以往 研究成果中有不少关于此参数的确定方法。刘兆昌[2] 使用标准曲线对比法成功求解，得到含水层的弥散 系数；李世钰等[5]通过蝙蝠算法求解了内蒙古自治 区呼伦贝尔草原某铜钼矿矿区的含水层弥散系数， 获得了良好的效果；S. Zou 等[6]采用地球物理学的 DPA 法求得了较高精度的弥散系数值。 本文将依据四川大学江安校区天然流场下的弥 散试验数据，运用粒子群算法求解其弥散系数，并 与标准曲线法和最小二乘法的计算结果作比对，比 较 3 种方法的优劣性，为地下水污染防治中弥散系 数的求取提供可靠理论依据。 1 江安校区弥散试验场基本情况 四川大学江安校区地下水弥散试验场位于成都 市双流区四川大学江安校区金工实验厂与水利水电 实验楼之间的绿化带内，下伏地层为粉质黏土夹卵 砾石，地下水属于砂卵石孔隙潜水，含水层厚度约 15 m，地下水埋深为 3.254.81 m。根据前期水文地 质勘察结果，试验场地下水流向为从北向南流动。 试验场布置 2 个抽水孔和 8 个观测孔图 1， 将编号 为 GE4、W4、Q1、Q2、Q3 和 W5 的 6 口井，两两 间距分别为 30、5、10、15、30 m，并依次沿地下 水流向布设。Q4、Q5 和 Q6 按与主流向成 7方向布 设， Q7 与 W4 的距离为 20 m， 按垂直于地下水流向 方向布置。其中 W4 和 W5 为抽水孔，GE4、Q1Q7 为观测孔， 每个孔均可作为示踪剂投源孔和取样孔。 本次试验仅采用抽水井 W4 与观测孔 Q1。 图 1 江安试验场点位布置图 Fig.1 Layout of test sites in Jiang’an Campus 2018 年 3 月 22 日 11 时 30 分，以 W4 井为投 源孔，投放 30 kg 食盐，以距离 W4 井 4.6 m 的 Q1 为取样观测孔，测定电导率；通过不定期在 Q1 观 测孔采取水样分析测定 NaCl 浓度，进行电导率率 定。 电导率与 NaCl 浓度测定至 2018 年 9 月 9 日 14 时 40 分结束，共获得 419 个电导率观测数据，换算 成 NaCl 浓度后的时间变化关系如图 2。 图 2 Cl–浓度随时间的变化关系 Fig.2 The relationship between Cl– concentration and time 2 野外单井一维流动二维弥散模型 在地下水流场中，确定一个二维 xy 平面。假设 投源井的坐标为0,0，在投源井中单位厚度含水层 投放质量为 m 的示踪剂。地下水流线方向为 x 轴， 沿流线方向距离投源井 x 的观测井位置则为x,0。 在不同时间 t 情况下，测得观测井中的示踪剂质量 浓度为 c。得到数学模型[7] 22 22 LT 22 ,,0 , ,00 , d d 0 , , |0 0 xy cccc DDvx yt txxy c x yx y nc x ymt c x y tt                          1 当地下水流速较大时可以忽略分子扩散系数， 同时，假设弥散系数与孔隙平均流速呈线性关系， 则 LL Dv ， TT Dv ，数学模型的解为  22 LT LT , ,exp[ / ] 444π mxvty c x y t vtvtvt n     2 式中 cx,y,t为 t 时刻点x,y处示踪剂的质量浓度， kg/L；m 为单位含水层厚度投放的示踪剂质量， kg/m；n 为岩层孔隙度；v 为地下水的实际流速， m/h；t 为时间从开始投放示踪剂起算，h；αL为 纵向弥散度，m；αT为横向弥散度，m；DL为纵向 弥散系数，m2/h；DT为横向弥散系数，m2/h；x 为 平行地下水流向观测井与投源井之间的距离，m； y 为垂直地下水流向上观测井与投源井之间的距 离，m。 陈崇希等[8]在地下水溶质运移理论及模型 中提出纵向弥散系数通常可取为横向和垂向弥散系 数的 10 倍。且当投源井为0,0点，x 轴为地下水流 向方向，则 y0，简化2式得 ChaoXing 100 煤田地质与勘探 第 47 卷  2 2 L L , ,exp[] 4 4π/10 /xvmt c x y t vt vt n     3 式3中，可根据弥散试验所得质量浓度随时 间变化曲线，采用标准曲线法、最小二乘法[9]等 解析方法求得弥散度和流速。上述方法求解参数 的基本原理均属于曲线拟合，而曲线拟合实质为 最优化算法问题。近 20 年来，在最优化算法领域 兴起了基于仿生学原理的优化算法，如蝙蝠算法、 粒子群算法、布谷鸟算法、蛙跳算法、蚁群算法 等。为丰富求解弥散系数方法，本文尝试采用粒 子群优化算法求解弥散试验中的弥散系数和地下 水流速。 3 粒子群优化算法 粒子群优化算法Particle Swarm Optimization， 简称 PSO，又称鸟群觅食算法，由 Eberhart 和 kennedy[10-11]于 1995 年提出，属于模仿生物进化算 法的一种。目前已广泛应用于函数优化、神经网络 训练、 模糊系统控制以及其他遗传算法的应用领域； 在解决实际问题方面具有精度高、收敛快、容易实 现等优点[12]。 粒子群优化算法是一种模仿生物进化的智能寻 优算法。最初是受到鸟群集体觅食行为的启发，在 对鸟群捕食行为的研究基础上，利用各个飞鸟的觅 食方法和群体之间的相互协作，使得整个鸟群的行 为产生从无序到有序的演化过程，实现群体在问题 求解空间的寻优。 假定鸟群觅食情景为鸟群随机分散在固定 区域空间搜索食物，不知道食物的具体位置，但 可以感知自己与食物的距离，也可以感知到其他 个体位置及其他个体与食物的距离，鸟群找到食 物最简单有效的方法就是搜索当前情况下离食物 最近的鸟的位置，然后只搜索那只鸟的领域空间， 从而可以缩小搜索范围。在将上述鸟群觅食行为 抽象为数学模型时，认为每只鸟都遵循以下 3 条 理想化规则 ① 避免和空间其他同伴发生碰撞； ② 每只鸟飞行的速度都尽量和周围同伴保持 一致； ③ 飞行时都尽量朝着自己认为较好的位置靠近。 粒子群优化算法在求解问题时，问题的解对 应于区域中的食物，种群中的鸟均为没有质量和 体积的粒子。用空间位置、速度和适应度函数 3 项指标表示每个粒子特性。粒子的适应度函数值 通过目标函数计算得到，其大小决定粒子所处空 间位置的好坏；此时，粒子个体的速度决定了自 己的移动距离，速度随自身以及其他个体的移动 经验进行动态调整，粒子搜索食物的过程则模拟 为寻优过程，粒子在迭代过程中，通过两个极值 来更新自己的空间位置，一个是粒子自身所经历 位置中计算得到适应度最优的解，即个体极值 pbest，另一个是目前整个种群最优解，即全局 极值gbest。迭代过程不是完全随机的，每个粒 子根据式4和式5来更新自身速度和位置，每更 新一次位置，计算一次适应度函数值，通过比较 新粒子适应度值的 pbest 和 gbest，更新个体极值 和全局极值，直至找到最优解。 , 1 ,,,,1 12 2 i j tttttt i ji ji jji j vvc r pxc r gx   4 , 11 ,, i j ttt i ji j xxv   5 式中 , t i j v表示第 i 个粒子在第 t 次迭代过程中的第 j 维速度，其中 i[1,n]；n 为粒子种群数；j[1,d]， d 为粒子所处空间维数；ω 为惯性权重因子，反映 粒子当前状态对下一状态的影响，平衡全局和局部 搜索；c1和 c2为正常数，c1叫认识学习因子，主要 影响局部搜索，提高搜索精度，c2叫社会学习因子， 反映粒子间协同合作与知识共享的群体历史经验， 有利于全局搜索；r1、r2为[0,1]间服从均匀分布的随 机数； , t i j p为第 i 个粒子在第 j 维空间第 t 次迭代过 程中个体最优位置； t j g为所有粒子在第 t 次迭代过 程中全局最优位置的第 j 维坐标。 , t i j x表示第 i 个粒 子在 t 次迭代过程中第 j 维坐标。 在迭代过程中，为 保证算法的收敛性，一般要设置速度的上限 max v， 取待优化参数变化范围的 1020； 并设速度下限 minmax vv ，其中负号表示方向；当粒子速度 v超 出规定范围时， 即 max vv或 min vv， 则令 max vv或 者 min vv。 所有粒子的位置好坏均由适应度函数值来评 价，对于求解 fx最小值问题时，个体极值 pbest 根 据下式进行更新 11 1 if otherwise ttt iii t i t i xf xp p p        6 全局极值根据所有个体极值中最好的位置进行 更新。 4 水动力弥散系数的求解 4.1 基于粒子群算法求解水动力弥散系数 根据上述原理，可以将粒子群算法中鸟群觅食 的食物位置 x 概化成αL,v的函数，由鸟群个体的飞 行速度来更新飞鸟的位置即 αL和 v 的解，再根据 ChaoXing 第 6 期 梅杰等 基于仿生原理的粒子群算法求解江安校区水动力弥散系数 101 适应度函数判断解的优劣。本文以模型中相应时间 t 下的示踪剂浓度 c 与实际观测值的均方根差值 fitness 作为目标函数，来求解弥散系数，也作为粒 子群算法中个体的适应度函数，其数学表达式为  2 LL 1 , , / N ii i fitnessvccvN    min 7 式中 L , fitnessv为关于 αL和 v 的目标函数；ci为 实测浓度值； L , i cv为算法计算浓度值；N 为样本 个数，这里为浓度观测次数。适应度函数 fitness 值 越小，优化解越好。 利用粒子群优化算法求解本次弥散试验的弥散 系数时，粒子数量为 40；最大迭代次数为 10 000 次； 学习因子 c1、 c2均为 2； 最大速度 vmax0.1xmax– xmin；最小速度取其负值 vmin –vmax；惯性权重 ω 0.6。主要参数赋值见表 1。 表 1 主要参数赋值及范围 Table 1 Values and ranges of major parameters 主要参数 v/md–1 αL/m N ω c1 c2 赋值 0.000 40.04 110 40 0.6 2 2 待求参数纵向弥散度和地下水实际流速的取值 范围分别为 αL∈[1,10]和 v∈[0.000 4,0.04]。使用粒 子群算法计算得到 αL6.42 m， v0.024 m/d。 代入计 算得 DLαLv0.154 m2/d。 4.2 标准曲线法和最小二乘法求解 为检验粒子群算法计算结果的有效性和参数的 精度， 采用标准曲线法和最小二乘法计算弥散系数， 并进行结果对比。 根据 Sauty 通过计算得到式1问题的解为[13] 2 rrr 0.5 r r ,exp1 4 KP ct Pt tt     8 其中，   2 0.5 rmaxrmax rmax 0.5 21 rmax rmaxrc L1 exp1 4 1 / / / P Ktt t tPP cc ctt t rP         9 式中 P 为 Peclet 数；tc为纯对流时间；trmax是无量 纲时间的最大值；cmax是最大浓度；r1是观测井到 投料井的距离。以 P 为参数，无量纲时间 tr和无量 纲浓度 cr为横纵坐标，画出标准曲线。将测得的浓 度换算为无量纲量 cr，并与相应观测时间 t 值绘在 尺度与标准曲线图相同的半对数坐标轴上，然后将 该图置于标准曲线图上，使两图横坐标重合，直到 大部分实测点落在某一条标准曲线上时为配好，求 得参数 P1，拟合效果见图 3。最后根据式9求得 纵向弥散度为 4.6 m。 图 3 实测数据点与标准曲线拟合关系 Fig.3 Fitting relationship between measured data points and standard curve 最小二乘法是一种数学优化技术，其原理来 自于统计学中的极大似然估计原理；可保证各实 测点到回归线纵向距离的平方和最小，并使计算 出的回归方程最能代表实测数据所反映出的趋 势，最小二乘法是一种重要的曲线拟合方法。可 以将实测数据导入 Origin 软件中，根据式3设置 回归方程的基本形式，通过自带的最小二乘优化 回归，即可得到纵向弥散度为 6.34 m。3 种方法的 计算结果见表 2。 将 3 种方法所得参数和实际流速代入式3，计 算浓度随时间变化曲线，并与实际浓度变化曲线图 4对比分析，分析结果见表 2。 5 结 论 a. 采用粒子群算法求解弥散系数和地下水流 表 2 不同求参方法的计算结果比较 Table 2 Comparisons of parameters calculated results from different calculation s 计算方法 纵向弥散度 αL/m 实际流速 v/md–1 标准误差/ 相关系数 R2/ 粒子群算法 6.42 0.024 1.14 95.08 标准曲线法 4.60 0.019 1.63 90.07 最小二乘法 6.34 0.028 1.49 91.66 注相关系数 R2表示计算曲线与实测数据点的相关性，越接近于 1，拟合越好。 ChaoXing 102 煤田地质与勘探 第 47 卷 图 4 3 种计算方法曲线与实测数据拟合关系 Fig.4 Fitting of the three calculated curves with the measured data 速，应用于四川大学江安校区试验场，求得弥散系 数为 6.42 m，地下水流速为 0.024 m/d。 b. 对比最小二乘法、 标准曲线法和粒子群算法 的计算结果发现，粒子群算法计算结果的标准误差 最小，相关系数最大，表明该方法在天然流场下弥 散度的计算中有较强的适用性。 c. 粒子群算法可以同时求解多个参数，对水文 地质求参问题具有应用价值，但其自身存在稳定性 等问题，需进一步讨论研究。 参考文献 [1] 吴耀国，田春声，李云峰，等. 判定地下水二维水动力弥散参 数的直线图解法[J]. 煤田地质与勘探，1997，25233–35. WU Yaoguo，TIAN Chunsheng，LI Yunfeng，et al. The linear graphic for discriminating two-dimensional hydrody- namic dispersion parameter of underground water[J]. Coal Ge- ology Exploration，1997，25233–35. [2] 刘兆昌. 地下水系统的污染与控制[M]. 北京中国环境科学 出版社，1991. [3] 张富仓，康绍忠，潘英华. 饱和–非饱和土壤中吸附性溶质水 动力弥散实验研究[J]. 水利学报，2002，16384–90. ZHANG Fucang，KANG Shaozhong，PAN Yinghua. Experi- mental study on hydrodynamic dispersion of adsorption solute in saturated-unsaturated soil[J]. Journal of Water Resources，2002， 16384–90. [4] 张宇峰，张雪英，徐炎华，等. 土壤中水动力弥散系数的研究 进展[J]. 环境污染治理技术与设备，200378–12. ZHANG Yufeng， ZHANG Xueying， XU Yanhua， et al. Progress in study of hydrodynamic dispersion coefficient in soil[J]. Tech- niques and Equipment for Environmental Pollution Control， 200378–12. [5] 李世钰，刘国东，王亮，等. 基于蝙蝠算法的适线法求解水动 力弥散系数[J]. 煤田地质与勘探，2016，44481–85. LI Shiyu，LIU Guodong，WANG Liang，et al. Computation of hydrodynamic dispersion coefficient with curve fitting based on bat algorithm[J]. Coal Geology Exploration，2016， 44481–85. [6] ZOU S，PARR A. Estimation of dispersion parameters for two-dimensional plumes[J]. Ground Water，1993，313 389–392. [7] 曹剑锋，迟宝明，王文科，等. 专门水文地质学[M]. 北京 科学出版社，2006. [8] 陈崇希，李国敏. 地下水溶质运移理论及模型[M]. 武汉中 国地质大学出版社，1996. [9] 王永峰. 求解半无限长多孔介质柱体水动力弥散系数的一种 最小二乘法[J]. 地下水，2010，3236–8. WANG Yongfeng. The least square gaining the hydro- dynamic dispersion coefficient of the semi-infinite pore medium cylinder[J]. Groundwater，2010，3236–8. [10] EBERHART R C， KENNEDY J. A new optimizer using particles swarm theory[C]//Proceed Sixth International Symposium on Micro Machine and Human Science，Nagoya，Japan，1995 87–92. [11] KENNEDY J，EBERHART R C. Particle swarm optimiza- tion[C]//IEEE International Conference on Neural Network， Perth，Australia，1995168–173. [12] 李雅琼. 基于粒子群算法的遗传算法优化研究[J]. 兰州文理 学院学报自然科学版，2017，31155–60. LI Yaqiong. Research of genetic algorithm based on particle swarm optimization[J]. Journal of Lanzhou College of Arts and SciencesNatural Science Edition，2017，31155–60. [13] 刘岩磊，王庆来. 不同流场情况下地下水弥散试验方法分析[J]. 地下水，2014，36210–13. LIU Yanlei，WANG Qinglai. Analysis of groundwater dispersion test s for different flow fields[J]. Groundwater，2014， 36210–13. 责任编辑 周建军 ChaoXing