基于瓦斯实时监测的炮掘工作面爆破自动识别与突出危险性预测

2022-07-20彭玉杰宋大钊李振雷何学秋王洪磊邱黎明

煤炭科学技术 2022年5期

彭玉杰，宋大钊，李振雷，何学秋，王洪磊，邱黎明

(北京科技大学土木与资源工程学院，北京 100083)

0 引言

煤与瓦斯突出是煤矿的主要灾害之一，是瓦斯、地应力、煤的物理力学性质等因素综合作用的结果[1-3]。我国是世界上煤与瓦斯突出最严重的国家之一，突出矿井多、分布广，突出事故多，重特大突出事故时有发生[4]。国内外开采实践表明，煤与瓦斯突出的发生呈区域性分布，突出区域只占整个开采区域的8%～20%，并且多数发生在煤巷掘进工作面[5]。炮掘是煤巷掘进的一种重要方式，爆破作业是一个瞬间剥落破碎煤体，释放游离瓦斯的过程。研究发现，在局部通风机出口风速相对一定及瓦斯传感器悬挂正确的前提下，瓦斯传感器所测得的瓦斯体积分数能反映煤巷炮掘工作面的瓦斯正常解吸能力、炮后扰动解吸能力、煤层瓦斯吸附与游离含量、与采动应力耦合的煤层瓦斯压力、煤层孔裂隙瓦斯补给能力等突出影响因素，炮后瓦斯体积分数异常增大说明该区域影响突出的因素异常增大，即此区域具有较高的突出危险性[6]。因此，采集炮后的瓦斯涌出信息并分析瓦斯涌出特征，可反映工作面前方的突出危险性。

当前普遍采用的突出预测技术，如钻屑指标法、钻孔瓦斯涌出初速度法[7]在突出预测实践中发挥了重要作用，尤其是钻屑瓦斯解吸指标K1值，K1值是煤巷掘进工作面突出危险性预测的一个重要指标，在大多数突出矿井得到了广泛应用，测试结果直接反应工作面前方的突出危险性。K1值的测定往往要经过打钻采样、解吸测定等步骤，单次测量受人为因素影响较大，因此，可将炮后瓦斯涌出特征与K1值结合，以此建立模型，辅助K1值进行突出危险性预测。近年来，国内外学者在炮后瓦斯异常涌出识别、应用瓦斯异常涌出预测突出危险性方面进行了深入研究。朱世松等[8-9]采取分段形态度量方法，建立了瓦斯报警时间序列形态特征库，提出了基于分段形态度量的瓦斯报警信号快速辨识算法，实现了对炮后瓦斯异常涌出的自动识别。杨禹华等[10]通过分析瓦斯含量监测数据和素描图，提取因爆破而引起瓦斯涌出异常的前兆特征，应用模糊模式识别技术对爆破瓦斯涌出异常进行动态预警。魏连江等[11]从K线图角度对瓦斯异常模式进行研究，确立瓦斯异常K线诊断准则，实现了对炮后瓦斯涌出异常模式的自动诊断。CHEN等[12]根据煤层突出强度随释放气体初始膨胀能的增大而增大的规律，提出了煤与瓦斯突出危险性快速预测技术。李希建等[6]在研究煤巷掘进工作面炮后瓦斯涌出规律的基础上，提出了用瓦斯峰谷比值作为预测预报煤巷炮掘工作面煤与瓦斯突出危险性的方法。邹云龙等[13]利用时间序列方法分离了炮掘工作面瓦斯涌出数据的趋势项、周期项以及噪声项，并利用趋势项、周期项以及噪声项之间存在的显著差异对考察巷道的突出危险性进行了预测。这些方法都对瓦斯异常涌出识别与突出预测具有一定的现实意义。

