地下水机器学习方法研究
——水位监测数据驱动的区域补排边界识别

齐永强 ，李文鹏 ，郑跃军 ，李 慧 ，王成见

（1.北京北水国际科技有限公司，北京 102206；2.中国地质环境监测院，北京 100081；3.青岛市水文局，山东 青岛 266071）

近年来我国地下水监测工作取得了长足的进步，获得了高频率高密度的水位监测数据，是极佳的地下水系统信息来源。但目前监测数据未得到充分利用[1-2]。在传统的地下水模型应用中，水位监测数据处于从属地位，与区域概化所得的水文地质边界以及水文地质参数共同决定模型的最终表现。随着水位监测的密度、频率和精度不断提升，监测数据中蕴含的信息可以更多更好地指导水文地质边界的识别和水文地质参数的确定[3]。本文应用了一种机器学习方法，用于提取区域地下水补排边界。

地下水工作者常借助地下水流模型系统整合已有的含水层信息，认识含水层行为，测试溶质运移假说[4-5]。地下水流模型模拟是使用解析方法或数值方法对饱和带水分运移控制方程求解的过程。水流模型的计算结果是在给定的初始条件和边界条件下，水头在空间各处和时间各点上的分布[6]。一般使用数值方法对控制方程求解，因为解析方法求解时对初始条件和边界条件有严格的限制，而客观世界中的情况极为复杂，只有数值方法能够处理与之相对应复杂性的模拟问题[7]。由于地下水数据存在长期而普遍的稀缺性，模型工作者需要对模型结构和参数进行插值和调整。对这一过程的校准通常通过试错法进行，即反复调整模型中指定的输入参数，直到水头的模拟值与监测值达到一致。在数据和参数量较大的情况下，也可以使用专门为参数识别而设计的计算程序实现自动反演，如李家兴[8]利用GAN实现了水头场与底层参数场间的双向预测，李竞生等[9]利用遗传算法尝试识别含水层参数。

在传统地下水模型研究中，地下水水位监测值常位于模型构建过程的下游，其作用主要体现在模型校正过程中对模型参数的调整。这主要是由于之前地下水水位监测数据无论从时间密度或空间密度均存在严重短缺[5]，其中蕴含的信息不足以勾画含水层的全貌。相应地，研究人员在实际利用地下水水位监测数据之前，一般需要针对系统补给和排泄边界和含水层参数分布作出一系列假设，即构建概念模型。但是假设中不符合实际情况的错误信息，将在模型中不断传导和放大，造成模型的多解性[10]。单纯使用地下水水位的实测数据修正和微调地下水模型不仅效率较低，而且容易陷入循环论证的困局。本文旨在于建立一套可行的技术方法，能够直接使用高密度、高频率的地下水水位监测数据对区域内地下水补给和排泄边界进行识别，帮助构建和修正含水层概念模型。

机器学习是人工智能的一个分支，其核心目的是从大量数据中自动分析获得规律[11]，近30年已发展成为一门多领域交叉学科，涉及概率论、统计学[12]、计算机图形学、逼近论、凸分析[13]、计算复杂性理论等多门学科。机器学习通常基于神经网络架构，以认知科学、神经科学等为依托，设计类脑智能的学习模型[14-15]。机器学习常分为监督学习和无监督学习。从高密水位数据中提取地下水补给和排泄边界属于模式识别，是一种形式的无监督学习。无监督学习指从现有的训练集出发，经过训练得出一定结果。与监督学习不同，无监督学习不需要人为标注，仅着重于捕捉和描述数据中存在的模式和规律[16]。机器学习与基于物理逻辑的传统模型不同，是以统计与算法为基础，将逻辑搭建交由神经网络自行生成，旨在挖掘数据中无法人为识别的信息与信号[17]。

