基于人工神经网络的晋祠泉水位模拟研究

邢立文，崔宁博

（1.山西省水利水电科学研究院，山西 太原030000；2.四川大学水利水电学院，四川成都610065）

我国北方地区分布有大片碳酸盐岩岩溶地层，这些地区的岩溶地下水资源非常丰富。与南方岩溶地区的岩溶地下水不同，北方的岩溶大泉具有汇水面积广，径、补、蓄、排独立，泉水排泄集中，流量稳定且集中，水质优良，地下调节库容巨大的特点，多数具备建立大型供水水源地的有利条件，在供水方面有很大优越性，是北方很多城市及工农业供水水源地。山西、陕西、山东、河南4省岩溶地下水资源开采量已达24.39亿m3/a[1]。改革开放以来，我国北方地区城镇化进程发展迅速、工农业及第三产业高速发展，岩溶地下水资源过量开发和煤炭产业粗放式扩张，对岩溶地下水循环条件造成了严重破坏和影响，很多岩溶大泉出现断流、水位下降、污染等问题，甚至引发环境地质灾害。我国北方岩溶地下水资源迫切需要从单纯的开发利用向开发利用与保护管理并重的可持续发展转变。

人工神经网络（ANN）是由大量神经元按一定方式连接而成用来模仿人脑神经网络结构和功能建立的灵活而复杂的数学网络。神经网络算法具备线性和非线性学习算法的能力，根据已经收集的信息来传输，收集和学习信息。人工神经网络根据已知的输入和输出数据进行训练，可以预测未来的输入对应的输出。

1 研究区概况

山西是我国北方岩溶泉水最多的省份，岩溶泉水占地下水资源总量的33.52%[2]，为山西省能源基地的建设提供了有力保障。据统计，流量大于0.10 m3/s的岩溶泉水有83处，其中流量大于1.00 m3/s的有这18处，这18处泉水多年平均总流量为98.84 m3/s。岩溶地下水在山西省生产、生活和文化建设中具有重要支撑作用，但20世纪80年代后，受自然环境变化、人类活动等因素影响，岩溶地下水日趋严重的环境问题已成为山西省国民经济可持续发展的制约因素。晋祠泉属汾河流域，是太原市排名第二的岩溶大泉，20世纪50年代平均流量为1.99 m3/s，年际泉水流量变化幅度为0.20 m3/s左右，属稳定型泉水。晋祠泉具有保障太原市经济建设供水的重任[3]，但是大量开采地下岩溶水资源导致晋祠泉水位、流量大幅度衰减，年均泉水流量由20世纪60年代的1.69 m3/s下降至70年代的1.13 m3/s，再由80年代的0.46 m3/s下降至90年代初的0.15 m3/s，最终于1994年4月30日断流，严重影响当地工农业生产，同时使旅游景观大为逊色，并危及千年文物古迹，引起了有关部门的高度重视。

2 研究方法

2.1 泉水位影响因素及其计算方法

晋祠岩溶水系统补给来源为岩溶裸露区降水入渗补给和汾河在岩溶河段的河道渗漏补给。晋祠岩溶水系统排泄主要分为人类活动排泄和自然排泄两类，其中自然因素包括晋祠泉天然排泄、以潜流形式向平泉和太原盆地孔隙水排泄，人类活动排泄包括晋祠泉域内岩溶水开采和矿坑排水。

（1）降水入渗补给。降水入渗是晋祠泉域岩溶地下水系统的主要补给来源，晋祠泉域面积为2 030 km2，其中裸露可溶岩面积391 km2，泉域内其他地区入渗条件较差，所以晋祠泉域入渗补给只考虑裸露区降水入渗补给量。

式中：Q为降水入渗补给量，万m3/a；α为降水入渗补给系数；P为裸露区降水量，mm；F为裸露可溶岩面积，km2。

（2）渗漏补给。渗漏补给量采用河道退水系数法[4]计算：

式中：Q渗为区间河道渗漏量；m为河道渗漏系数；Q上为上游断面实测流量；η为河道退水系数，η=Q上/Q下，其中Q下为下游断面实测流量。