笔者在研究煤巷掘进工作面炮后瓦斯涌出规律的基础上，基于瓦斯异常涌出识别研究现状，首先对瓦斯体积分数监测数据进行处理，提出了爆破事件自动识别的方法；之后分析了炮后瓦斯体积分数增长速率、峰值和衰减速率3个指标与K1值的相关性，并基于这3个指标，建立了炮掘工作面突出危险性预测模型，对工作面的突出危险性进行预测。以此作为现有预测方法的补充，对提高煤与瓦斯突出事故的预测、预防能力有着重要的意义。

1 炮掘工作面爆破事件自动识别

瓦斯涌出异常是重要的煤与瓦斯突出前兆，爆破作业使瓦斯异常涌出，瓦斯涌出异常可视化为瓦斯体积分数曲线中的异常波动。因此，若要自动分析炮后瓦斯体积分数变化特征，挖掘煤与瓦斯突出前兆，首先要从瓦斯体积分数监测曲线中准确识别爆破事件。爆破登记表记录了爆破时间、瓦斯体积分数等信息，但爆破登记受人为因素的影响较大，且爆破登记时间与实际爆破时间存在一定误差。因此，仅靠爆破登记表等现有资料不能精确反映爆破信息。笔者提出一种基于瓦斯体积分数实时监测数据，对爆破事件进行识别提取的方法，可实现爆破事件的自动识别与信息提取。

1.1 爆破事件自动识别方法

1.1.1 自动识别步骤

瓦斯监测系统可实现对掘进工作面瓦斯体积分数的实时监测，瓦斯体积分数曲线的变化能够反映井下通风、异常地质构造、人为作业等环境的变化，同时，不同爆破事件的炮后瓦斯体积分数变化表现出相似性。因此，基于瓦斯体积分数监测数据，提出了一种爆破事件自动识别的方法。爆破事件自动识别的具体步骤如图1所示。

图1 爆破事件自动识别步骤Fig.1 Steps of auto-identification of blasting

1.1.2 自动识别算法

1)数据预处理。

①差分运算。瓦斯体积分数数据是非平稳的离散时间序列，对瓦斯体积分数数据进行差分运算，可以消除一些数据波动，使瓦斯体积分数变化特征规律更加明显。去除因仪器校检等原因产生的异常数据后，取瓦斯体积分数向量为

C=(c(1),c(2),…,c(N))T

(1)

式中：c(t)为t时刻的瓦斯体积分数；t=1,2,…,N，N为瓦斯体积分数监测时间。

取向量C的一阶差分向量

C′=(c′(1),c′(2),…,c′(N-1))T

(2)

c′(t)=c(t+1)-c(t)

(3)

式中：c′(t)为t时刻瓦斯体积分数一阶差分值，t=1,2,…,N。

②卷积运算。将式(2)与向量

V=(5,4,3,2,1,0,-1,-2,-3,-4,-5)T

(4)

进行卷积运算得到向量

W=(w1,w2,…,wN+9)T

(5)

(6)

从向量W中截取与向量C相同维数列向量为

S=(s1,s2,…,sN)T

(7)

si=wi+4,i=1,2,…,N

(8)

向量C与向量S的典型对应关系如图2所示。

图2 差分卷积向量S与瓦斯体积分数向量CFig.2 Difference convolution vector S and gas volume fraction vector C

由图2可知，S中较为明显的波峰波谷与瓦斯体积分数向量C中的爆破与炮后瓦斯体积分数峰值位置表现出良好的对应性。因此，要识别瓦斯体积分数监测曲线中的爆破事件及其特征，需对S中符合一定规则的波峰波谷进行精准有效的识别。

2)识别方法。应用MATLAB中findpeaks函数识别S中的波峰位置，得到对应的爆破时刻及此时的瓦斯体积分数；取S的相反数-S，同样应用findpeaks函数识别-S中的波峰位置，得到对应的炮后瓦斯体积分数峰值及其时刻。findpeaks函数参数设置如下，最小峰高为0.2，最小峰距为90.0，最小峰突为0.3，最小峰宽为3.0。

1.1.3 识别结果评价