本文介绍了一种地下水系统补排边界的识别方法，其核心是在三维空间中构建研究区域，针对区域上的任意点判断其落入补给区或排泄区的概率，引入机器学习处理优化问题中的梯度下降[18-19]。梯度下降法由法国数学家Cauchy[20]发明，由Robbins等[21]重新提出，后由Rumelhart等[22]发扬并运用于机器学习领域。近十几年来，机器学习受到越来越多的学者与业界关注[23-24]，梯度下降算法也得到了广泛的研究与长足的发展[25-26]，被应用于约束优化[27]、正则优化[28-29]和函数复合[30]等不同的场景中。机器学习实践中，常用损失函数描述模拟结果和实测数据的差异，在指定数据集上，损失函数越小，模型预测越准确。为了确定何种参数组合下的损失函数最小，可通过参数在损失函数的“场”中向损失函数值减小的方向移动。其中损失函数值下降最快的方向称为负梯度方向，使用的算法称为梯度下降法，是机器学习中最常见的优化算法。

梯度下降法与地下水的流动具有天然的相似性。在地下水流场中的任意质点，如果允许其沿地下水流动，在下一时刻这一质点位于流场排泄区的几率增大；反之如果质点向地下水上游流动，下一时刻位于流场补给区的几率增大。这一理念在地下水动力学领域已有应用，最为典型的是美国地质调查局（USGS）开发的MODPATH质点追踪程序包[31]。然而MODPATH必须应用在已经求解的MODFLOW流场中，对于仅有水位监测信息的情形并不适用。本文描述的方法将尝试解决此问题。

1 研究方法

1.1 地下水监测网空间剖分

地下水水位监测网络的核心是三维空间中的一系列散点类水头值，位于每口监测井筛管的中心位置。若要通过内插法得到模型范围内任意点的水头值，需首先为散点建立空间网格，描述它们之间的拓扑关系。由于在地下水区域监测网络中，水平方向的长度（X、Y方向）一般远大于垂直方向（Z方向），将三维散点沿Z轴投影到大地平面上，使用德劳内三角化规则为二维平面上的点集建立三角网格。德劳内三角化是计算几何中常见的三角网格建立规程，其核心为不应有任何一个顶点在三角形网的任一三角形外接圆内部，从而避免产生狭长三角形（图1）。

图1 地下水监测网空间剖分Fig.1 Spatial discretization of the groundwater monitoring domain

1.2 求取水力梯度

h∈R|V|。当集合V中的h值均已知时，可以在集合

设三角网格M上的点集合为V，三角集合为F，函数h代表三角网格节点（监测井位）上的水位值，即F上求取函数h的梯度。定义一个针对三角集合F的水力梯度变换算子gradF，用来将集合V上的水位值映射为集合F上的水力梯度值：

式中：gradF——水力梯度变换算子。

对于F集合中的一个三角形j，gradF变换为：

式中：A(j)——三角形j面积；

hi——三角形j的第i个顶点上的水位值；

T——90°旋转矩阵；

——向量沿三角形平面旋转90°后的向量，

综上，可以将gradF变换的矩阵形式写为：

假设网格M中以某点N为顶点的三角形有k个，顶点N处的水力梯度为此k个三角形中水力梯度最大值。定义一个针对点集合V的水力梯度变换算子gradV，用来为顶点V集合上的任一顶点N生成唯一的水力梯度值：

式中：gradV——针对网格顶点N的水力梯度变换；

ji——经过顶点N的所有三角形。

1.3 补给区和排泄区的定义与识别

地下水补给区和排泄区做如下定义：在同一含水层内部指定随机分布的一系列质点，质点沿水位场梯度向水位较低的方向运动，在给定时间后，质点密度较高的区域称为排泄区；反之，质点沿梯度向水位较高的方向运动指定时间后，质点密度较高的区域称为补给区。

在三角网格M内定义任意质点，其空间位置属性p随时间t的变化可以表示为：

式中：p——质点的空间坐标；

t——质点时间。

此处时间为计算质点轨迹过程中使用的时间项，应与地下水监测网的实际监测时间进行区别。点位置为p0。对于给定质点时间步长 Δt＞0，Runge-

