煤矿掘进工作面是主要的透水地点，为防止掘进期间工作面透水事故发生，需要对巷道前方隐伏水源和导水通道进行超前探测^[1-2]。目前矿井掘进工作面水害隐患超前探测中应用最多的是瞬变电磁法和直流电阻率法等对低阻体更为敏感的电法类探测方法^[3-6]。其中，直流电阻率法探测因其施工便捷，在矿井掘进工作面水害超前预测中得到广泛应用。

20世纪末，直流电阻率法开始用于巷道掘进工作面超前探测^[2]，最初采用单电极供电观测系统，后来发展成为三电极供电观测系统^[7-8]，并形成了以几何交汇原理为基础三电极法探测信号处理和解释方法。此后多位学者通过工程实例及其揭露验证情况对直流电阻率法超前探测方法进行了有效性分析，认为矿井直流电阻率法超前探测法可应用于水害超前预测^[5,9]。

随着应用的不断深入，对矿井直流电阻率法超前探测精度提出了更高的要求。诸多学者对该技术进行了更为深入的研究，认为直流电阻率法超前探测的工业技术研究超前于理论研究，虽然在工业应用中取得了诸多良好的效果，但是其带来一个严重的问题，即该技术的基础理论研究明显滞后于工业方法技术研究，对提升探测精度带来不利影响^[10-15]。其中李飞等^[15]从“球壳理论”入手，其成果表明数值解无法与“理论”解相互验证，且视电阻率异常位置与异常体实际位置也未能对应，认为探测结果主要为巷道底板而不是巷道前方地质情况的反映，对现有探测方法的有效性提出了质疑。为了进一步完善该技术，王鹏等^[16]通过数值模拟方法计算了点电源在全空间的电位分布，拟合了异常距离与异常特征参数极小值之间的关系式，给直流电阻率法超前探测数据处理提供了新思路，由于文中对不同规则几何体的特征参数选择差异较大，导致该思路在实测数据处理中较难开展。

随着计算机技术的不断发展，基于神经网络算法的电法地球物理响应反演方法越来越受关注。围绕神经网络在直流电阻率法探测领域，多位学者进行了系统研究，分别采用ANN神经网络算法^[17]、免疫遗传算法优化的BP神经网络算法^[18]和自适应神经网络算法^[19]实现了地面直流电阻率探测反演，以及多种群遗传优化的ANN神经网络算法实现井孔直流电阻率法探测反演^[20]，通过与最小二乘法比较，认为人工神经网络算法在反演效率和精度方面均有明显优势。

然而，基于神经网络方法反演矿井掘进工作面直流电阻率法超前探测实测数据尚未进行深入研究。在使用神经网络对直流电阻率法反演时，良好的正演计算方法是实现精准反演的前提。其中有限元法网格剖分灵活，在复杂模型数值计算中优势明显，在高精度直流电阻率法正演中使用较为广泛^[21-23]。

本文基于三维非结构有限元法实现了正演计算，并针对目标地质电性参数建立了矿井直流电阻率法超前探测响应信息库，采用L-M人工神经网络对实测数据进行处理，研究成果对煤矿井下掘进工作面直流电阻率法超前探测方法的完善与推广应用有重要意义，使其更好地服务于矿井防治水工作。

1   矿井直流电阻率法超前探测基本原理

矿井掘进工作面直流电阻率超前探测技术是以含水目标体与围岩的电性差异为物理基础，通过在掘进工作面附近布设点源建立全空间电流场，在已掘巷道中采集电位差异，通过信号处理技术实现含水异常体的超前探测目的。直流电阻率法超前探测装置形式多样，在实际的工程应用中多采用单极−偶极(Pole-Dipole)及其变种组合测量方式，几何示意如图1所示。

图  1  直流电阻率超前探测几何示意

Figure  1.  Geometric of the advanced detection using the DC resistivity method

下载: 全尺寸图片 幻灯片

依据直流电场理论，点电源场的电位公式为：

式中：U为观测点测量电位；I为点源激励电流；$ \rho $为均匀全空间的电阻率；R为观测点与点电源的距离。

由此可得测量电极M和N的电位：