（3）岩溶水资源开发利用。晋祠泉域2012年有149口岩溶水井，取水量为1 605万m3[5]。根据晋祠泉域岩溶水开发利用情况，1994—2012年岩溶地下水开采量见表1。

表1 晋祠泉域1994—2012年地下水位与各均衡要素汇总

（4）煤矿矿井排水量。1994年泉域排泄区岩溶水的矿井排水量为0.14 m3/s，随着煤层开采深度的增加排水量逐渐增大，2006年达到最大值[2]。2008年山西省煤矿兼并重组后，矿井排水量逐年减小，见表1。

（5）侧向排泄量分析。侧向排泄量是泉域岩溶地下水循环的一个重要因素。根据晋祠泉域的水文地质条件，除了以泉水形式排泄的岩溶水外，还有部分岩溶水横向排泄变为盆地孔隙水和平泉。晋祠泉域侧向排泄量与水力梯度和断面渗透系数等有关，因此采用Darcy公式计算侧向排泄量，结果见表1。

为了晋祠泉水复流和改善晋祠泉生态环境，从2013年开始，山西省委、省政府及有关部门积极采取了多种有效治理措施，诸如关停小煤矿和实施煤矿兼并重组，减少采矿排水；水源置换和关井压采；汾河清水复流等。晋祠泉水复流工程实施前后地下水要素见表2。

表2 晋祠泉2013—2017年地下水位与各均衡要素汇总

2.2 泉水水位预测模型

基于前馈神经网络Feed-forward-net、动态递归神经网络Elman-net、时延神经网络Time-delay-net、非线性动态自回归神经网络Narx-net、级联神经网络Cascade-net共5种人工神经网络构建晋祠泉水位预测模型。

（1）前馈神经网络。前馈神经网络是一种方便、直观、有效的基于误差反向传播算法的比较成熟的人工神经网络模型。在这种网络中，学习过程由正向和反向传播组成。在正向传播过程中，输入信号从输入层经隐含层单元传向输出层，每一层神经元的状态只影响下一层神经元的状态。若在输出层不能得到期望输出时，则转入反向传播，将输出信号的误差沿原来的连接通路返回，通过修改各层神经元的权值和阈值，使得网络全局误差最小[6]。常见的3层前馈神经网络模型结构见图1（wij为输入层第i个神经元与隐含层第j个神经元之间的连接权值，vjt为隐含层第j个神经元与输出层第t个神经元之间的连接权值）。

图1 前馈神经网络结构

设隐含层第j个神经元的阈值为θj，输出层第t个神经元的阈值为γt，则隐含层第j个神经元的输入sj、输出bj为

式中：xi为输入层第i个神经元的输入；n为输入层神经元个数；f（）为模拟神经元非线性特性的激励函数。

输出层第t个神经元的输入Lt、输出Ct为

式中：p为隐含层神经元个数。

（2）动态递归神经网络。动态递归神经网络作为一阶时滞算子，达到了记忆的目的[7]，因此具有适应时变特性的能力，增强了网络和网络本身的全局稳定性。其处理动态信息的能力强于前馈神经网络，可以用来解决快速优化的问题。动态递归神经网络结构见图2。

图2 动态递归神经网络结构

动态递归神经网络的数学表达式：

式中：y（t）为 m 维输出节点向量；x（t）为 n维中间层节点单元向量；u（t）为 r维输入向量；xζ（t）为 n 维反馈状态向量；W3为隐含层到输出层的连接权值；W2为输入层到隐含层的连接权值；W1为承接层到隐含层的连接权值；g[]为输出神经元的传递函数，是隐含层输出的线性组合；f｛｝为隐含层神经元的传递函数。

（3）时延神经网络。时延神经网络在结构上与普通神经网络相似，但是每个神经元可以处理一个输入子集，并且对于具有不同时滞的输入数据具有不同的权重[8]。时延神经元TDN是时延神经网络的基本单元。 每一个 TDN 有 M 个输入 I1（t）、I2（t）、…、IM（t）和一个输出O（t）。每个输入Ii（t）包含N步时延，用以保存当前时刻之前的N个时刻的输入Ii（t-d）（d=1，

