查恩爽，肖 霄

（1.吉林大学建设工程学院，吉林 长春 130026；2.河北地质大学水资源与环境学院，河北 石家庄 050031）

吉林省西部地区地处松嫩平原西南部，是我国重要粮食产区。受半干旱大陆性季风气候控制，区内年均降水量约为390 mm，年内分布不均，主要集中于5—7月，春季具有十年九旱的特点。区内灌溉用水主要取自地下潜水及承压含水层[1-2]。为了满足农业灌溉的需要，戴长雷等[3]提出通过修建地下水库以丰补欠，实现区内水资源的人工调蓄。然而，吉林省西部地势低缓，地质历史时期潜水水位埋藏较浅，蒸发强烈，加之嫩江将大兴安岭地区大量矿物携带至松嫩平原西南部沉积，导致严重的土壤次生盐碱化[4]。为了改善土壤质量，区内开始实施旱田改水田工程[5]，即通过大水灌溉压盐的物理治理方式，降低表层土壤含盐量。但旱改水工程的实施，导致地下水位抬升。在停灌期间，受区内强蒸发条件的影响，可能加剧土壤次生盐碱化。综合考虑区内水资源调配及土壤次生盐碱化生态环境问题，吉林省西部地区开始实施节水灌溉[6]。

近年来，范思雨等[7]基于1980—2015年吉林省西部土地利用数据，论述了吉林省西部地区水田演化趋势以及对应的盐碱土分布范围变化。结果显示，旱改水工程对盐碱土的演化具有区域性的抑制作用，但在工程外围，土壤盐碱化程度呈现加剧的趋势[8]。杨沁瑜等[9]基于蚁群模型分析了吉林省西部土地利用安全格局。才源等[10]分析了西部盐碱土近年空间分布状态。卞建民等[11-12]系统研究了冻融影响下盐碱土中盐分迁移转化机理；张瑜等[13-14]对比分析了不同改良剂对盐碱土肥力改善效果。尽管关于吉林省西部土壤次生盐碱化的机理和改良的研究较多，但缺少不同方案实施后的效果预测与对比。

本文同时考虑水资源和土壤生态两方面因素，通过随机模拟的方法，对比分析目前区内典型的水资源利用方案，包括人工调蓄、旱改水工程、节水灌溉等实施5年后，潜水资源与生态环境的风险。此次研究建立的神经网络—随机模拟方法体系，更好地预测地下水位变化驱动因子，并遍历性考虑不同驱动因子影响下地下水中长期变化趋势，为我国干旱半干旱地区地下水资源可持续利用和管理提供了新方法。

1 研究区概况

选取吉林省西部相对独立的水文地质单元洮儿河流域山前倾斜平原作为研究区 （图1），面积为1 960 km2。区内第四系潜水主流向为自西北向东南，含水层厚度小于50 m，沿地下水流向含水层厚度逐渐减小[15-16]。含水层上游主要由中粗砂和砂砾石构成，渗透性好；下游由中细砂构成并含黏土夹层，渗透性变差。总体而言，潜水资源较为丰富，5 m 降深条件下单井涌水量可达1 000 m3/d[17-18]。

图1 研究区位置图、边界条件以及参数分区Fig.1 Location of the study area,boundary conditions and parameter zonation

区内潜水含水层接受大气降水、河道、灌溉回渗入渗补给以及侧向径流补给。排泄以灌溉开采为主，伴有蒸发和侧向径流排泄。区内洮儿河上游布设镇西和务本两处水文站，下游布设洮南水文站，记录了近50年洮儿河的流量和当地大气降水量。据此估算大气降水入渗量、河道入渗量总和超过5×108m3/a。2014—2018年，区内平均农业灌溉用水量约为 4×108m3/a，接近大气降水入渗补给量和河道入渗补给量的总和。

2 潜水资源与环境风险等级划分

