基于遗传算法的水源热泵系统抽灌量优化配置

王 麒，李松青，王心义，3，徐流洋，姬红英，夏大平

（1. 华北水利水电大学地球科学与工程学院，河南郑州 450000；

2. 河南理工大学资源环境学院，河南焦作 454000；3. 中原经济区煤层（页岩）气河南省协同创新中心，河南焦作 454000；4. 河南理工大学能源科学与工程学院，河南焦作 454000）

1 研究背景

水源热泵技术是一种利用清洁可再生能源的技术，即将地球水体所储藏的太阳能或地热能作为冷热源，采用热泵原理进行能量转换的供暖空调系统，其中可以利用的水体包括浅层地下水、河流、湖泊以及海洋［1-3］。若以浅层地下水为水源，为保证持续稳定的冷热能供应并避免环境水文地质负效应，必须考虑地下水抽水和回灌的问题，且应结合当地的水文地质条件及建筑物布局来配置抽水井和回灌井（简称抽灌井）及其运行模式［4-5］。

作为水源热泵系统设计中的关键环节之一，抽灌井的配置问题受到众多学者的关注。Wu 等［6］、Busby 等［7］认为场地水文地质条件是影响抽灌井设计的首要因素，并提出了抽灌井设计的水文地质要求。Russo等［8］利用神经网络方法对水源热泵系统运行中抽灌井周围温度进行分析预测，为抽灌井的合理设计提供理论参考。Andrews［9］建立地下水渗流场及温度场的耦合模型，模拟了水源热泵系统运行时对地下水温度场的影响，有效地指导了抽灌井井位的布置。Gunawardhana 等［10］研究了气候、水文地质和水源热泵运行等因素影响下温度的传播规律；结果表明，气候变化的影响主要集中在地下浅地表地带，地下水补给和水平流动阻滞着热量的向上运移，含水层注入速率、注入深度和各向异性是避免浅层地下临界热干扰的重要参数。Zhou等［11］和Gao等［12］利用水热运移耦合方程进行数值模拟，根据模拟结果提出了抽灌井布局的合理方式。刘立才［13］采用地下热传递数值法、能量均衡法和热量储存法评价了北京平原地区单位面积含水层利用地下水热泵系统可承载的能量负荷，为浅层地热能的开发和保护提供了技术支持。徐玉良等［14］建立地下水源热泵抽灌系统的三维模型，根据“热贯通”发生的时间和程度，对抽灌井的布局进行了优化选择。显然，为防止地下水回灌造成的“热贯通”，利用现代数学方法优化配置抽灌井布局的成果较多；而对于开采井和回灌井位置确定后，如何合理配置抽水量和回灌量（简称抽灌量），以使水源热泵系统水泵运行能耗最小方面的研究较少。

遗传算法作为一种随机搜索与优化算法，具有良好的适应性和全局优化能力，被广泛用于组合优化、模式识别、规划策略等领域。王占武［15］利用遗传算法进行多目标最佳路径的选取，赵大兴等［16］基于高适应度值遗传算法进行最优路径规划。袁成林［17］提出一种基于对应连通子图交叉混合遗传算法，对解决单目标旅行商问题有很好的效果。徐旭等［18］以利用河套灌区曙光实验站的土壤剖面分层含水率和溶液浓度、表土含水率、实际腾发量等观测资料，基于遗传算法开展了田间尺度土壤水力参数和溶质运移参数优化研究。总之，遗传算法在求取最优解时具有明显的优势。

针对抽灌量的配置不同其水源热泵系统运行能耗也不同的问题，本文基于水文地质条件和抽灌井布局及需水量情况，建立了水源热泵系统水泵运行能耗最小的控制模型，并利用遗传算法对模型进行了求解，确定了抽灌系统的抽灌量，为水源热泵系统的合理运行提供了技术支撑，同时也为水文地质条件类似地区水源热泵系统运行模式的科学调配提供了参考。

2 基本模型

水源热泵系统中抽灌量的配置涉及地下水运动的基本理论和水泵能耗原理。

2.1 地下水运动数学模型针对地下水向承压完整井的平面二维非稳定流运动，对应的数学模型为［19］：

