煤炭工业是关系国家经济命脉的重要基础产业，支撑着国民经济持续快速健康发展^[1]。矿井水害是煤矿的主要灾害之一，是制约煤炭资源开发和影响煤炭企业安全的重要因素^[2]。随着中国煤炭开采深度和强度的不断增大，在开采过程中产生的剧烈扰动使地下涌水通道相互贯通，含水层的完整性遭到严重地破坏，增加了突水的危险性。一旦发生矿井突水，必须快速、准确地查明突水来源，为后续展开防治水工作提供理论依据和指导。由于地下水的水化学特征能够准确、有效地反映含水层的物化特性，所以，根据水化学特征来判断突水来源具有快速、准确、经济的特点^[3]。但地下涌水通道错综复杂，水化学特征随时间和空间会发生一定的变化，所以地下水的水化学特征表现出高度的非线性关系，采用传统的数学模型难以准确判别。

近年来，随着人工智能在数据挖掘方面的飞速发展，已被广泛应用于各科领域，尤以神经网络^[4-5]、支持向量机(SVM, support vector machine)^[6]、遗传算法(GA, genetic algorithm)^[6-7]和微粒群算法(PSO, particle swarm optimization)^[8]等模型为研究热点。目前，已有不少学者将其应用于矿井突水水源的判别，并取得了一些成果^{[6-7, 9]}。回声状态网络(ESN, echo state network)是近年来新兴的一种递归神经网络，独特而简单的训练方式以及高精度的训练结果已使其成为当前研究的热点之一^[10]。罗轶^[11]将ESN应用于交通流的预测并与Elman神经网络进行对比，结果表明：二者在预测精度相当的情况下，前者的训练速度较后者有了极大的提高。乔俊飞等^[12]采用PSO算法对ESN神经网络输出权值进行训练，并对污水的生化需氧量(BOD, biochemical oxygen demand)进行预测，结果证明了该模型能够有效预测污水处理关键水质参数。然而，ESN神经网络在矿井突水水源判别方面的研究较为匮乏。由于模型的判别精度一定程度上取决于样本的选取，而大多数有关水源判别的文献并没有对水源样本进行分析，所选取的待测样本特征较为明显，判别相对容易，虽然得到较高的判别精度，但不能充分检验模型的判别能力。因此，笔者对某矿含水层的水化学成分进行了统计分析和聚类分析，选取了具有代表性的样本作为训练样本，水化学特征与之近似的样本作为待测样本，将其载入基于阻尼最小二乘奇异分解法(DSVD, damped singular value decomposition)的GA_ESN判别模型，并与融合了全局动量因子和Dropout技术的GA_BP神经网络进行对比，以检验改进的GA_ESN模型的学习能力和预测能力，旨在针对水化学特征相近的水源，寻求一种能够快速有效判别突水来源的新方法。

1 样本的选取 1.1 样本特征分析

本文中所选取的水源样本来自文献[13]中某矿26个典型突水样本，该样本包含了11个水化学特征组分：Ca²⁺，Mg²⁺，K⁺，Na⁺，HCO₃^-，Cl^-，SO₄^2-，总硬度、碱度、pH、总矿化度(TDS, total degree of salinity)。其中1~6和X1~X3号样本为4-6煤顶板水，属于第1类水源；7~13和X4~X5号样本为奥灰水，属于第2类水源；14~20和X6号样本为6煤底至奥灰砂岩水，属于第3类水源。采用Rock Ware AqQa软件对该样本的6大常规离子和pH，TDS绘制了Durov图，如图 1所示。

图 1中2个三角区域为阴阳离子相对质量浓度分布，中间方形区域为6大离子质量浓度综合分布情况，2个矩形区域分别为pH和TDS分布。从方形区域可以看出，第2类水源分布较为集中，结合三角区域阴阳离子分布可初步得出该水源的水质类型为SO₄·Cl—K+Na型，从pH和TDS矩形区域可以看出7号样本出现偏离现象。第1，3类水源分布较为分散且相近，水质类型难以确定，通过观察可以发现4，17，18和X3，7，16号样本的空间分布较为相近，但它们分属于不同的水源类型，表明了其水质在一定程度上受到了外界的干扰，仅从空间分布角度分析，难以断定样本类型。