潜水资源风险根据区内潜水含水介质厚度确定。含水介质厚度是指由潜水含水层隔水底板算起至地表之间的距离，是一个静止变量，代表了地下储水空间。含水介质厚度有别于含水层厚度 （隔水底板至潜水面的距离）。研究区内潜水含水介质厚度1～45 m （图2）。75%的区域含水介质厚度小于12 m，50%的区域含水介质厚度小于5 m。若未来5年内潜水水位埋深超过12 m，意味着超过75%的区域存在含水层疏干的风险，因此将水位埋深超过12 m 定义为高资源风险警戒水位埋深 （WRII）；水位埋深5～12 m 定义为中度资源风险水位埋深 （WRI）。当水位埋深过浅，小于1 m时，容易导致次生盐碱化[19]。因此，当水位埋深处于1 m 以浅 （包括埋深为负值，代表水位高于地表），认为区内潜水存在高环境风险 （ERI）。综上，水位埋深1～5 m 视为低风险，作为该区域潜水资源开发利用的水位控制标准。

图2 基于含水介质厚度概率分布的资源风险分级Fig.2 Resource risk classification based on probabilistic distribution of aquifer thickness

3 循环神经网络及地下水位变化天然驱动因子预测

河道入渗补给和大气降水入渗补给是该地区潜水含水层水位变化的主要驱动力。基于已掌握的河道入渗量及大气降水入渗量历史时间序列，合理预测未来5年的变化趋势，是准确模拟地下水位动态变化规律的关键。目前关于降水量、入渗量等时间序列数据拓展的方法有很多：统计学方法、物理过程模拟方法以及人工智能神经网络方法等[19-23]。不同分析方法各有利弊，统计学方法 （如回归分析、模糊综合评判、小波分析等）通常要求时间序列统计学关系相对稳定、服从正态分布或均匀分布等特定概率分布，难以处理异常气候变化过程，如干旱或暴雨等。物理过程模拟方法充分考虑降雨等变量的影响因素，如气温、下垫面、大气环流等。理论上该方法是一种准确的模拟方法，但在现实应用过程中，由于气候变量和影响因素过于复杂，难以量化，因此物理过程模拟结果准确度也难以保证。人工智能神经网络方法 （如经典的BP 神经网络）早在20世纪90年代已被用于气候因子时间序列预测，但由于过度拟合，一度被摒弃。近年来，人工智能算法飞速发展，出现诸多优秀的神经网络结构及训练算法，如卷积神经网络、图神经网络、循环神经网络等。其中，循环神经网络的专长是寻找不同周期历史数据的关联性，并预测数据的未来演化趋势，适合于时间序列的数据处理和信息挖掘。近年来，该方法已成功应用于降水量时间序列预测[19]。

本文采用循环神经网络方法 （RNN，Recurrent Neural Network），基于1970—2018年间大气降水量及地表水径流量的实测水文时间序列，分别推测2019—2023年大气降水量和地表径流量数值。RNN 计算方程如下[24]：

式中：Xt—输入数据；

Yt—输出数据；

Ht—中间过渡层；

W、U、B—权重矩阵。

此次循环神经网络输入和输出层之间设计3个隐层。输入层和前两个隐层的输出结果经过Sigmoid 激活函数处理：

经过激活函数处理后，可保证输出结果数值介于0～1 之间，利于不同层位权重的公平、稳定训练。

选择1970—2010年水文变量序列作为输入数据，通过式 （1）和（2）循环迭代计算确定2011—2018年水文变量，并与实测值进行对比，根据计算误差采用ADAM 算法 （Adaptive Moment Estimation Algorithm）[25]更新权重矩阵，直到计算值与实测值拟合为止。1970—2010年为神经网络训练阶段，2011—2018年为神经网络权重验证阶段。计算值与实测值的拟合程度采用均方差损失函数表示：

式中：t—时间，t=1～8分别对应2011—2018年；

N—总预测年数，N=8；

~Yt—第t年水文变量实测值；

Yt—第t年水文变量计算值。

确定权重矩阵后，通过式 （1）和式（2）预测2019—2023年大气降水量及镇西站、务本站和洮南站河流径流量。通过降水入渗系数法，根据大气降水量预测值计算获得大气降水入渗补给量。以镇西站和务本站河流径流量之和作为所选取的水文地质单元上游地表径流汇入量，以洮南站径流量作为下游地表径流流出量，二者之差作为河道入渗对地下水补给量。以年为时间单位用RNN 计算各水文站径流量。计算河道入渗补给量及大气降水入渗补给量时，忽略水文站之间径流过程的时间滞后性以及大气降水量和地表径流量之间的关联。