式中：T为含水层导水系数，m2/d；Q为单井抽水量，m3/h； μ*为弹性释水系数；r为计算点至抽水井的距离，m；t为抽水时间，min；S为抽水所引起的r处水位降深，m。

式（1）已由Theis C.V.［19］给出其解析解，即t时刻距离抽水井r处的水位降深S为：

式中：S为水位降深，m；W ( u )为井函数；u为中间函数。

回灌可以看作是抽水的逆过程，只是水的运动方向不同而已，因此式（2）同样适用于回灌过程水位升值的计算［20］。在抽水与回灌同时发生的情况下，根据叠加原理，某点的水位降深等于抽水产生的水位降深与回灌产生的水位升值之差，计算式为：

式中：Sw为某点的水位降深，m；S抽为抽水引起的水位降深，m；S回为回灌引起的水位抬升值，m。

2.2 优化配置模型以单井回灌量为变量，以水泵运行能耗最小为目标函数，可建立水源热泵系统抽灌量的优化配置模型。

抽灌过程中水泵的能耗与抽灌量及地下水埋深紧密相关；水运移过程中的水头损失除受抽灌量影响外，还与输水管道布局及直径有关。因此，目标函数中由三方面构成：（1）抽水过程中克服地下水位至地面位置水头的能耗；（2）水运移过程中产生的能耗，包括沿程水头损失和局部水头损失；（3）回灌过程中克服含水层水头差的能耗。

目标函数：

约束条件包括：（1）单井回灌量须小于单井回灌能力；（2）单井水位降深在设计降深范围内；（3）回灌总量与抽水量相等，以实现对地下水资源的可持续开发利用。即：

约束条件Ⅰ：

约束条件Ⅱ：

约束条件Ⅲ：

式中：Q抽为取水井的抽水量，m3/h；Q回为单井的回灌量，m3/h；H1为初始水位埋深，m；H0为初始水位，m；Srr为抽水井对自身产生的降深，m；Sij为回灌井的降深，i=j 时为回灌对自身产生的降深，i≠j时为第i回灌井对第j回灌井产生的降深，m；η 为水泵的工作效率；γ 为水的容重，kg/m3；F 为水泵总能耗，kW·h；F抽为抽水能耗，kW·h；F回为回灌能耗，kW·h；F损为水运移过程中造成的水泵能耗，kW·h。

3 应用模型建立

3.1 抽灌井布局研究区位于开封市东郊，地貌上属黄河冲积扇区。场地基本平整，高程约为70 m。取水目的层为埋深70 ～210 m承压含水层，平均厚度为76 m，天然条件下水位埋深为16.7 m。含水层属第四系松散类孔隙含水层，岩性为粉土、中细砂、中砂、泥质细砂、细砂、砂质黏土和黏土的互层，富水性较强，适合水源热泵系统的开发应用。

研究区内共布设5眼管井（图1），井间距离最小为21.10 m，最大为57.24 m，分布面积约为1193 m2。由于5 眼井间的距离较近、目的含水层均属松散岩层且无断层的存在，天然条件下5 眼井的水位基本相同，相互之间水力联系密切，地下水位之间干扰明显。

管井直径280 mm，设计井深207 ～211 m之间，取水段75 ～200 m，开采井稍深于回灌井，管井结构见图2。井3和井5与水源热泵机组连接管道的直径为150 mm，回灌井连接管道的直径为100 mm。

水源热泵系统中设计抽灌井有两种工作模式：（1）以井3 作为抽水井，以井1、井2、井4、井5为回灌井；（2）以井5作为抽水井，以井1、井2、井3、井4作为回灌井。

图1 抽灌井位及管道布局

图2 管井结构图

3.2 应用模型构建钻探资料表明，1193 m2的有限研究区内，埋深70 ～210 m的取水段为承压水，且含水层水平、岩性一致，为均质各向同性、无垂向补给的承压含水层。由于抽灌井为井底直接坐落于隔水层之上的承压完整井，因此，地下水的运动可用式（1）来表征，水位降深可用式（2）和式（5）来计算。