1)分类评价。从连续的瓦斯体积分数监测曲线中，识别出爆破事件，可以看作一个二分类问题，即把异常波动的曲线段分为爆破与非爆破2类。为了评价识别方法的有效性与准确性，应用模型评价指标召回率Re与精准率Pr对提出方法进行评价。

召回率Re表示真正爆破的样本中，爆破事件被识别的比例。此时有2种可能，一种是把真正的爆破事件识别为爆破，另一种就是把真正的爆破事件识别为非爆破事件。召回率Re衡量的是识别方法的查全率。召回率Re[14-15]定义为：

(9)

式中：TP为把真正的爆破事件识别为爆破的样本数；FN为把真正的爆破事件识别为非爆破事件的样本数。

精准率Pr表示在识别为爆破的样本中，真正爆破事件的比例。识别为爆破有2种可能，一种是把真正的爆破事件识别为爆破，另一种是把非爆破事件识别为爆破。精准率Pr衡量的是识别方法的查准率。精准率Pr[16-17]定义为：

(10)

式中：FP为把非爆破事件识别为爆破的样本数。

2)数值评价。对识别出的爆破事件，进行爆破时刻、炮后瓦斯体积分数峰值等特征提取，提取准确性的评价指标为平均绝对误差M：

(11)

1.2 爆破事件自动识别应用实例

笔者以某矿炮掘巷道1138和11227运输巷2019年1—4月的瓦斯体积分数监测数据为分析对象。此矿井是煤与瓦斯突出矿井，1138运输巷所属的3号煤层及11227运输巷所属的22号煤层的平均瓦斯含量分别为17.55、10.79 m3/t，平均瓦斯压力分别为0.97、1.80 MPa，瓦斯含量高，瓦斯压力大。掘进过程中易引起局部瓦斯涌出量过大，煤与瓦斯突出危险性高。在掘进巷道需布置瓦斯传感器T1、T2、T3实时监测掘进过程的瓦斯体积分数，T1、T2、T3的布置如图3所示。

图3 掘进工作面瓦斯传感器布置Fig.3 Placement of gas sensors in driving face

由图3可知，T1距掘进工作面最近，可实时反映掘进工作面的瓦斯体积分数变化情况，因此爆破事件的识别采用瓦斯传感器T1的监测数据。应用提出的识别方法对爆破事件进行识别，典型结果如图4所示。图4为从1138运输巷4月9—10日瓦斯体积分数监测曲线中识别出的2次爆破事件，识别出的爆破时刻分别为2019年4月9日15:30和22:35，爆破时刻的瓦斯体积分数分别为0.05%、0.04%；炮后瓦斯体积分数达到峰值的时刻分别为2019年4月9日15:32和22:38，瓦斯体积分数峰值分别为0.34%、0.4%。爆破登记表中的爆破登记时间分别为2019年4月9日15:15和22:45。此方法是基于瓦斯体积分数实时监测数据对爆破时间进行识别提取的，与人工登记的爆破时间相比，提取的爆破时间与实际爆破时间更为接近，此外，该方法还可识别提取出未曾记录的炮后瓦斯体积分数峰值及时刻。将所有识别结果与实际爆破时间和炮后瓦斯体积分数峰值进行比较，应用评价指标对识别结果进行评价。评价结果见表1。

图4 部分爆破事件识别结果Fig.4 Part of auto-identification results of blasting

由表1可知，此识别方法对瓦斯体积分数曲线中爆破事件的识别召回率平均为84.13%、精准率平均为77.21%。将人工标定的实际爆破时间与自动识别出的爆破时间，实际炮后瓦斯体积分数峰值与识别提取出的瓦斯体积分数峰值对比计算其平均绝对误差。炮后瓦斯体积分数峰值平均绝对误差平均为0.018 9%；爆破时刻平均绝对误差平均为2.323 7 min。因此，此方法对于爆破时刻及炮后瓦斯体积分数峰值的识别提取较为精确。

表1 识别结果与方法评价Table 1 Auto-identification results and method evaluation

2 炮掘工作面突出危险性预测