式中：bi为第i个输入的偏置；f｛｝为激励函数。

时延神经网络结构见图3，其中隐含层有J个神经元，输出层有R个神经元。输出层各个神经元输出为

式中：N1为输出层节点的时延步数；r为输出层节点编号；为第r个输出层节点与隐含层节点间的连接权重；为第r个输出层节点的偏置；（t）为第 j个隐含层的输出；N2为隐含层节点的时延步数；j为隐含层节点编号；Ii（t-d）为第 i个输入层节点；为第j个隐含层节点与输入层节点的连接权重；为第j个隐含层节点的偏置。

图3 时延神经网络结构

（4）非线性动态自回归神经网络。Narx神经网络将前馈神经网络的输出向量引入到输入向量中，在保持时延后通过外部反馈。本文采用Narx神经网络代替三层前馈神经网络进行长期预测，给出了一种非线性动态自回归神经网络，Narx神经网络的结构见图4。其中：R=nu+ny，表示两层Narx神经网络的输入个数，nu为非线性系统的输入阶数，ny为非线性系统的输出阶数，nu、ny的值可以通过模型阶次辨识方法估计得出；j（图4 中 i=1）为 Narx神经网络隐含层中第 i个神经元和输入向量中第j个元素之间的连接权值；S1为Narx神经网络隐含层中神经元的个数；为Narx神经网络隐含层中第i个神经元的偏置；为Narx神经网络隐含层中第i个神经元的净输入；f1为Narx神经网络隐含层神经元的传输函数；为Narx神经网络隐含层中第i个神经元的输出值；为 Narx 神经网络的输出神经元和隐含层中第i个神经元之间的连接权值；为Narx神经网络输出神经元的偏置；为Narx神经网络输出神经元的净输入；f2为Narx神经网络输出神经元的传输函数，在此使用纯线性函数；既表示Narx神经网络输出神经元的输出值，又表示Narx神经网络在k+1时刻的预测输出值ym（k+1）。

图4 非线性动态自回归神经网络结构

（5）级联神经网络。在分类和预测方法的基础上，选择BP神经网络作为级联神经网络的前级网络，RBF径向基网络作为后级网络。充分利用BP神经网络学习时间短、RBF网络训练速度快的优点，建立了基于BP神经网络和RBF网络的级联神经网络预测模型。将自然因素作为BP神经网络的输入变量，利用自然因素预测BP神经网络的输出。以BP神经网络的输出和人类因素为输入，通过RBF网络输出自然和人为联合作用下的泉水位预测值。该方法简化了BP神经网络和RBF网络的拓扑结构，同时可以取得满意的预测精度和计算效率。因篇幅有限，不再赘述RBF网络的原理和结构，级联神经网络结构见图5。

图5 级联神经网络结构

2.3 降水量预测模型

LSTM深度学习模型是Hochreiter＆Schmidhuber为了避免长依赖问题结合适当的梯度学习算法设计的一种特殊的循环神经网络（RNN），该模型具有删除或添加信息到元胞（cell）状态的门（Gate），可以学习长期依赖的历史信息。LSTM深度学习模型与RNN模型最主要的区别在于其算法中加入了一个判断信息有用与否的元胞结构，而LSTM的关键正是该元胞状态——元胞的状态类似于输送带，元胞的状态在整个链上运行，只有一些小的线性操作作用其上，信息很容易保持不变的流过整个链，即激活函数的结构均为线性的，LSTM神经网络只是对其进行传递，没有做复杂的变换，不会更改信息。LSTM深度学习模型的每个判断元胞当中放置了4扇阀门，分别为输入阀门、遗忘阀门、更新阀门和输出阀门，这些阀门可以打开或关闭，用于判断模型网络的记忆态在该层输出的结果是否达到阈值从而加入到当前该层的计算中。当一个输入信息进入LSTM的判断元胞当中，该元胞即可根据既定规则来判断该值是否有用。只有符合算法认证的信息才会留下并被记忆，不符合的信息则通过遗忘阀门被遗忘，从而有效地解决循环神经网络（RNN）长期依赖的问题。LSTM深度学习模型结构见图6。