图 2是结合硬度和碱度特征分别对4，17，18和X3，7，16号样本绘制Schoeller图。从图 2中可以看出，4，17，18号样本的11个特征指标极为相近，尤其是HCO₃^-，Cl^-，总硬度、碱度、pH特征指标几乎没有差异；X3，7，16号样本相似度有所降低，但3条样本曲线均有不同程度的交叉重叠，因此采用上述统计的方法很难区分水源类型。为了进一步将样本的近似程度量化，采用模糊C均值聚类法对26个样本进行聚类分析。

1.2 FCM聚类分析

模糊C均值聚类^[14](FCM，fuzzy C-means)是用隶属度确定每个数据点属于某个聚类程度的一种聚类算法，其思想是使得被划分到同一类的对象之间相似度最大，而不同类之间相似度最小，是一种柔性的模糊聚类方法，其结果会给出每个样本点归属为某一类别的概率，通常选取概率较大者作为该样本的类别归属。

使用Matlab软件中fcm函数对26个样本数据进行FCM聚类分析，迭代次数为1 000，分类树为3类，其他参数采用默认值，计算得到的隶属矩阵如表 1所示，分类结果如图 3所示。

由图 3可以看出，FCM对第2类水源的划分与原始分类基本吻合，只有7号样本被划分为C1类，与1.1节的分析相一致。而第1，3类水源中有近50%的样本出现了与原始分类不相符的情况，表明了从空间距离角度分析样本的相似性与实际结果相违背。以4和17号样本为例，它们的隶属度最为接近，表明从空间距离角度来看，二者具有极高的相似性，但实际二者来自不同的含水层，出现这种差异，一方面表明了水源之间发生了一定的水力联系；另一方面表明采用图示法和聚类分析法无法找到样本与目标之间更为深层的映射关系。空间相似但分属不同含水层的样本是最能体现这种映射关系的代表性样本；同时，由于空间的相似性，使模型无法从距离相似这一特征来直接判别水源类型，增加了模型的判别难度。因此，应选取具有代表性的特征样本作为模型的训练样本，供模型学习其内在特征，选取与之空间近似样本作为待测样本，最大限度地检验模型的学习能力和预测能力。结合表 1和图 3的分析结果，在同一簇类中，选取不同突水来源但隶属度最为相近的样本作为训练样本，与之近似的样本作为待测样本。样本筛选结果见表 2。

2 改进的ESN神经网络 2.1 ESN神经网络

ESN是一种新型递归神经网络^[15]，其网络结构如图 4所示。

储备池状态更新方程：

$ \mathit{\boldsymbol{x}}\left( {\mathit{n}{\rm{ + 1}}} \right){\rm{ = }}{\mathit{f}_{{\rm{inter}}}}{\rm{[}}{\mathit{\boldsymbol{W}}_{{\rm{in}}}}\mathit{\boldsymbol{u}}\left( {\mathit{n}{\rm{ + 1}}} \right){\rm{ + }}\mathit{\boldsymbol{Wx}}\left( \mathit{n} \right){\rm{ + }}{\mathit{\boldsymbol{W}}_{{\rm{back}}}}\mathit{\boldsymbol{y}}\left( \mathit{n} \right){\rm{], }} $

(1)

式中：u(n)∈R^K为输入变量；K为样本特征维度；y(n)∈R^L为输出变量，L为样本输出维度；x(n)∈R^N为状态变量，N为储备池神经元个数；n为样本数量；f_inter为储备池激活函数；W_in∈R^K·N, W∈R^N·N, W_back∈R^N·L分别为输入层、中间层、反馈层连接权重。

收集每个样本[u(n)，x(n)]到状态矩阵M∈R^n·(K+N)，并计算与之对应的输出矩阵Y∈R^n·L, 最后根据式(2)计算出输出权重W_out∈R^(K+N)·L。