式中：$ {R_{AM}} $、$ {R_{BM}} $、$ {R_{AN}} $和$ {R_{BN}} $分别为对应测量电极与点源之间的距离。由于供电电极B置于无穷远，因此，测量电极M和N的电位差可表示为：

式(3)可改写为：

式中：K为装置系数，式(4)即为均匀介质电阻率的计算公式。

在实际的直流电阻率法探测中，地下介质电性分布不均匀，仍按式(4)计算的结果为不均匀体的综合反映，称之为视电阻率。在直流电阻率法探测中，通常将式(4)中$ \rho $记为测量电极M和N中点处视电阻率。在直流电阻率法超前探测中，M和N电极均位于已掘巷道中，显然依据此方法，现有的超前探测中计算所得视电阻率对应已掘巷道底板区域的地质电性信息更为恰当。

直流电阻率法探测响应是地质区域内介质电性的综合反映，已有文献成果表明直流电阻率法超前探测响应中包含未掘区域地质信息，同时证明超前区域对响应信号的贡献弱于已掘巷道区域^[16]。尽管从“强信号”中识别“弱信号”较为困难，但是在直流电阻率法超前探测信号处理中，通过建立适应的测量响应与未掘区域地质体信息的映射关系，仍有可能实现超前探测目的。

2   基于有限元法的直流电阻率法超前探测数值计算

2.1   直流电阻率法超前探测正演几何模型

在数值建模中，建立包含顶板围岩、煤层和底板围岩的三层地质模型，模型外形是边长为1 000 m的正方体。其中煤层厚度为5 m，在煤层中设置长度为200 m断面为5 m×5 m的空气巷道，电极布设于巷道底板，点源位于距离工作面4 m位置，模型最外层设置为0电势面，接收电极是由距离工作面8~100 m且以4 m等间隔布设的24个电极阵列组成。设置异常体为边长为10 m的正方体且位于掘进区域正前方，其位置可根据数值计算需要在巷道掘进前方移动。

2.2   有限元数值建模

2.2.1   控制方程

尽管在实际的工程应用中，直流电阻率法超前探测使用的是低频交流电，但因为在低频条件下电容和电感效应近似可以忽略，可认为是直流供电。由于观测场的空间范围远大于供电电极本身的大小，采用点电流源作为激励。其控制方程可表示为：

式中：V为电势；$ \sigma $为介质电导率(电阻率的倒数)；$ \delta $为狄拉克函数；R_a为观测点与参考点的距离。

2.2.2   非结构化网格剖分

网格剖分是有限元计算的重要环节，相较于结构化网格剖分，非结构化网格剖分的每个单元之间没有隐含的连通性，剖分方式更加灵活。四面体网格是三维非结构化网格剖分中最常见的一种，理论上可以对任何复杂的几何模型进行剖分，本次采用自由四面体网格剖分方式。依据矿井直流电阻率法超前探测理论，直流电阻率法超前探测需要设置无穷远电极，因此，必须建立足够大的模型才能保证计算精度。为了提高计算效率的同时保证计算精度，对部分区域网格进行加密处理，加密的网格节点主要集中在：① 供电电极与测量电极之间区域；② 低阻异常体与巷道区域；③ 不同介质之间分界面附近区域。

在实际的建模过程中采用控制物理域网格单元尺寸和非结构化网格增长率的方式实现上述网格加密策略。其中，低阻异常体与巷道网格单元大小在0.2~1.0 m之间，最大单元增长率为1.25；煤层区域网格单元大小在0.5~20 m之间，最大单元增长率为1.5；顶、底板围岩区域网格单元大小在5~80 m之间，最大单元增长率为1.8；加密后的直流电阻率法超前探测模型网格剖分如图2所示，共包含2 032 400个单元。

图  2  有限元直流电阻率超前探测模型网格剖分

Figure  2.  Finite element mesh generation for the DC resistivity method-based advance detection model

下载: 全尺寸图片 幻灯片

2.2.3   迭代方程求解