图6 LSTM深度学习模型元胞结构示意

（1）遗忘阀门Ft。遗忘阀门通过sigmoid激活函数读取旧元胞的输出值Yt-1和当前元胞的输入值Xt，并结合偏置矢量和相应的权重来判断是否保留当前元胞的状态 St，即“1”表示“保留“，“0”表示“删除”。 该阀门的计算原理：

式中：σ为sigmoid激活函数；Wf为遗忘阀门的权重；Yt-1为上一个元胞的输出值；Xt为当前元胞的输入值；bf为遗忘阀门的偏置矢量。

（2）输入阀门It。输入阀门用来判断输入值Xt是否应该用于更新当前元胞的状态St。输入阀门一方面由sigmoid激活函数负责决定当前元胞的输出值Yt，另一方面由tanh激活函数负责创建一个当前元胞的候选矢量Ut，该矢量由更新阀门判断是否应将其加入到当前元胞的状态St中。该阀门的计算原理：

式中：Wt为输入阀门的权重；Wc为候选矢量的权重；bt为输入阀门的偏置矢量；bc为候选矢量的偏置矢量；tanh表示双曲正切激活函数。

（3）更新阀门Ut。该阀门通过将遗忘阀门和输入阀门相乘来判断当前元胞是否应该将旧元胞的状态由St-1更新到St。该阀门的原理：

式中：St-1为旧元胞的状态。

（4）输出阀门Ot。输出阀门负责计算当前元胞的输出值Yt。输出阀门通过运行sigmoid激活函数，判断旧元胞的输出值和当前元胞的输入值来确定是否输出当前的元胞状态。接着把当前元胞状态St通过tanh激活函数进行处理，并将它和sigmoid激活函数的输出Ot相乘，最终确定当前元胞的输出值Yt。输出阀门的原理：

式中：Wo为输出阀门的权重；bo为输出阀门的偏置矢量。

2.4 模型训练算法

晋祠泉水位预测模型训练算法选择了traingd（梯度下降法）、traingdm（有动量的梯度下降法）、traingda（自适应lr梯度下降法）、traingdx（自适应lr动量梯度下降法）、trainrp（弹性梯度下降法）、trainbr（贝叶斯正则化算法）、traincgf（Fletcher-Reeves共轭梯度法）、traincgp（Ploak-Ribiere共轭梯度法）、traincgb（Powell-Beale共轭梯度法）、trainscg（量化共轭梯度法）、trainbfg（拟牛顿算法）、trainosss（一步正割算法）、trainlm（Levenberg-Marquardt莱文贝格-马夸特算法）、trainr（学习功能的随机顺序增量算法）共14种算法对5种人工神经网络进行训练。

2.5 模型精度评价

晋祠泉水位模型精度采用均方根误差（RMSE）、效率系数（NSE）、平均绝对误差（MAE）、确定系数（R2）、整体评价指标（Gpi）[9]进行评价，其计算方法：

式中：Xi、Yi分别为第i个时段泉水位的实测值和估算值；X′、Y′分别为泉水位实测值和估算值的平均值；n为时段数；Ti为 RMSE、MAE、R2、NSE 归一化后的值；T′i为对应参数的中位数；当Tj为RMSE和 MAE时，αj取-1，其他情况取1。

3 结果与讨论

3.1 训练算法效率

通过5种神经网络类型和14种训练算法组合，建立了70种不同的网络模型来预测晋祠泉水位。通过分析发现：在晋祠泉水位估算中，traincgb、trainrp、traincgf、traincgp等训练算法效果比较理想；trainbfg、traingd、trainlm、trainr等算法无法模拟得到晋祠泉水位；其他算法大多能有效估算晋祠泉水位，但存在训练次数多、隐含层节点数多、占用CPU大等问题，见表3。

表3 神经网络类型、训练算法运行效率比较

3.2 模型精度