RNN 训练阶段，大气降水量计算误差＜5 mm（相对误差5%），各水文站年径流量计算误差低于104m3（相对误差＜1‰），8年验证期年降水量计算误差为3～20 mm（相对误差小于＜10%）、年地表径流量计算误差为1 320 ～3.2×104m3（相对误差小于＜5%）。

由图3可知，未来5年区内降水量370～450 mm/a，平均值为410 mm/a。大气降水入渗系数按0.23（已扣除蒸发损失）计算，可获得区内未来5年大气降水入渗补给量为1.5×108～ 2.0×108m3/a；而河道入渗补给量波动于3.2×108～ 6.0×108m3/a。按2014—2018年灌溉用水量估计，未来5年灌溉用水量（潜水）波动于3.8×108～ 6.9×108m3/a，灌溉回渗量为4.0×107～1.0×108m3/a。

4 数值模拟与潜水资源利用方案讨论

潜水水位埋深的预测基于二维多孔介质非稳定水流运动方程[26]：

式中：T—导水系数/ （m2·d-1）；

h—潜水水位标高 /m；

b—含水层底板标高 /m；

x，y—空间坐标 /m；

光镜可见青年对照组中大鼠心肌纤维排列整齐紧凑，结构清晰，心肌细胞形态形态、胞质及间质均正常；与青年对照组相比，自然衰老模型组中大鼠心肌细胞排列紊乱，细胞间隙增大，心肌细胞变性增多，间质纤维增生，肌纤维断裂。与 24月龄相比，TSPJ低、高剂量组大鼠心肌纤维形态得到明显改善，纤维排列整齐，心肌细胞炎性浸润减轻，间质纤维增生减少。结果如图1所示。

图3 年降水量和年河道入渗补给量实测值和计算值图Fig.3 Average annual rainfall and annual average river infiltration recharge

μ—给水度；

Q—单位面积补给量（包括大气降水入渗补给量、河道入渗补给量和灌溉回渗量）或排泄量（开采量和蒸发量）/（m·d-1）。

在研究区进行网格离散后（此次剖分为10 000个三角形单元格），应用水流运动方程，采用FEFLOW求解[27]。

根据区内多年地下水位等值线和流场分布，洮儿河流域外围以地下水流线作为隔水边界。上游两处断面为流入边界，下游一处断面为流出边界 （图1）。根据多年平均水力梯度，通过达西定律计算得到地下水侧向流入量约为1.0×107～5.0×107m3/a，侧向流出量约为3.0×107～ 4.0×107m3/a。

参数分区根据陈社明[28]对吉林省西部地区数值模型识别和验证结果确定，并结合潜水含水层岩性特征，定义渗透系数和给水度置信区间。由于研究区潜水含水层为山前洪积扇相沉积、下游接续辫状河湖相沉积，因此假设渗透系数和给水度分别服从对数正态分布和正态分布[29]。不同参数分区，参数取值见表1[28]。

模拟首先设置初始水位为10 m，通过FEFLOW 模拟1970—2018年水头变化，获取2018年地下水位空间分布。获得的2018年年底的地下水位作为2019—2023年水位预测的初始水位。

此次研究对比6种不同水资源开发利用模式下（表2），研究区内地下水位空间分布规律及风险等级。根据RNN 计算结果和目前研究区内灌溉用水状况，确定维持现状条件下，区内补给量未来5年为4.7×108～8.0×108m3/a，开 采 量3.8×108～ 6.9×108m3/a，假设二者服从均匀分布，作为此次模拟的情形1。情形2 假设干旱条件下，开采量与现状条件相同，但区内补给量由现状条件降低至0～4.0×108m3/a。情形3 假设区内补给量与现状条件相同，开采量的均值与现状条件相同，均为4.5×108m3/a，但由于无序开采导致开采量波动幅度增加一倍，开采量取值4.0×108～11.0×108m3/a。情形4 考虑地下水库调蓄功能，使得区内年均补给量和开采量不变，但以丰补欠，旱季增加开采，雨季降低开采并增加补给；情形5 在节水灌溉条件下，区内开采量降低至2.0×108～ 3.0×108m3/a，对应补给量随着灌溉回渗量减小略有下降；情形6是在旱田改水田工程实施条件下，开采量维持不变（灌溉用水主要取自地表水和第四系承压水），区内补给量增加50%。

表1 不同分区内渗透系数、给水度概率分布[28]Table1 Probability distribution of coefficient of permeability and specific yield in different zones