利用非结构化网格对式(5)进行离散化后，得到形式为Ax=b的大型线性方程组，其中，A为包含介质电阻率信息的系数矩阵，b为包含点电源信息的右端项，x为待求解的电场。对该线性方程组的求解而言，稳定双共轭梯度(BiCGSTAB)求解过程比传统的共轭梯度更快速平滑，并且可实现面向非对称系数矩阵的求解^[19,24]。本文采用BiCGSTAB方法对线性方程组Ax=b进行求解，对应的伪代码如图3所示。

图  3  BiCGSTAB算法伪代码

Figure  3.  Pseudocode of the BiCGSTAB algorithm

下载: 全尺寸图片 幻灯片

2.2.4   数值计算有效性验证

为了验证数值计算方法的有效性，建立全空间含球状低阻体模型，采用二极测深装置比较数值解与解析解。模型如图4所示，球状体半径为10 m且位置固定，供电电极A与测量电极M相对于异常体对称分布，A与M的间距从1 m到100 m变化，球心与AM连线的距离为30 m。电性参数中，围岩电阻率为100 Ω·m，球状体电阻率为1 Ω·m。分别采用本文的有限元方法和解析方法进行计算，电流归一化后的接收电势如图4b所示，结果显示数值解与解析解^[25]最大相对误差小于0.2%，表明数值计算方法正确且具有高精度。

图  4  模型验证

Figure  4.  Verification model

下载: 全尺寸图片 幻灯片

3   数据预处理

3.1   数值计算响应与实测响应匹配方法

由于非均质性以及矿井中存在导电介质干扰和探测电极接地不良等因素影响，通常在井下实际采集的响应与数值响应有一定的差异，为了使数值计算响应建立的神经网络算法能够用于实测数据处理，首要解决的是建立数值计算响应与实测响应的匹配方法。

本文采用的数据预处理过程如图5所示，包括获取实际含煤地层电阻率、三维正演响应计算、数值响应与实测响应误差匹配、建立合成响应库。

图  5  数据准备过程主要流程

Figure  5.  Flow chart of data preparation

下载: 全尺寸图片 幻灯片

根据NLH矿区测井资料显示(图6)，3^−1上煤层的电阻率主要为200~320 Ω·m，该煤层厚度约为5 m，顶板范围和底板岩层电阻率约为150 Ω·m。在数值计算中，参考上述地层电阻率变化规律，设置煤层电阻率为300 Ω·m，顶底围岩电阻率均为150 Ω·m，含水异常区的电阻率以1 Ω·m间隔从1 Ω·m增至300 Ω·m，此外设置巷道空间电阻率为10⁷ Ω·m。

图  6  3^-1上煤层附近测井曲线

Figure  6.  Log curves near of coal seam 3^{-1 upper}

下载: 全尺寸图片 幻灯片

经验模态分解(EMD)是新型的序列信号处理方法，在机械振动噪声分离等领域信号处理中取得良好效果^[26]。为了实现直流电阻率法数值响应与实测响应的误差匹配，在巷道底板岩层中布设电极，采用直流电阻率超前探测装置测量了80组样本，标准化后采用EMD方法对实测信号进行模态分解，分解为4个不同的模态，按照矿井电阻率法响应特征，将测量响应归为误差项和趋势项，其中误差项由3个模态组成，如图7所示。

图  7  测量信号模态分解

Figure  7.  Modal decomposition of measured signals

下载: 全尺寸图片 幻灯片

由图7可知，分解得到的趋势项趋于稳定，主要由理想测量响应成分组成，将数值响应与误差项进行响应重构，得到的合成响应与实测响应对比(图8)，合成响应在一定程度上可以与实测响应匹配，本文在构建响应信息库时采用该方法对数值响应进行匹配，以期达到模拟实测数据的目的。

图  8  实测响应与合成响应对比

Figure  8.  Comparison between measured responses and synthetic responses

下载: 全尺寸图片 幻灯片

3.2   神经网络预测模型变量定义

常规的矿井直流电阻率法探测中，通过计算视电阻率推测地层介质相对富水情况。但是在矿井直流电阻率法超前探测中，采集的响应电压与超前探测地层电阻率的数学关系尚未厘清，因此，无法根据测量响应和装置结构直接计算视电阻率，采用视电阻率方式预测含水异常区的条件尚不满足。