在计算质点轨迹时，本方法假定流场为稳定态，即在计算期间不发生变化。在初始质点时间t0，质Kutta[32-33]给出了质点在未来时间位置的近似值：

如果p(tn+1) 是质点在下一质点时间步的位置，pn+1即为p(tn+1)的四阶近似值。pn+1完全由当前质点位置pn和其位置函数f共同决定。

地下水中沿流场运行的质点，其位置函数由达西定律确定：

式中：v——质点流速；

q——达西流速；

ne——介质有效孔隙度；

K——含水层渗透系数；

I——本含水层的水力梯度。

在饱和渗流介质理论体系中，K值和ne值均为不随时间变化的固定值。在含水层的任一瞬时，一个质点的运动轨迹仅与其轨迹上的I值有关。假设一个流场的I值在给定时刻为已知，可以根据以上推导求取质点迹线，只要使用的 Δt值足够小，这条迹线的形态即为固定值，与K值和ne值无关。若要求取地下水流场中任意质点的迹线形态，理论上仅需流场的水力梯度（I）的分布函数，不需要含水层的渗透系数和孔隙度参数。

1.4 应用实现

对于同一含水层中任意给定的一组观测井，可以使用观测井筛管的中点生成一组点集合V，在XY平面上为点集V构造符合德劳内规则的三角网格，将网格重新投影至三维空间中即形成三角网格曲面M，M由点集V和三角集合F组成。在给定的时间切片上，点集V上的水头h值均为已知。在曲面M上任意点的水力梯度I，均可以通过gradF（三角形内）或gradV（三角形顶点处）计算得到。根据水力梯度在M上的分布函数，可以计算曲面上任意质点的近似运动迹线，监测井的空间分布越密集，质点运动的质点时间步长越短，此近似运动迹线就越接近真实迹线。并将计算逻辑在环境地学计算平台EnviFusion-CGS中实现。

EnviFusion-CGS是中国地质环境监测院和北水国际联合开发的环境地学可视化计算平台，用于环境地学领域基于数据的融合建模和知识发现。可针对多种环境要素（土壤、地下水、地表水、大气等）的不同维度数据（现场分析、仪器测试、数值模拟、地球物理等）进行融合分析。EnviFusion-CGS将本学科多源异构数据统一到四维时空的多组分属性数据体中，帮助研究人员以所见即所得的形式对海量数据进行可视化分析和统计挖掘。

2 示范区概况

大沽河是山东半岛最大的河流。向南流经烟台市的3个区县和青岛市的5个区县，长约179 km，最后注入山东半岛南部沿海岸的胶州湾。大沽河含水层位于大沽河流域的中下游，距离青岛市城区约50 km，面积约420 km2。大沽河含水层由大沽河流域冲积的松散沉积地层形成。在大沽河流域，农业灌溉是用水量最大的项目，同时也是区域经济可持续性发展的重要因素[34]。农业和轻工业对于大沽河含水层的经济发展尤为关键。

2.1 水文地质条件

大沽河流域位于山东半岛东部地区，地势较低，地表向南倾斜延伸至海岸线，古岘镇是该流域的分水岭。古岘镇北部是由花岗岩和变质岩构成的低矮山丘，海拔高度为50～200 m；大泽山海拔高度达736 m，为该地区最高的山。古岘镇的南部是胶莱盆地平原地区，该地区的平均海拔高度为20～50 m。大沽河附近的平原曾是一片广阔的山谷，山谷宽度达6 km，后经古代冲积物填充。山谷的底部由白垩系砂岩和页岩基岩组成，为承压含水层底板。山谷东西边界被冲积物覆盖，目前很难鉴定地貌。山谷中充满了砂质冲积物，是大沽河含水层的主要组成部分[35-36]。