基于5种人工神经网络构建了晋祠泉水位预测模型，采用1994—2012年实测降水入渗补给量、河道渗漏补给量、岩溶水开采量、煤矿矿井排水量、侧向排泄量和岩溶水位数据为训练集对神经网络模型进行训练，以2013—2017年为神经网络模型测试集验证模型的预测精度，预测精度指标包括效率系数NSE（该值越接近1，表明模型可信度越高）、均方根误差RMSE（该值越小，表明模型偏差越小）、平均绝对误差MAE（该值越小，表明模型精度越好）、确定系数R2（该值越接近1，表明模型吻合度越高）。为了统一模型精度评价体系，引入整体评价指标（Gpi）评价模型适用性，结果见表4、表5。

表4 人工神经网络泉水位模拟结果预测模型适应性评价

表5 人工神经网络泉水位预测模型适应性评价

由表5可以看出，本文选用的5种人工神经网络模型的各项适应性评价指标差异明显：NSE、RMSE值动态递归神经网络模型最优，前馈神经网络模型最差；MAE值级联神经网络模型最优，前馈神经网络模型最差；R2值前馈神经网络模型最优，非线性动态自回归神经网络模型最差。模型的适应性排名并不一致，无法单靠某一指标判断模型的精度和适应性，从而选择最合适的泉水位预测模型。因此引入了整体评价指标（Gpi）来综合评价模型的适用性，5种人工神经网络模型的Gpi值动态递归神经网络模型＞级联神经网络模型＞时延神经网络模型＞非线性动态自回归神经网络模型＞前馈神经网络模型，就此可以判断动态递归神经网络模型为预测晋祠泉水位变化的最合适模型。

由图7可以看出，晋祠泉在1994—2012年的近20 a间水位呈现谷状曲线，在2008年左右晋祠泉水位达到最低值774.94 m后，2009—2012年呈现逐渐上升趋势，最高值达到788.21 m。降水入渗量在1994—2008年虽然变化幅度不大但仍呈现整体上升的态势，这与晋祠泉水位下降的趋势相反，可见降水入渗并非影响晋祠泉水位变化的关键因素。同时，晋祠泉开采量呈现1994—2008年逐年增加继而到2012年逐年减少的态势，这与晋祠泉水位变化一致，可以判断出晋祠泉开采是晋祠泉水位变化的关键影响因素。

图7 晋祠泉域水位人工神经网络模型训练、验证及预测

2013年后，为晋祠泉水复流和改善晋祠泉生态环境，山西省政府采取了多种积极有效的措施，即使2013—2017年降水量下降明显，晋祠泉水位仍然保持了逐年上升的良好态势，这也印证了晋祠泉开采等人类因素才是影响泉水位变化的关键因素。

3.3 未来10 a晋祠泉水位预测

为了预测2018—2027年晋祠泉域地下水位变化趋势，采用动态递归神经网络模型预测晋祠泉域2018—2027年地下水位的变化趋势并测定了晋祠泉的复流时间（见图7）。其中基于晋祠泉域1956—2017年共62 a的降水系列资料，采用LSTM人工神经网络模型预报未来2018—2027年（未来10 a）晋祠泉域的降水量（见图7），并根据所预测的降水量采用式（1）计算降水入渗补给量，采用式（2）计算河道渗漏补给量，采用Darcy公式计算侧向排泄量，岩溶水资源开发利用量和煤矿矿井排水量则采用《太原市晋祠泉域水资源保护条例》中的限定值。

从图7可以看出，由于人类积极影响因素的介入，晋祠泉水位保持平稳上升的良好态势，泉水位在2019年可上升到802.66 m，超过了晋祠泉复流的最低水位802.59 m[10]，泉水位在未来10 a最高可上升至803.74 m，表明多种积极有效措施的实施可以取得良好的泉水复流效果。

4 结 论

人类因素是影响晋祠泉水位的关键因素，人类消极干预严重是造成晋祠泉水位下降、断流的主要原因。在人类积极干预介入下，泉水位完全具备复流的条件，晋祠泉在2018年左右可以复流。各种训练算法并非均能有效适用于晋祠泉水位预测，其中 trainbfg、traingd、trainlm、trainr等算法无法模拟晋祠泉水位，traincgb、trainrp、traincgf、traincgp 等训练算法效果较好，其他算法大多存在训练次数多、占用CPU大等问题。动态递归神经网络模型为预测晋祠泉水位变化的最合适模型。