针对水源热泵系统抽灌井管道布局（图1），优化配置抽灌量的数学模型由式（6）—式（9）构成。其中以井3为抽水井的抽灌模式下，水流在管道运输中水头损失F损的计算方程［21］为：

式中：λ为沿程阻力系数；ε为局部阻力系数；L0为抽水井3至水源热泵机组的管道长度，m；L1为水源热泵机组到A点（图1）的管道长度，m；L2为A点到D点的管长，m；L3为D点到井1的管长，m；L4为D点到井5的管长，m；L5为A点到C点的管长，m；L6为C点到井2的管长，m；L7为C点到井4的管长，m；Q1、Q2、Q3、Q4、Q5分别对应各井的流量，m3/h；d1、d2为管道的直径，m。

同理，以井5为抽水井的抽灌模式下，F损的计算方程为：

式中：L0为抽水井5到水源热泵机组的管道长度，m；L1为水源热泵机组到A点的管道长度，m；L2为A 点至回灌井1 的管道长度，m；L3为A 点至B 点的管道长度，m；L4为B 点至回灌井3 的管道长度，m；L5为B点至C点的管道长度，m；L6为C点至回灌井2的管道长度，m；L7为C点至回灌井4的管道长度，m。由此可见，水泵的能耗与抽灌量的配置紧密相关。

4 参数确定

4.1 水文地质参数为确定模型计算所需要的水文地质参数，在水源热泵供暖前后共开展了两次抽水试验。第1 次抽水试验于2016 年11 月27 日进行，以井3 为抽水井，其余4 眼井为观测井，抽水流量为90 m3/h，试验进行1560 min 后，水位降深6.5 m，基本达到稳定状态。

表1 井4观测时间和降深

井4 的初始水位为16.9 m，距离抽水井距离为21.72 m，水位降深数据如表1。利用井4 观测数据，基于AquiferTest 软件中的Theis 配线法可求得水文地质参数：导水系数T为1050 m2/d，渗透系数K为13.8 m/d，弹性释水系数μ*为2.05×10-4。拟合结果如图3，计算公式见式（2）—式（4），计算参数见表2。

表2 各观测孔数据求得的水文地质参数

图3 Theis配线法拟合曲线

表3 第2次抽水试验求得的水文地质参数

同理，利用其它各观测井水位监测数据可求得水文地质参数，以各参数的平均值作为第1次抽水试验确定的水文地质参数，见表2。

为验证第1次抽水试验计算的水文地质参数准确性，并探究水源热泵系统对水文地质参数的影响，在一个供暖期后，于2017年5月9日开展了第2次抽水试验。此次抽水试验以井5为抽水井，其余井为观测井，抽水流量为110 m3/h，试验进行1380 min后，各井水位基本达到稳定状态。与第1次抽水试验计算类似，求得的水文地质参数统计结果如表3。

对比两次的计算结果可以发现，第2次抽水试验计算的水文地质参数大于第1次的计算结果，原因是经过一个取暖期的抽水，承压水位降低，含水层透水性变好，导致导水系数和渗透系数有所增大；同时承压水位的降低，导致上部含水层部分被疏干，进而使得弹性释水系数减少。由于两次抽水试验的计算结果差别不大，因此选用其平均值作为研究区的水文地质参数，即导水系数T为1052.5 m2/d，渗透系数K为13.85 m/d，弹性释水系数μ*为7.52×10-4。

4.2 抽水量和单井最大回灌量抽水量和进水温度是参考开封市其他供暖小区水源热泵运行情况确定的。设计抽水量基本稳定在90 m3/h，进水温度维持在19.6 ℃，回水温度约14 ℃，可满足研究区供暖负荷的需求。

为了确定单井的最大回灌量，共进行了3眼井的回灌试验。每次回灌试验均是在抽水水位达到基本稳定后才开始实施的，直至抽水井和回灌井的水位均基本稳定后再停止试验。

回灌试验于2016年11月30日实施，单井回灌时间均在36 h以上。以井3为抽水井，井1、井2、井5分别为回灌井。试验结果表明，井3的抽水量为90 m3/h，井1、井2、井5的稳定回灌量及水位抬升值见表4。