潜水含水层由2个含水层组成：上部含水层（由黏质砂土和砂质黏土组成，渗透系数较低）和下部含水层（由渗透系数和孔隙度比较大的砂卵砾石组成）。地下水的总体流向为由东北向西南。大沽河含水层的上中下游分布有一系列的横断面，见图2。包气带厚度大小不一，为1～6 m。包气带由粉质黏土，黏质粉土以及砂土组成。南村水文监测站的北部地区（该监测站位于大沽河含水层的中游）包气带为较薄的砂层，是降雨入渗的天窗，入渗面积达9.2 km2。大沽河含水层渗透系数为2.3～159.3 m/d，给水度为0.05～0.18[37-38]。

图2 示范区水文地质条件Fig.2 Site map and hydrogeological cross sections

2.2 地下水位波动

青岛市大沽河中下游地下水系统监测始于1995年。2008—2010年青岛水文局建立了地下水水位自动监测网，在大沽河含水层范围内共建设147口自动监测井，对于研究农业灌溉抽水与含水层储变量关系至关重要。区域内春旱一般和春季灌溉同时发生，结果导致地下水水位严重下降。在灌溉的峰值季节，地下水水位呈现下降趋势，当8—9月降水丰富时，水位逐渐回升。入秋之后，水位逐步下降。12月至翌年2月农业灌溉较少，地下水水位相对稳定。全年地下水水位波动一般在3 m之内。

地下水资源的长期动态取决于气候变化周期和人类活动。根据1975—1995年地下水位连续观测结果，地下水资源经历了以下几个阶段：

（1）小幅下降期（1975—1980年）

该时期是地下水开采的初期，地下水水位开始逐年下降，降幅约为0.1 m/a。

（2）大幅下降期（1980—1984年）

该时期是地下水开采的高峰期，地下水水位大幅下降，降幅约为0.5～1.0 m/a。

（3）迅速补给期（1985年）

第九号台风登陆，青岛市连续4 d暴雨，地下水水位迅速上升，平均升幅达到了2.05 m。

（4）大幅下降期（1986—1989年）

由于持续的干旱天气，这一时期的地下水水位又开始急剧下降，降幅约为0.5 m/a。

（5）补给和稳定期（1990—1995年）

1990年青岛市降水量达到了900 mm，超出该地区多年平均降水量的50%。地下水也在这一年得到了迅速补给。随着1989年引黄济青工程正式通水，青岛市对大沽河含水层的依赖明显减小。

2.3 地下水补给

大沽河含水层的补给来自降水、农业灌溉、河流入渗以及少量的侧向补给。7—9月的降水量约占全年降水量70%，因此这一时期也是地下水补给的高峰期。贾立华[36]研究了大沽河地下水污染脆弱性计算得到了大沽河补给地下水的数值。利用DRASTIC绘制的大沽河含水层显示约有10%的含水层在这一时期的补给量达到约200 mm，部分地区高达400 mm。该区有一个高入渗天窗，面积为9.2 km2，也是地下水补给的重要来源之一[36-38]。

2.4 地下水排泄

大沽河含水层是地下水埋深较浅的潜水含水层。其砂质含水层的特性使得地下水和地表水相互作用密切。在20世纪80年代，地下水被大规模开采之前，由于淤泥层和黏土层的存在，局部地区可能是半承压含水层。随着季节交替，河流补给地下水或地下水排泄至河流。地下水开采时，含水层变成典型的潜水含水层。随着地下水水位下降，河流下渗补给地下水成为地下水和地表水之间交换的主导作用。从此，大沽河含水层以垂向循环为主，通过降雨或者灌溉补给；通过抽水和蒸散排泄。当地地下水动态存在季节性和年际变化[39-40]。

3 模型构建与讨论