$ {\mathit{\boldsymbol{W}}_{{\rm{out}}}}{\rm{ = }}{\mathit{\boldsymbol{M}}^{\rm{ + }}}\mathit{\boldsymbol{Y}}{\rm{ = (}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{M}}{{\rm{)}}^{{\rm{ - 1}}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{Y}}{\rm{, }} $

(2)

式中：M⁺为矩阵M的伪逆。

式(2)为标准ESN计算W_out的方法，在实际应用中，M的法矩阵条件数通常比较大，直接采用伪逆法会产生不稳定解，将会导致待测样本的计算结果被放大。即模型学习精度高，预测能力差，出现过拟合。另外，由于矿区所保留的突水样本相对有限，因此所建立的模型参数不宜过多，否则会使过拟合问题更为严重。故引入了正则化方法来解决这一问题。

2.2 阻尼最小二乘奇异分解法

阻尼最小二乘法^[16]又称Levenberg-Marquardt法，是解算病态问题较为有效的一种正则化方法。马奎特提出给式(2)系数矩阵M^TM的对角矩阵加上一个正常数，以改善方程的条件，见式(3)。

$ {\mathit{\boldsymbol{W}}_{{\rm{out}}}}{\rm{ = (}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{M}}{\rm{ + }}\mathit{\alpha I}{{\rm{)}}^{{\rm{ - 1}}}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{Y}}。$

(3)

根据式(4)对M进行奇异值分解：

$ \mathit{\boldsymbol{M}}{\rm{ = }}\mathit{\boldsymbol{U \boldsymbol{\varSigma} }}{\mathit{\boldsymbol{V}}^{\rm{T}}}{\rm{ = }}\sum\limits_{\mathit{i}{\rm{ = 1}}}^\mathit{q} {{\mathit{\boldsymbol{u}}_\mathit{i}}{\mathit{\sigma }_\mathit{i}}\mathit{\boldsymbol{v}}_\mathit{i}^{\rm{T}}, } $

(4)

式中：U∈R^n·I和V∈R^(K+N)·I为正交阵；Σ=diag(σ₁, σ₂, …, σ_I)∈R^I·I, σ₁≥σ₂≥…≥σ_I＞0为矩阵M的奇异值；u_i, v_i分别表示U，V的列向量。

联立(3)、(4)，整理后可得：

$ {\mathit{\boldsymbol{W}}_{{\rm{out}}}}{\rm{ = }}\sum\limits_{\mathit{i}{\rm{ = 1}}}^\mathit{q} {{\mathit{f}_\mathit{i}}\frac{{\mathit{\boldsymbol{u}}_\mathit{i}^{\rm{T}}\mathit{\boldsymbol{Y}}}}{{{\mathit{\sigma }_\mathit{i}}}}{\mathit{\boldsymbol{v}}_\mathit{i}}, } $

(5)

式中：${\mathit{f}_\mathit{i}} = \frac{{{\mathit{\sigma }_\mathit{i}}}}{{{\mathit{\sigma }_\mathit{i}} + \mathit{\alpha }}}$。

由式(5)可以看出引入的阻尼因子α可以对奇异值的大小进行调整，来削弱奇异值对解的影响程度，从而减弱病态问题。

2.3 改进的ESN模型算法流程

改进的ESN模型搭建步骤如下：

Step 1  样本预处理。为了消除量纲、降低样本维度及噪声，将原始样本数据归一化至[0, 1]，然后对其进行主成分分析(PCA，principal component analysis)，并保留95%的原始信息。

Step 2  构建DSVD_ESN神经网络框架。采用DSVD正则化方法计算W_out，并由L曲线法确定阻尼因子α；预留ESN网络中需要寻优的关键参数接口(比例因子、激活函数、稀疏度D和储备池规模N)，以备遗传算法调用。

Step 3  遗传算法接口设定。初始化种群规模、最大迭代次数和遗传操作等参数；调整相关接口与DSVD_ESN模型进行对接。将K折交叉验证法嵌入GA算法中，以确保进化方向的正确性^[7]。

Step 4  将Step1得到的样本数据载入Step 2和Step 3所建立的GA_DSVD_ESN模型中，通过全局搜索得到最优的DSVD_ESN模型。此时，该模型可用于判别该矿区其他突水样本。

GA_DSVD_ESN模型流程如图 5所示。

3 实例分析

以表 2中的水源样本为例，将3种突水类型分别编码为“1”，“2”，“3”，并建立了4种水源判别模型，如表 3所示。

3.1 基于BP的判别模型

GA和BP算法参数设定见表 4和表 5。

由表 5中的网络结构可以看出，经过对样本归一化和PCA处理后，样本维度由11降为5，即模型输入维度为5；N_BP模型参数根据经验和反复测试得到；GA_N_BP模型参数通过GA算法对BP神经网络初始权值、动量因子、学习率和Dropout分数进行全局搜索得到。2种模型的判别结果见表 6。

从表 6可以看出，N_BP模型回判误判率为22.73%，将部分第3类水源错判为第1类水源，待测样本中，将7和18号样本错判为第1类水源；而GA_N_BP模型的预测效果有了较大的提高，回判误判率降低为4.55%，但仍对3，7号样本出现误判，由1.2节分析可以看出3号样本特征分布与17号样本极其相似，这是由于GA_N_BP模型没有很好学习4与17号样本间除空间相似之外其他内在特征，故模型的学习能力有待提高。

3.2 基于ESN的判别模型

采用GA算法分别对标准ESN和DSVD_ESN模型的关键参数进行全局寻优，结果见表 7所示。

由表 7可以看出，GA_ESN和GA_DSVD_ESN模型最大的不同在于网络结构，GA_DSVD_ESN模型储备池规模和稀疏度均明显小于GA_ESN模型，表明采用DSVD正则化方法可以有效的简化模型结构，减少隐层神经元个数和连接权值，从而降低了模型的连接关系；从模型所选取的激活函数可以看出，无论采用哪种权值训练方法对激活函数的选择没有影响，表明该激活函数能很好的表达这种输入—输出的映射关系。它们的判别结果如表 8所示。

从表 8可以看出，2种判别模型均能正确的判别突水来源，GA_DSVD_ESN模型的残差分布相对较低，维持在1E-2量级。

图 6为模型判别结果与目标值的对比曲线。可以看出，GA_DSVD_ESN模型曲线与目标曲线更为接近，上下浮动较小；而GA_ESN模型对7号样本的预测出现了较大的波动，拟合效果不如前者。

表 9是对2种模型判别结果的评估，不难看出，二者的回判正确率均达到了100%，但从回判结果的标准均方根误差(NRMSE, normalized root mean squared error)来看，GA_ESN模型达到了1E-14量级，GA_DSVD_ESN模型只有1E-2量级，这是由于GA_ESN模型采用伪逆法计算W_out，是一种无偏估计，导致模型的回判精度极高，但这种方法计算出的权值浮动较大，模型敏感度增强，从而导致模型的泛化能力降低。因此，该模型虽能正确预测突水来源，但其精度远不如前者，然而，GA_DSVD_ESN模型能将回判NRMSE控制在1E-2量级，且准确率为100%，表明了DSVD正则化方法能够将ESN模型的学习能力控制在有效范围内，以牺牲不必要的训练精度来提高模型的泛化能力，有效解决了过拟合问题；同时结合表 7和第1节的分析结果可以得出，在待测样本极为相似的情况下，GA_DSVD_ESN模型能够大幅简化储备池神经元个数及其连接权值，优化了网络的映射关系，提高了模型的训练速度，同时获得满意的预测效果，充分证明了该模型对于空间相似的水源样本具有较强的学习能力和预测能力。

4 结论

1) 采用Durov图和Schoeller图对某矿26个典型突水样本中的11种特征组分进行成分分析，结果发现：该矿区奥灰水含水层的水化学特征组分分布较为集中，而4-6煤顶板水和6煤底至奥灰水砂岩水特征组分分布较为分散且相似，给水源判别带来很大的困难。

2) 采用FCM法对水源样本进行聚类分析，得到了与结论1)一致的结果，并根据隶属矩阵确定了每个样本分属不同类别的概率，量化了样本间的相似程度。

3) 选取隶属度较高的样本作为待测样本，用以评估模型的判别能力。实验结果表明：融合动量因子和Dropout技术的N_BP模型判别效果最差；而GA_N_BP模型的判别效果较N_BP模型有了很大的提高，但仍不能满足实际要求；而GA_DSVD_ESN模型能够准确判别相似水源样本，较GA_ESN模型具有建模速度快，准确度高，泛化能力强的特点，同时该模型耦合了PCA主成分分析法和K折交叉验证法，使建模过程更具合理性，判别结果更具可靠性。因此，针对水化学特征相似的水源样本，GA_DSVD_ESN模型可作为一种快速有效判别方法。