本文建立的神经网络超前预测模型中，建立23组归一化电势差与超前水害预测结果的映射关系，并提出用异常体超前距离和异常率实现超前预测。定义异常体的异常率如下：

式中：$ {\rho _{\text{s}}} $为背景地层电阻率；$ {\rho _{\text{w}}} $为含水异常体电阻率。

4   基于人工神经网络直流电阻率法水害超前预测实现

4.1   人工神经网络的基本原理

人工神经网络中最具代表性的是BP神经网络，它是一种多层前馈神经网络，由输入层、隐含层和输出层组成，其中隐含层可由一层或多层组成。BP神经网络的拓扑结构如图9所示^[27-29]。图中X_n是BP神经网络的输入值，Y_m是神经网络的预测值，w_ij和w_jk是网络权值。

图  9  BP神经网络拓扑图

Figure  9.  Topology of BP neural network

下载: 全尺寸图片 幻灯片

标准BP神经网络是一种简单的最速下降静态寻优算法，存在收敛速度慢，容易陷入局部最优解的不足。Levenberg-Marquardt(L-M)算法是一种基于标准数值优化技术的快速算法，它将梯度下降法和高斯−牛顿法相结合，因此，算法本身既具有高斯−牛顿法的局部收敛性，又具有梯度下降法的全局特性^[30]。

设x⁽ⁱ⁾为第i次迭代后的向量，更新后的向量可表示为：

其中，$\Delta {\boldsymbol{x}}$可由牛顿法则求取：

式中：${\nabla ^{\text{2}}}$E(x)为误差指标函数E(x)的Hessian矩阵；$\nabla$E(x)为E(x)的梯度。

设误差指标函数为：

则有：

式中：$e_k({\boldsymbol{x}})$为误差；m为输出变量个数；${{{\boldsymbol{J}}}}({{{\boldsymbol{x}}}}) = {\displaystyle\sum\nolimits_{{1}}^{{m}} {{e_k}({{{\boldsymbol{x}}}})\nabla } ^2}{e_k}({{{\boldsymbol{x}}}})$为Jacobian矩阵，即：

对于高斯−牛顿法则有

L-M算法是高斯−牛顿法的改进，即

式中：$ \mu $为大于0的常数；I为单位矩阵。

式(14)中，当$ \mu $为0时，为高斯−牛顿法；当$ \mu $为很大值时，近似于梯度下降法，因此L-M算法兼顾高斯−牛顿法的收敛性和梯度下降法的全局特性。

4.2   L-M算法神经网络建模

根据前述，将采集到的标准化电势差响应作为神经网络模型的输入信号，将异常体超前距离和异常率作为输出信号，即输入层有23个节点，输出层有2个节点，并采用10个节点的隐含层。从由数值计算得到的3 000组数据中，随机选择2 700组数据作为训练组，250组数据作为验证组，50组数据作为测试组用以测试网络的预测能力。

网络训练性能如图10所示，训练1 000轮次后网络均方误差为0.002 47。使用50组测试样本测试该网络，结果如图11所示，从图中可以看出，对于直流电阻率法超前探测，L-M神经网络对超前探测异常体位置预测的绝对误差小于±0.6 m，对超前探测异常体异常率预测的绝对误差小于±0.05，表明建立的神经网络对于数值响应具有较好的预测能力。

图  10  网络训练性能

Figure  10.  Network performance in training

下载: 全尺寸图片 幻灯片

图  11  L-M神经网络对样本测试结果

Figure  11.  L-M neural network-based sample test results

下载: 全尺寸图片 幻灯片

5   实例应用

5.1   工程概况

NLH矿井位于矿区的最南端，根据以往勘探资料，本区影响3^−1上煤巷道掘进及工作面回采的主要充水含水层为侏罗系中统延安组(J₂y)碎屑岩类承压水层和侏罗系中统直罗组(J₂z)碎屑岩类承压含水层。延安组(J₂y)砂岩含水层为3^−1上煤层的直接充水含水层，在巷道掘进期间对巷道的影响表现为巷道顶板滴淋水；直罗组(J₂z)砂岩含水层为3^−1上煤层的间接充水含水层，富水性较强，若存在断层或裂隙等导水构造，则在巷道掘进期间存在突水的隐患。由于3^−1上煤层巷道掘进工作面前方的含/导水构造尚未查清，巷道前方顶板含水层水对巷道的安全掘进构成一定威胁，因此，需在掘进工作面进行超前探测，以便提前对巷道掘进工作面前方的低阻异常区段采取适当的水文技术措施，保证巷道的掘进安全。