表4 回灌试验分析数据

根据稳定回灌量与水位抬升数据，利用最小二乘法可建立线性回归方程：

式中：y为稳定回灌量，m3/h；x表示此状态下的稳定水位抬升值，m；相关系数R2为0.9998。

《地源热泵系统工程技术规范》（GB50366-2005）要求回灌井的水位抬升值不应大于5 m，因此估算出单井最大回灌量Q回max为42.5 m3/h。

4.3 水头损失系数水头损失系数包括沿程阻力系数和局部阻力系数。对于钢管及铸铁管管道，沿程阻力系数λ 的计算方式为：

当水流速度v≥1.2 m/s时，

式中d为管道内径，m。

当水流速度v＜1.2 m/s时，

水流在转弯时会产生局部水头损失。根据场地抽灌井管道布局，参见《给水排水手册》，确定局部阻力系数ε 为0.6。

4.4 其它参数模型中计算所需的参数还包括水泵效率、供暖期时间。依据购置的水泵型号，水泵效率η 取0.8。开封市供暖时间起于11月15日，止于次年3月15日，供暖时间为120 d。

5 模型求解与结果分析

5.1 遗传算法针对上述建立的水泵运行能耗最小的目标函数及其约束条件，可利用遗传算法求解以实现抽灌量的优化配置。遗传算法计算路线见图4。

遗传算法是基于生物进化原理的一种具有鲁棒性的随机搜索算法，具有良好的适应性和全局优化能力，在优化模型求解中应用广泛［22-23］。其具体的计算过程是：①参数设置。主要包括遗传算法种群数目、遗传代数、交叉率、变异率。②初始化种群。在全局空间中随机生成一组解构建染色体，作为初始种群。③适应度函数确定。为了能够选出优良的子代，使得搜索朝着最优解的方向进行，需要确定合适的适应度函数。以目标函数方程为基础，根据求解目的进行适当变换，作为遗传算法求解的适应度函数。④选择、交叉、变异。根据适应度值，采用轮盘赌方法选择个体重组为新的群体，适应度值越大，被选择的几率就越大；通过交叉方式，将两个个体重组产生新的个体，生成［0，1］之间的随机数与交叉率比较，小于交叉率执行交叉操作；通过变异算子，改变个体中某个字符，使种群进行变异，生成［0，1］之间的随机数与变异率进行比较，小于变异率执行变异操作。⑤当达到最大遗传代数或适应度值趋于稳定时，算法终止，选取适应度最好的结果为算法运行结果。

图4 遗传算法计算路线

在遗传算法中，种群数目、遗传代数、交叉率和变异率是影响遗传算法性能的重要参数。种群数目越多，获得优良种群的概率越大，但是运行速度就越慢；遗传代数过小则无法完成优化搜索任务，过大则导致计算时间过长；交叉率越大，种群丰富程度越高，但优良种群被破坏的概率也越大；变异率越大易跳出局部极值找到最优解，但是过大易沦为随机搜索算法［24］。参考前人成果［25-27］，本文设置的参数为：种群数目为100，最大遗传代数为100，交叉率为0.8，变异率为0.1。

基于水泵能耗最小的目的，构建的适应度函数f ( x )为：

式中：F为目标函数方程。用F2是为了增大个体之间的差异，以便于选择更优的个体，加速最优解方向的搜索。

5.2 模型求解针对以目标函数为基础建立的适应度函数，应用遗传算法来求解。模型计算所利用的主要参数见表5，表5中参数M为含水层厚度，m。

表5 模型主要参数设定

图5 井3为抽水井时遗传算法运行过程

图6 井5为抽水井时遗传算法运行过程

当以井3为抽水井和以井1、井2、井4、井5为回灌井时，遗传算法的迭代运行结果如图5。图中黑色点代表每代种群中目标函数最小值，红色点代表每代种群中平均目标函数值。显然，遗传代数超过8次时，目标函数值（适应度值）即趋于稳定。计算结果为：Q1为24.9 m3/h，Q2为15.5 m3/h，Q4为12.3 m3/h，Q5为37.3 m3/h时，水泵在供暖期内的能耗最小，其值为184 588 kW·h。