2.1 炮后瓦斯体积分数特征及与K1值的相关性规律

目前普遍采用的煤与瓦斯突出预测指标K1值能够有效地反映突出的危险程度。K1值是煤层每克煤样采集后在第1分钟内的瓦斯解吸量，K1值越大，煤与瓦斯突出危险程度越高[16]。炮后瓦斯体积分数随K1值表现出不同的变化特征，将炮后瓦斯体积分数变化特征与K1值结合，对比不同K1值条件下，炮后瓦斯体积分数曲线拟合特征。

应用前文1138和11227运输巷爆破事件自动识别提取的结果，分析炮后瓦斯体积分数变化特征。可以发现，炮后瓦斯体积分数的变化可以分为2个阶段：第1阶段为直线上升阶段，在此阶段，瓦斯体积分数呈现直线上升的趋势，瓦斯涌出主要受落煤瓦斯涌出及煤壁瓦斯涌出的影响；第2阶段为缓慢下降阶段，此阶段的瓦斯体积分数呈现缓慢下降的趋势，瓦斯涌出主要受煤壁瓦斯涌出的影响。

为分析对比不同K1值条件下，炮后瓦斯体积分数变化的具体特征。对识别出的每一爆破事件，以提取的爆破时刻为起始点，第1阶段每隔10 s提取一个瓦斯体积分数值，第2阶段每隔5 min提取一个瓦斯体积分数值。为了降低偶然误差，对K1值相同的爆破事件，取炮后相同间隔时间瓦斯体积分数的平均值。得到炮后2个阶段不同K1值条件下，炮后瓦斯体积分数均值随时间Δt的变化规律，如图5、图6所示。由图5可知，在炮后瓦斯体积分数变化的第1阶段，即炮后2 min内，瓦斯体积分数变化呈现相同的直线上升变化特征，且K1值越大，炮后瓦斯体积分数增长速率和峰值均越大；由图6可知，在炮后瓦斯体积分数变化的第2阶段，瓦斯体积分数随时间表现出统一的对数衰减特征，且K1值越大，炮后瓦斯体积分数衰减的速率越快，恢复至平稳值所需时间越长。为定量分析炮后瓦斯体积分数增长速率、瓦斯体积分数峰值和瓦斯体积分数衰减速率与K1值的相关关系，对其进行相关性分析。设2个长度为n的序列x(x1,x2,…,xn)、y(y1,y2,…,yn)，两者的相关系数rxy可以用下式表示。

图5 不同K1值条件下炮后瓦斯体积分数第1阶段变化特征Fig.5 Variation characteristics of gas volume fraction in the first stage after blasting under different K1 values

图6 不同K1值条件下炮后瓦斯体积分数第2阶段变化特征Fig.6 Variation characteristics of gas volume fraction in second stage after blasting under different K1 values

(12)

根据式(12)对炮后瓦斯体积分数增长速率、瓦斯体积分数峰值和瓦斯体积分数衰减速率与K1值两两组合求取相关系数，结果如图7所示。可知，炮后瓦斯体积分数增长速率、瓦斯体积分数峰值和瓦斯体积分数衰减速率与K1值的相关系数分别为0.85、0.92、0.79，都与K1值强相关；同时，炮后瓦斯体积分数增长速率、瓦斯体积分数峰值和瓦斯体积分数衰减速率之间也具有很强的相关性。

图7 炮后瓦斯增长速率、峰值、衰减速率和K1值相关性矩阵图Fig.7 Correlation matrix diagram about growth rate,peak value,decay rate of gas after blasing and K1 value

2.2 突出危险性预测模型