2019年1月16日至12月25日期间，共进行了96次直流电阻率法超前探测，并对巷道掘进过程中揭露的异常区域顶板滴淋水情况做了记录。

5.2   L-M神经网络模型水害超前预测

基于已训练的L-M神经网络模型，对 NLH矿井直流电阻率超前探测进行预测，针对96次实际探测，共得到96组预测异常位置及异常率，预测结果如图12所示。其中，图12a为70组在预测异常位置处揭露有滴淋水时的预测结果，图12b为26组在预测异常位置处揭露无滴淋水时的预测结果，图12c为96组预测结果与实际揭露情况对比。由于施工过程中仅记录了超前100 m范围的滴淋水情况，因此，对于预测距离超过100 m的含水体并未进行验证。

图  12  L-M神经网络预测与揭露滴淋水情况

Figure  12.  L-M Neural Network-based Prediction and revealed results of dripping water

下载: 全尺寸图片 幻灯片

图12a中，70组在预测异常位置处揭露有滴淋水时的预测异常位置在12.1~82.4 m区间，其中65组预测异常体位置在15~80 m之间；图12b中，26组在预测异常位置处揭露无滴淋水时的预测异常位置在9.9~139.3 m之间，其中仅5组预测异常体位置在15~80 m之间。

图12a中，70组在预测异常位置处揭露有滴淋水时的预测异常率值在0.11~0.94之间，且平均值为0.53；图12b中，26组在预测异常位置处揭露无滴淋水时的预测异常率值在0.07~0.43之间，且平均值为0.17，揭露有滴淋水处的预测异常率显著大于无滴淋水处的预测异常率；此外，图12b中，有2个预测异常位置处的预测异常率较大，分别为0.43和0.32，其中0.43对应的预测异常位置为133.2 m，0.32对应的预测异常位置为9.9 m，二者均在15~80 m之外。

综合分析图12中各预测结果与揭露情况，可发现预测结果基本反映了顶板滴淋水的实际情况，全区段预测准确率为67%；当预测异常位置在15~80 m之间时有着更高的准确率，该区间内的预测准确率超过85%，基于L-M算法的神经网络可有效应用于直流电阻率法超前预测。

6   结 论

a. 采用非结构网格剖分和BiCGSTAB迭代有限元方法建立了矿井直流电阻率法超前探测模型，所得数值解与解析解相对误差小于0.2%，该数值计算方法具有较高精度。

b. 针对数值响应与实测数据存在误差干扰情况，采用EMD算法将测量信号分解为趋势项和误差项，并将数值响应与误差项进行响应重构，实现数值计算响应与实测响应的匹配，为神经网络建模奠定数据基础。

c. 通过对数值响应测试发现，基于L-M算法的矿井直流电阻率法超前预测神经网络训练1000轮次后均方误差为0.002 47，对50组测试样本异常体位置预测的绝对误差小于±0.6 m，对异常率预测的绝对误差小于±0.05，该神经网络对于100 m范围内的异常数值响应具有较好的预测能力。

d. 将训练好的神经网络用于实测数据处理发现，该神经网络在100 m范围内预测准确率为67%，当预测异常位置在15~80 m之间时预测准确率超过85%，具有较高的准确率，该神经网络可用于矿井直流电阻率法超前15~80 m范围内异常体预测。

e. 综合分析认为L-M人工神经网络可有效应用于矿井直流电阻率法超前预测，但是其预测结果依赖于建模的响应库。尽管采用了数值响应与实测响应匹配方法，但无法还原井下实际情况。为了能更好地服务于矿井防治水工作，下一步将收集更多的实测数据，开展基于数值响应和实测响应联合构建训练模型，提高预测方法的准确率和普遍适用性。