当以井5 为抽水井和以井1、井2、井3、井4 为回灌井时，遗传算法的运行结果如图6。显然，遗传代数超过13次时，目标函数值（适应度值）即趋于稳定。计算结果为：Q1为42.5 m3/h，Q2为33.4 m3/h，Q3为0 m3/h，Q4为14.1 m3/h时，水泵能耗最小，其值为183 020 kW·h。

5.3 结果与分析在以井3为抽水井的抽灌模式下，经历一个供暖期后，与初始水位相比，回灌井1的水位抬升值为+0.364 m，井2抬升值为+0.280 m，井3水位降深为1.743 m，井4水位降深为0.007 m，井5 水位抬升值为+0.769 m，供暖期内的水泵能耗为184588 kW·h。在以井5 为抽水井的抽灌模式下，计算得出回灌井1的水位抬升值为+0.687 m，井2抬升值为+0.632 m，井3抬升值为+0.110 m，井4抬升值为+0.383 m，井5水位降深为1.774 m，供暖期内的水泵能耗为183 020 kW·h。结果如表6。

表6 两种抽灌模式下计算结果

比较两种抽灌模式下的计算结果可以看出，抽水井水位降深相差不大，各回灌井水位降深均在允许设计降深5 m范围内；在水泵能耗方面，一个供暖期内，以井5为抽水井的抽灌模式比以井3为抽水井的抽灌模式节省了1568 kW·h，减少了水源热泵运行的成本。此外，根据以井5为抽水井的抽灌模式下抽水和回灌量的最优配置结果来看，井3的回灌量为0，仅由井1、井2、井4进行回灌即可满足最优抽水和回灌量相等的配置要求，也即井3可作为其他回灌井轮替检修时以及为避免堵塞而进行“回扬”时的备用回灌井。

本文研究表明：（1）水源热泵系统建立前，应根据水文地质条件并结合建筑物布局来优化抽灌井位置和管道布局，以避免抽灌井工程的浪费，并保证抽灌井和输水管道布局的合理性；（2）因建筑物位置限制而无法按优化方案布局抽灌井和管道时，应从水源热泵运行能耗最小角度出发来优化配置抽灌量；（3）水源热泵系统运行过程中，应严格按照优化方案进行抽灌量的科学调度，克服过去靠经验而行的缺点；（4）为保证抽水量全部回灌至地下，避免因故障或堵塞而致回灌水溢出地面，必须定期对回灌井进行检修和“回扬”，因此在满足优化配置回灌井数的前提下应增设1 眼回灌井以备用；（5）由于现场调节回灌井阀门开启度（使其处于优化回灌量）较为费时费力，而回灌井检修或“回扬”时间一般较短，故某回灌井检修或“回扬”期间，可暂将其回灌量转移至备用井，其他回灌井的水量可不调整，以便检修或“回扬”结束后迅速恢复优化抽灌状态。

6 结论

（1）基于水文地质条件和抽灌井位置及管道布局，建立了以水泵运行能耗最小为目标函数，以回灌量小于单井回灌能力、水位降深在设计允许降深范围内、回灌总量与抽水量相等为约束条件的数学模型，该模型是符合实际的，为水源热泵系统抽灌量优化配置奠定了基础。

（2）针对水源热泵系统抽灌量优化配置数学模型，在确定种群数目、遗传代数、交叉率和变异率等参数的基础上，构建了适应度函数，应用遗传算法进行了求解，为多约束条件下优化模型的求解提供了可借鉴的方法。

（3）优化计算结果显示，以井5为抽水井，井1、井2、井4为回灌井，井3为回灌备用井时的水源热泵系统运行能耗最低，属最优配置；水源热泵系统可按井5 抽水量为90 m3/h，井1、井2、井4回灌量分别为42.5 m3/h、33.4 m3/h、14.1 m3/h 的模式进行调度，克服了过去靠经验调配回灌量的缺点，为科学快速调配抽灌模式提供了依据。