表2 不同工况条件下研究区域补给量与开采量取值范围Table2 Ranges of recharge and exploration rate under different scenarios

根据表1和表2中的参数取值范围以及所服从的概率分布，采用蒙特卡洛方法随机选取区域补给量、开采量、渗透系数及给水度。对应不同参数取值，分别进行100次水流模拟，确定5年后的水位埋深空间分布。取水位埋深平均值，进行潜水资源与环境风险分析。如图4（a）所示，在维持现状条件下（情形1），区内地下水位埋深（5～12 m）普遍存在中度资源风险（WRI），特别是洮儿河下游约30 km2范围内地下水位埋深超过12 m，存在高资源风险 （WRII），现状开采条件可能导致潜水含水层疏干。同时，局部范围潜水水位埋深低于1m，存在较高的土壤盐渍化风险 （ERI）。在多年干旱条件下（情形2），研究区内水位埋深普遍增加，存在高潜水资源风险的区域面积较情形1 条件下增加将近1倍（图4b）。若区内继续无序开采（情形3），则导致区内潜水风险加剧，将近一半区域存在高度资源风险（图4c）。

图4 不同方案条件下水位埋深分布图及风险等级区划Fig.4 Burial depth of water table and risk classification under six different scenarios

后三种不同水资源规划方案的实施，将不同程度地缓解区内潜水资源风险。如图4（d）所示，修建地下水库开展人工调蓄 （情形4），洮儿河高度资源风险区域面积降低至14 km2，但人工调蓄条件下，由于开采与补给具有反相关特点，旱季开采量较大、而丰水期未能得到充分的补给，导致研究区的西南部出现了水资源高风险区域。在现状条件下，采取节水灌溉措施可使区内潜水资源风险和环境风险程度普遍降低，全区仅有不到10 km2的区域存在高度资源风险 （图4e）。而旱田改水田工程的实施（图4f），可使区内潜水水位明显抬升，部分位置甚至出现地下水补给河水（水位埋深低于0 m）。尽管旱改水工程有效抑制了水资源风险，但沿河两岸存在次生盐碱化生态环境风险区域面积超过700 km2。

不同情形下潜水水位埋深的概率分布显示，潜水资源无序开采时 （图5a），潜水水位埋深较少处于低风险区域，甚至超过40%的潜水水位埋深处于高度资源风险状态。因此认为：地下水资源的无序开采是导致研究区潜水资源枯竭最主要的因素。采取水资源管理措施后 （图5b），可使区内地下水位埋深普遍减小，即旱田改水田 （情形6）条件下水位恢复效果最为明显，但却普遍存在高环境风险区域。相比而言，节水灌溉可使潜水水位埋深较大概率存在于低风险区域，兼顾资源需求和环境效应，较其他方法优势明显。

图5 不同方案条件下地下水位埋深概率分布（水位埋深为负值代表水位高于地面）Fig.5 Probabilistic distribution of the buried depth of groundwater table under six scenarios (Negative values indicates the water level higher than the land surface)

5 结论

此次研究以洮儿河流域为例，定义了潜水资源风险和生态环境风险，并通过随机模拟的方法，对比3种现状条件、3种水资源管理措施条件下区内潜水资源与环境风险，得到主要结论：

（1）潜水资源的无序开采，是导致研究区资源风险加剧的主要诱因。

（2）地下水库人工调蓄、旱田改水田、节水灌溉均可使区内潜水水位得到明显恢复，其中旱田改水田工程实施后潜水水位恢复效果最为明显，同时也加剧了次生盐碱化造成的生态环境风险；节水灌溉在满足区内灌溉需求的前提下，可将区内净开采强度降低2.0×108～3.0×108m3/a，保障区内潜水水位埋深处于合理范围之内 （1～12 m 之间）。从潜水资源风险与生态环境风险两方面考虑，节水灌溉是目前吉林省西部水资源利用最有效的管理方法。

（3）循环神经网络成功预测降水和河流入渗补给量，确定了中长期地下水位变化驱动因子；随机过程模拟计算中长期地下水位变化的同时，可明确地下水位变化的不确定性。神经网络和随机过程模拟相结合，提高了地下水位预测的可靠度，为干旱半干旱地区地下水资源利用与管理提供了新的方法。