大沽河含水层中地下水资源的补给和排泄受到人类经济生产活动的深刻影响。区域内大范围的农业种植区域依赖地下水进行灌溉，而且大沽河沿线存在大量水利工程构筑物。到2003年，沿河道建有8个中型水库，90个小型水库和1 223座闸门，对地下水的循环路径与模式产生了巨大影响。2013年1月1日——12月31日为研究期，在大沽河含水层范围内的147口监测井中，选取了78口具有连续逐日水位监测数据的监测井作为研究对象（其余监测井均存在不同程度的数据缺失）。仅使用监测井的空间位置信息和连续监测水位数据，在不作任何概念模型假设的前提下，使用机器学习算法对含水层的补给和排泄区域进行识别。

3.1 网格构建

根据78口监测井的X、Y坐标，基于德劳内规则，构建大沽河含水层的三角网格，如图3所示。生成的三角网格范围，见图3（b），与传统水文地质勘察工作所确定的“含水层边界”，见图3（a），略有不同。这是因为本算法仅具有内插功能，对于没有落在三角网格范围内的区域无法进行预测。当三角形网格过于粗糙时，会导致质点运动迹线不够平滑难以收敛，对图3（b）中的三角网格进行加密，见图3（c）。三角网格的加密仅为增加空间分辨率，不会增删原始网格中所携带的信息，其算法较为普遍，在此不再赘述。

图3 示范区地下水监测网的空间剖分Fig.3 Spatial discretization of the site

3.2 水力梯度场计算

在图3（c）所示的三角网格中，根据网格节点处的水位数值，计算网格三角形内部和节点上的水力梯度值。由于输入数据为逐日监测数据，所以分别计算365 d的水力梯度场，见图3（c）。含水层北部靠近山区，受地形坡度影响，水力梯度普遍较南部大。此外，靠近大沽河河床的位置由于含水层透水性较好，水力梯度普遍较小。

3.3 迹点与迹线计算

在图3（c）所示的三角网格中，首先为一定数量的质点随机给定初始位置，按照追踪手段依次绘制质点在下一个时间步长的位置并连成迹线，见图4（a）。初始质点的数量并没有具体的规定，但在计算能力允许的前提下尽量提高质点密度有助于更好地利用机器学习算法的优势。迹线的分布定性展示了流场的排泄方向。为了排泄区域定量化，在XY二维空间中划分了均匀方形网格，见图4（b）。依次计算方形网格中每一个网格周边迹点数量，并以此计算迹点的相对密度。按照排泄迹点密度大小得到迹点密度图，见图4（b）。

图4 地下水排泄迹线和排泄区识别（2013-01-01）Fig.4 Discharge path lines and discharge regions of the study site (2013-01-01)

在普遍意义上，此密度分布并不能直接等同于含水层的排泄区，因为水力梯度I的分布函数仅反映潜在的渗流方向，而实际流速受渗透系数K和介质有效孔隙度ne的共同约束。当含水层内存在强烈影响渗透系数或有效孔隙度分布的地质异常（如阻水构造）时，排泄迹点的分布不能用来推断排泄区的分布。反之，在大沽河含水层此类含水介质分布较为稳定，渗流特征较为均一的地区，此密度分布图可以看作含水层排泄区域的近似表达。可以看出含水层的地下水排泄基本以大沽河河道为主轴。偶有偏移的原因可能是：（1）傍河开采井的汇流作用；（2）古河道的地下渗流作用；（3）监测井密度不足造成排泄区形态失真。

类似地，将水力梯度方向反转，可以求取地下水的补给迹线和迹点密度分布图，见图5。含水层内的地下水补给作用明显集中在北侧的上游地区，含水层西北侧古岘镇——高岚村一线存在明确的分水岭，此线以西的地下水向西流动，以东的地下水向东流动。含水层东部边界存在明显的地下水流入，最终被大沽河吸纳。此外，含水层的南部底端存在局部的补给区域，与棘洪滩水库位置重合。

图5 地下水补给迹线和补给区识别（2013-01-01）Fig.5 Recharge path lines and recharge regions of the study site (2013-01-01)

3.4 地下水补排规律的时间变化