炮后瓦斯体积分数增长速率、瓦斯体积分数峰值和瓦斯浓度衰减速率与K1值具有很强的相关性，K1值通常采用钻孔法测定，一次测定往往受人为因素影响较大。因此，利用炮后瓦斯体积分数特征与K1值的强相关性，建立K1值与炮后瓦斯体积分数增长速率、峰值或衰减速率的关系方程，通过每次爆破作业炮后的瓦斯体积分数变化特征，对前方工作面的有突出危险性进行预测。为确定最优回归模型的指标数量，将K1值与炮后瓦斯体积分数增长速率、峰值和衰减速率3个指标组合，建立7个回归模型，以拟合优度R2和残差平方和SE反映模型优劣，拟合优度R2越大，残差平方和SE越小，建立的模型越优。回归模型及参数见表2，拟合优度R2、残差平方和SE的计算公式如下。

表2 回归模型及参数Table 2 Regression model and parameters

(13)

(14)

由表3对比可知，应用炮后瓦斯体积分数增长速率、瓦斯体积分数峰值和瓦斯体积分数衰减速率3个指标建立的回归模型拟合优度R2最大，残差平方和SE最小。由此选择回归模型

K1=0.016r+0.488p+0.004d-0.022

(15)

对掘进工作面前方的突出危险性进行预测。《防治煤与瓦斯突出细则》中规定[17]：反映干煤样突出危险性的指标K1值，其临界值为0.5 cm3/(g·min1/2)，即当测定的K1值超过0.5 cm3/(g·min1/2) 时，预测前方工作面具有突出危险性。因此，反映突出危险性的指标K1值的临界值设为0.5 cm3/(g·min1/2)。

1.3 突出危险性预测实践

1.3.1 11224运输巷及掘进概况

将突出危险性预测模型应用于该矿11224运输巷，11224运输巷所属的22号煤层平均厚度1.2 m，平均原煤瓦斯含量10.8 m3/t，煤层坚固性系数f平均为0.5～0.6。11224运输巷沿22号煤层走向掘进，支护形式为锚网索支护，运输巷设计施工总长度400.18 m，实际施工长度451.5 m。11224运输巷的位置如图8所示。

图8 11224运输巷位置Fig.8 Location of No.11224 headentry

1.3.2 11224运输巷突出危险性预测

喷孔是在突出煤层尤其是松软的高地应力、高瓦斯突出煤层的钻孔施工过程中，从钻孔中喷出煤和瓦斯的现象。喷孔可认为是小型的煤与瓦斯突出[18-19]。首先应用提出的识别方法对11224运输巷2017年8—12月的瓦斯体积分数监测数据进行爆破事件自动识别，之后根据式(15)对识别出的每一爆破事件进行K1预测值的计算，成功地预测出11224工作面前方的25次喷孔和1次小型压出事件。图9选取了11224运输巷2017年8月和12月的瓦斯体积分数监测曲线及对应事件。如图8a所示，在8月9日八点班和零点班的2次爆破事件中，K1预测值分别为0.528 90、0.581 49，均超过临界值0.5，8月10日四点班即发生了喷孔事件，而8月8日、9日K1实测值均为0.27，并未超限；在8月13日4点班的一次爆破事件中，K1预测值为0.662 71，超过临界值0.5，8月15日即实测到K1值为0.52，并且在打钻过程中发生了喷孔。如图8b所示，12月8日，在11224掘进工作面进行了3次爆破作业，第3次爆破事件K1预测值为0.701 12，远超过临界值0.5；12月9—10日即发生了3次喷孔事件，且10日K1实测值为0.52，并伴有1次小型压出。

图9 11224运输巷2017年喷孔、压出事件Fig.9 Jet hole and small outburst of No.11224 headentry in 2017

2017年8—12月，在11224运输巷中共发生32次喷孔事件，应用建立的预测模型，K1预测值超限17次。设定预测周期为2 d，即如果K1预测值超限后的2 d内，有喷孔或其他煤与瓦斯突出现象发生，则认为此次预测成功。在K1预测值的17次超限中，预测成功14次，虚报3次；在发生的32次喷孔事件中，报准25次，漏报7次。为客观反映预测模型的效能，综合考虑报准率、虚报和漏报的影响，采用R评分法对模型效能进行评价，R越大，则效果越好。预报效能R[20]定义为

(16)