为考察本含水层补给区和排泄区随时间的变化规律，使用EnviFusion-CGS中内置的自动化参数设置，计算2013年每月1日计算补给区和排泄区的分布图，见图6。2013年含水层补给区与排泄区分布总体变化不大，其显著特征（西北部分水岭、大沽河道排泄带、南端水库补给等）均维持稳定状态。含水层与东侧边界的沟通方式在雨季来临后（8月）发生一定变化，其补给区与排泄区的分布位置有一定调整。这一变化与大沽河沿线分布的橡胶坝（图2）有关，雨季河道蓄水后，扰动了地下水流场，从而改变了补给和排泄的形态。

3.5 含水层边界识别与细化

图2、图3和图6中的外侧轮廓线为以往研究人员为大沽河含水层划定的边界，此边界是融合了钻孔数据、地质信息、水位统测和模型研究的综合成果，长期以来指导了大沽河含水层的各项水资源研究，具有重要的历史意义。但本次研究进一步改进和细化：

图6 大沽河含水层补给区、排泄区识别Fig.6 Recharge and discharge areas of the study site (first days of each month, 2013)

（1）含水层西北部存在较为明确的分水岭，见图5（b），其形态显著且稳定。在2013年全年未见其东西两侧的含水系统存在沟通，结合更长期的地下水水文监测数据判定后，考虑将含水层边界向东调整至分水岭处。

（2）含水层东北侧边界在历史研究中常考虑为隔水边界，但本研究显示其可能存在侧向径流补给。应考虑在此区域收集更多监测井数据，首先排除由监测点稀疏造成的潜在误判；若确有外源流入，应向东查找稳定的地下水分水岭，扩展含水层东侧边界，或者在模型研究中将其定义为侧向流入边界。

（3）大沽河上橡胶坝的蓄水状态对区域流场的扰动在雨季和旱季存在差别，后续研究中应对其进行重点考察。

3.6 方法局限与适用条件

本方法的应用前提是空间上高密度的水位监测点分布，但含水层空间异质程度与监测点密度需求之间尚无经验，亟待大量的实证研究，进行深入的探索和进一步细化。

本方法将地下水流场f(t,p)简化为梯度场I，并假定质点沿梯度最大的方向前进。这一假设是建立在的数值分布较为均匀的前提下，在渗透系数连续分布的单一含水层中大多可以满足。但若含水层渗透系数分布存在强烈的各向异性时（如裂隙含水层），这一前提将不复存在。假设2个相邻监测点分别属于2个垂向分隔的含水层，当其间存在水头差时，虽然水力梯度客观存在，但由于2点间渗透系数的分布存在强烈的各向异性（水平渗透系数远大于垂向渗透系数），其间不会存在真实的地下水流动。又如在隔水断层两侧水力梯度固然存在，但质点不会沿梯度运行到断层对侧，所以由此推断的补给区和排泄区有可能存在偏差。

如果含水层中存在渗透系数强烈各向异性的情形，其必然结果是在这类边界附近的质点运行速率（水力梯度）远大于含水层的平均值。EnviFusion-CGS中记录了所有位置的水力梯度和全部质点的运行速率，可以据此首先对渗透系数突变进行筛查，如果确实存在此情形，可以将含水层拆分后分别进行分析，或者与传统水文地质勘察和模型研究配合进行进一步的信息迭代。

4 结论

（1）高频率高密度水位监测数据是极佳的地下水系统信息来源。在不建立地下水数值模型的前提下，以监测井空间位置为节点，通过建立德劳内三角网格、求取网格中任意位置的水力梯度数值、借鉴使用机器学习领域中的优化算法，推断和识别区域内地下水补给和排泄边界。

（2）使用环境地学计算平台EnviFusion-CGS构建了识别研究区域补给区和排泄区空间信息的详细工作流程，帮助地下水科学家构建和修订已有的含水层概念模型。

（3）以山东省青岛市大沽河中下游含水层作为示范区，使用本文提出的方法对含水系统补给区和排泄区的空间分布及其动态变化进行分析，取得了良好效果。