冲击地压“双驱动”智能预警架构与工程应用

陈 结,杜俊生,蒲源源 ，姜德义,齐庆新

(1.重庆大学 煤矿灾害动力学与控制国家重点实验室，重庆 400044; 2. 重庆大学 资源与安全学院，重庆 400044; 3. 煤炭科学研究总院 深部开采与冲击地压研究院，北京 100013)

煤炭是我国的主要能源，即使在“双碳”战略的大背景下，直到2030年煤炭在我国一次能源结构中占比仍将达到46%。浅部煤炭资源日益枯竭使我国煤矿开采逐渐转入深部。目前超过一半的煤炭资源储存在1 000 m以深的地层中。为保证持续发展所需的能源供给，势必加快向深部开采的步伐。深部开采过程中，煤岩体应力环境比浅部更加恶劣，具有冲击危险性的矿井数量增多，严重威胁煤矿安全生产。近年来，我国主要产煤基地均有严重的冲击地压事故发生，造成了巨大的人员财产损失。因此，冲击地压的防灾减灾问题研究仍是煤炭行业刚性需求及科研热点。

冲击地压预测预警是防灾减灾的重要环节，可以为冲击地压现场治理提供指导，其核心目标是在灾害发生之前准确判定灾害的时间、区域与危险性等级。目前绝大多数冲击地压预测方法是对采动区域的危险等级进行评估。冲击地压危险性预测初期，主要是通过采动区域状态评价指标参数对区域危险等级进行静态评价。目前主要方法为综合指数法、可能性指标法、地质动力区划法、应力集中评价法及多因素耦合法等。通过液压支架压力监测、巷道锚杆监测、钻屑法监测及煤岩体应力监测对冲击地压危险性进行动态分析。随后，微震、地音、电磁辐射、声发射等地球物理方法广泛用于冲击地压危险性动态评价。相比于静态评价法，地球物理法可以实时捕捉煤岩体破裂、应力集中等现场信息，具有更好的现场工程应用价值。由于冲击地压灾害的致灾因素具有多元性和不确定性，通常将多参量指标结合起来进行更为全面的监测预警。例如，窦林名等综合采场应力监测与微震监测，形成了“震动场-应力状态”耦合一体化监测预警系统，利用获得的波速(梯度)异常系数与煤体应力增量对冲击危险性进行实时预警。姜福兴等以巷道掘进期间的微震、地音和应力信息以及获取的真实地质资料构建了回采工作面冲击危险性动态评估方法，对区域危险性进行评价。刘少虹等建立了一种基于电磁波和地震波CT联合探测冲击危险性的评价方法，实现了掘进工作面及回采工作面大面积冲击危险性等级评价。目前冲击地压危险性的动、静态这种物理模型的评价方法已有相关应用，但通过动-静态协同的多参量物理评价方法进行冲击危险性联合监测预警的工程应用并不常见。

冲击地压预测预警更为关键的问题是灾害时间与区域的量化预测。目前存在的主要问题是缺乏现场监测参量与灾害事件、区域之间的映射关系，导致时间与区域的预测很难量化。学者对于冲击地压时间与区域的研究主要集中在灾害发生以后，对孕灾过程中能量、应力、微震等参量进行回溯分析，研究各参量的时空演化规律，其中绝大部分属于定性或者规律性描述。例如，李康等研究了临空巷道冲击地压耦合致灾因素的时空演化过程，提出了巷道静载应力、动载应力在孕灾过程中的时空演化规律。赵志鹏等研究了近直立煤层冲击地压发生前微震事件的时空演化规律，发现微震频次与能量在冲击地压发生前的异常变化，定性描述了微震能量偏差值与灾害发生的关系。张月征等利用钻孔应变数据的分形盒维数变化表征了冲击地压时空演化过程，发现冲击地压发生前，相关区域应变的分形盒维数会出现“升高→降低”变化。HE等研究了某急倾斜工作面发生冲击地压前的微震事件空间分布，揭示了孕灾过程中灾害状态量的变化情况。然而，上述研究很难直接使用孕灾过程中的现场参量前溯预测灾害的时间与区域，即无法构建监测参数与灾害时间区域的定量映射模型，实现真正意义上的灾害预测。

综上所述，传统的数学物理模型很难量化监测参数与灾害时间区域的定量映射关系，主要原因包括：① 模型泛化性差。冲击地压的数学物理预测模型建立在灾害机理(假说)上，而目前尚无统一的冲击地压灾害机理。基于机理的数学物理模型仅适用于特定条件下的冲击地压预测；② 模型鲁棒性差。数学物理预测模型通常考虑有限几个影响因素，模型中参数所占权重大，一旦有参数获取不精确，将直接影响到模型的预测效果；③ 预测动态性差。采掘过程中冲击地压主控因素会动态变化，而显式的数学物理预测模型很难实时反映参数变化，模型的预测动态性无法保证。但是，以深度学习、知识图谱等为代表的数据驱动方法可以很好地克服数学物理模型的以上弱点，其不关注先验机理知识，不关注参数间的耦合关系，模型可以动态更新，非常适用于冲击地压时间与区域的定量预测。近2年来，已有不少学者将机器学习应用于冲击地压的研究，包括冲击危险性评价、冲击地压微震信号处理等方面。然而，目前还没有直接利用数据驱动的方法进行冲击地压时间区域定量预测的报道。

鉴于此，笔者提出一种全新的冲击地压“双驱动”(物理驱动及数据驱动)智能预警架构：在物理驱动框架下，建立基于贝叶斯概率的动静态协同物理评价模型，对冲击地压危险性等级进行实时评价；在数据驱动框架下，利用深度学习模型实现未来微震事件时间置的精准预测，继而量化冲击地压的区域与时间。基于“双驱动”架构，开发相应的工作面冲击地压智能预警平台，实现预测结果的三维动态可视化，更加直观地展示冲击地压孕灾期间危险特性信息，方便工程技术人员准确、快速制订防灾减灾措施。

1 冲击地压“双驱动”智能预警技术架构

“双驱动”智能预警架构基于集成综合预警思想，利用煤岩体特性、开采参数、实时监测数据及历史灾害资料建立了冲击地压危险等级评价的物理驱动模型和灾害时间、区域预测的数据驱动模型。其预警的整体流程如图1所示。

图1 冲击地压“双驱动”智能预警构架

在物理驱动方面，利用动静态协同算法，采用采区现场勘测得到煤岩体结构特性及采场空间特征信息及煤岩体力学参数，计算冲击地压静态综合指数和可能性指数，形成冲击地压的静态评价体系；同时利用煤矿井下安全监测系统获取的应力、微震参数及反演的波速异常系数，形成动态评价体系。然后利用贝叶斯概率模型对动静态各指标进行危险等级的后验概率计算。最后，根据冲击地压设置的等级阈值确定风险等级。

在数据驱动方面，利用现场微震在线监测系统获取的微震事件属性参数(时间、能量、坐标)作为深度学习模型的训练集。充分考虑微震事件的长短程时序依赖规律，建立具有线性与非线性多模块耦合的深度学习模型MSNet。将工作面连续的微震事件时间序列作为模型输入，预测未来微震事件发生的时间及位置，继而对冲击地压潜在的时间与区域进行定量预测。最终，通过物理和数据“双驱动”模型实现冲击地压综合智能预警。

最后，基于提出的冲击地压“双驱动”预警架构，开发了相应的智能预警平台，实现冲击危险性等级、冲击危险性区域及潜在灾害时间的实时预警和三维动态展示。预警平台主要分为5个模块：现场数据获取与传输、数据存储与前处理、物理与数据驱动模型、实时预警与定量预测、三维显示及警情发布。

图2展示了双驱动模型的冲击地压灾害预警平台的总体框架。从以下4个层面进行简要描述：① 在现场数据获取层面，获取井下采场区域煤岩体赋存特征信息；利用物探及钻探方法获取采动区域及围岩体周围的相关结构特征；通过在线监测系统监测钻孔应力、微震等实时数据，结合煤矿地质资料及历史数据，形成完整的多源异构数据体系；② 在数据存储与前处理层面，利用物联网实时传输技术，将井下在线监测系统获取的数据、灾害区域的静态性指标参数、地质参数及历史数据等实时传输并汇集至指定的数据仓库，通过数据去噪、删除、清洗等形成可用数据；在物理与数据驱动模型层面，分别建立冲击地压危险性等级动静态协同评价及微震事件定量预测模型；③ 在实时预警层面，根据“双驱动”架构的输出结果，确定冲击地压的等级、区域及时间；④ 在三维显示及警情发布层面，借助可视化平台对井下采区三维构造、监测数据及预测预警结果进行平台显示。当预警判别结果表示存在警情时，通过预警信息发布系统对结果进行多终端实时发布，并采取对应的应急避险预案和防控措施。

图2 冲击地压综合预警平台总体框架

2 冲击地压危险性动静态协同评价模型

2.1 静态评价

综合指数法是《防治煤矿冲击地压细则》推荐采用的方法。本研究从地质因素和开采技术两方面进行综合指数评价计算，并从应力和煤层冲击倾向性进行可能性指标计算，基于综合指数评价计算和可能性指标计算获得区域静态评价结果。

2.1.1 综合指数评价计算

从地质因素和开采扰动2个方面的综合指数来讨论冲击危险性静态指标，地质因素的危险性指数的计算公式如下：

(1)

式中，为第个地质影响因素的实际指数；max为第个地质影响因素中最大的危险指数；为地质影响因素的数量。

同理，根据细则对冲击地压开采扰动因素危险性指数进行分析,计算公式如下：

(2)

式中，为第个开采扰动影响因素的实际指数；max为第个开采扰动影响因素中最大的危险指数；为开采扰动影响因素的数量。

最后，综合地质因素及开采扰动因素指标可描述为

=max(,)

(3)

其中，为冲击地压综合指数的分数值。冲击地压区域危险程度则通过进行分类，具体的等级划分方法见表1。需要注意的是，本文后续的钻孔应力、微震参数各指标及综合评价的风险评级按照表1的划分标准进行。

表1 静态指标综合指数评价标准

2.1.2 可能性指标计算

一般可能性指数分为2部分，第1部分为应力对冲击地压的贡献，第2部分为煤层本身冲击倾向性对冲击地压的贡献。

应力对冲击地压事件的贡献度为

(4)

式中，=[]，为应力集中系数，为岩石容重，N/m，为开采深度，m，[]为煤层单轴抗压强度，MPa。

煤层冲击倾向性对发生冲击地压的贡献度W可表述为

(5)

式中，为冲击倾向性的弹性能量指数。

冲击地压的可能性综合指数为

(6)

冲击地压可能性指标的等级评价可按照表2的标准进行等级划分。

表2 冲击地压可能性指数评价标准

2.2 基于动态应力监测的冲击危险性评价

2.2.1 钻孔应力指标

为掌握采场区域的应力状态，对工作面上下巷钻孔应力的数据进行分析。基于时间序列的某测点所反应的围岩应力等级状态通过式(7)表示。

=++

(7)

式中，为钻孔应力综合表征指标；为钻孔应力值的危险指标；为应力增幅值的危险指标，即当前测点应力值与最初设点应力值之差；为应力增速值的危险指标；，，分别为，，的权重比，根据现场应力特征及专家意见，本文分别取0.5，0.3，0.2。

由式(8)计算得出。

(8)

式中，为测点的应力，MPa。

由式(9)计算得出。

(9)

式中，为测点的当前应力与开始统计的应力之差，MPa。

由式(10)计算得出。

(10)

式中，为测点当前应力与上一刻的应力差。

若式(8)～(10)的，，取值不大于0时，，，取值为0。

2.2.2 波速异常系数

为了探测煤层区域应力场的相对大小，利用地震CT对煤层监测区域的波速场进行反演，并基于反演结果计算波速异常系数()。

(11)

2.3 微震多维信息动态评价

(1)微震值。

根据古登堡和里克特研究得到的地震频度和地震震级之间的幂率关系，学者们通过变换将地震震级和频度联系起来，建立如下方程：

lg(≥)=-

(12)

其中，和均为常数，为统计区域内地震活动水平的表征参数，为该区域内大小震级数量之间的一种比例关系参数，大震级事件相对较多时，则变小，可以一定程度上表征该区域的应力变化和煤岩体的破裂情况。

为研究微震震级与频度的关系，通过最小二乘法拟合计算可得到，，其中

(13)

式中，为第档的微震震级，在整个时间窗内，对每个事件能级进行大小排序，以最小的微震事件能级作为起始震级，震级档间距Δ设为0.2；为第档震级内的微震实际事件数；为事件窗口内震级的总档数。

(2)微震()值。

微震的频次、震级与微震活动具有正相关性。微震频次和震级越大，冲击危险性也越大。研究表明，()是一个能够反映某一区域微震活动程度的定量参数，可以较为直观反映微震活动的增强或平静趋势。()定义如下：

(14)

式中，为第个微震事件的能量；lg为第个微震事件的能级；为时间窗内微震事件总次数。

(3)累计差量指标。

微震累计差量是通过分析某时间段内微震事件率的累计差量变化来判定采区围岩体的稳定性，设累计差量指标为，大事件累计差量指标为，事件率累计差量指标为，能率累计差量指标为。

根据马雄忠等提出的声发射累计差量指标算法来计算微震事件率累计差量指标，设时间窗内的微震事件率数据序列为,,…,-1,(≥4)，其中，为时序为时的事件率。

对序列项前后相减，得

令=max{Δ,0},=min{Δ,0},

将正差量的和设为，负差量的绝对值和设为，即

(15)

(4)等效能级参数。

某段时间内发生微震事件的能级采用如下计算方式：

(16)

式中，为时间窗内微震事件的总能量；为时间窗内第个事件的能量；为微震事件总数。

等效能级参数()是反映新的微震事件相对于正常水平的偏离程度的指标，具体公式如下：

(17)

2.4 基于贝叶斯方法的动静态协同评价模型

根据前文提出的冲击地压静态评价指标(，)、应力动态评价指标(，)及微震多维信息指标(，()，，)，形成了协同评价冲击地压的动静态综合指标体系。由于在统计时间区间内，各指标的计算值存在未知性。因此，首先通过归一化的方法对各指标进行规范化处理，转换为统一的正向指标。然后利用贝叶斯方法建立动静态各指标参数的危险等级后验概率关系，最后利用常规的加权平均模糊算法对建立的动静态指标进行协同综合评价。

贝叶斯概率的基本思想是：假定所研究的对象具有一定的认识，并且这些对象存在设定的先验概率分布，然后通过新的样本对设定的先验分布进行不断更新，从而得到后验分布概率。本研究利用贝叶斯概率方法进行冲击地压多参量综合危险性判别。贝叶斯概率计算公式为

(18)

式中，()为事件的先验概率；(|)为事件的后验概率；(|)()为调整系数。

将贝叶斯概率公式引入冲击地压评价时，式(18)改写为

(19)

式中，为各评价指标计算值(=1,2,3,…,);为划分的危险等级(=1,2,3,4)；为各指标设置的等级阈值。

则利用贝叶斯概率综合评价冲击地压的计算步骤如下：

(1)计算()，即各评价指标参数值在冲击危险性各等级的先验概率。这里设定各指标计算结果处于每个风险等级的概率相等，即，(1)=(2)=(3)=(4)=14。

(2)计算(|)，概率采用几何概率中距离计算值法，将各评价指标计算值与设定的危险等级评价标准之间距离的绝对值倒数进行计算可得

(20)

其中，=|-|。越小，则表示所测指标属于对应等级的概率越大。

(3)根据式(19)，计算出(|)，即各评价指标参数值在冲击危险性各等级的后验概率。然后通过加权平均模糊综合算法计算动静态各指标在各危险等级下冲击地压综合后验概率，公式如下：

(21)

式中，为各动静态指标的权重。

笔者基于“应力场-震动波场”监测预警方法，选择微震多维信息、在线钻孔应力及地震-CT波速异常系数作为应力场-震动场动态综合指标体系，结合冲击地压静态评价体系，形成冲击地压危险性动静态协同综合评价模型，其综合评级流程如图3所示。

图3 冲击地压危险性等级综合评价流程

3 微震事件定量预测数据驱动模型

3.1 模型基础数据

笔者使用的微震数据来自某煤矿13230工作面。选取在2019-01-01—2019-12-31一年内工作面推进过程监测到的10 196次微震事件作为本文的基础数据。微震事件预测是时间序列预测任务，模型的输入是连续的微震事件序列，其中序列长度为模型超参数，笔者选取序列事件长度为12，即利用12个连续的微震事件预测未来的事件。模型的输出是即将发生的微震事件序列，选择模型要预测的事件序列长度为6。由于监测的微震事件在时间轴上的分布是不均匀的，因此，本研究为每个微震事件引入一个时移参量(单位为min)，表征其与前一个事件的时间差。MSNet模型输入数据为长度为12的微震事件序列，每个事件包含空间位置(,,坐标)、能量、时移5个属性，单个输入样本维度为12×5。MSNet模型输出数据为长度为6的微震事件序列，事件属性与输入数据相同，单个输出样本维度为6×5。表3展示了工作面输入微震事件时间、空间坐标、事件能量、时移的部分数据。

表3 工作面输入微震事件属性

3.2 数据驱动模型设计

微震事件在时间轴上的分布有长短程依赖特征，即一个微震事件或微震事件序列的发生不仅与之前短时间范围内(如几小时内)的微震事件有关，也会受到更长时间范围(数天到数月)的微震事件影响。例如，在冲击地压发生前的一小段时间范围会发生大量能量各异的微震事件，此为短程依赖；同时，顶板来压等长周期事件也会影响微震事件的发生，此为长程依赖。因此，本研究建立可以同时获取短程与长程的微震事件时序规律的MSNet模型，更精确地预测未来发生的微震事件。如图4所示，MSNet的数据流向有两大路径，其中线性路径是利用AR自回归模型，输入的微震事件序列直接预测未来的微震事件(序列)；非线性路径则通过卷积层、循环神经网络层、跳跃循环网络等模块，实现输入的微震事件序列到未来微震事件序列的非线性映射。MSNet的最终预测结果是线性路径结果与非线性结果的加和。下面依次介绍MSNet模型用到的几个关键模块，卷积模块，循环神经网络模块，跳跃循环网络模块、自回归模块。

图4 MSNet模型结构

3.2.1 MSNet非线性模块

MSNet模型的非线性模块主要包括卷积层、循环神经网络及跳跃循环网络3个模块，通过这些模块的综合训练实现微震事件的非线性预测。数据首先流入MSNet的卷积模块，用于解析12个(单个序列长度)微震事件间的短程依赖关系以及事件5个属性间的关系。如图4所示，卷积核尺寸确定为6×5，卷积核数量为32个。卷积层激活函数选择为ReLU(其表达式为ReLU()=max(0,))。式(22)展示了第个卷积核的运算过程，为卷积核矩阵，为输入的时间序列矩阵，为偏置，*为矩阵的哈密顿算子。运算结果采用zero-padding保证输入前后的维度一致。

=ReLU(*+)

(22)

(23)

MSNet设计一个跳跃循环网络(skip-LSTM)，解决潜在的超长程依赖问题。流入LSTM单元的卷积层输出不是连续时间序列，而是每隔一个特定周期(跳跃步长)选取一个数据点，即把当前隐藏单元与相邻周期处于同一位置的隐藏单元链接起来作为LSTM的输入，以扩展信息流的时间跨度。引入跳跃步长后，skip-LSTM的数据流向如式(24)所示。循环神经网络模块及跳跃循环网络模块的输出拼接在一起形成一个长向量，然后通过一个隐藏层的全连接网络得到最终的非线性路径输出。

(24)

式中，-和-分别为前时刻的隐藏状态和记忆单元。

3.2.2 MSNet线性模块

为了充分考虑事件之间的潜在线性关系，MSNet设计了一条平行于非线性路径的自回归线性路径。假设各微震事件属性向量之间具有线性关系，利用前次的微震事件预测第(+1)次事件。本研究中自回归模型如式(25)所示，和为模型参数；为待预测事件的属性向量。MSNet中的自回归模块利用单隐藏层的全连接网络实现。非线性路径结果与线性路径结果进行向量加和，得到最终的模型预测结果。

(25)

3.3 模型训练

深度学习模型通过训练迭代更新模型参数，使损失函数(预测结果与真实值之间的差异)达到最小值。本研究中，模型的预测输出是未来连续6次微震事件的属性(维度为6×5的向量)，因此定义MSNet的损失函数是五维空间中2个长度为6的点序列之间的平均距离，损失函数如下

(26)

10 196次微震事件数据集人为划分为训练集与测试集，前者用于训练模型使之达到最佳拟合效果，后者用于测试模型表现。考虑到数据样本有限，数据集中不再单独划分验证集用于模型超参数优化。选择数据集中2019年1—10月的9 863次微震事件进行模型训练，11月与12月的333次微震事件进行模型测试。为了尽可能多地利用微震数据，模型的输入微震序列之间取步长为1，即第1个训练/测试样本输入为～(+11)的12次连续微震事件，预测(+12)～(+17)的6次微震事件的坐标、能量、时移。本研究在构建训练/测试样本时，没有样本跨越划分时间点(11月1日0时)，保证训练样本和测试样本完全隔开，无数据无泄漏。

3.4 模型结果分析

在整个训练过程中，采用小批量梯度下降法更新模型参数，训练损失与测试损失最后都达到收敛，模型没有出现过拟合现象。把样本输出的连续6个微震事件属性分别看作6个长度为5的向量，利用余弦相似度来衡量属性的预测值与真实值之间的相似性(值越接近于1，表示2个向量越相似，预测精度越高)。其计算方式如式(27)所示，其中为微震事件真实的属性向量；为模型预测的属性向量；为向量长度。

(27)

模型结果表明，模型对微震事件坐标与坐标的预测效果最好，测试样本的平均余弦相似度分别为0.997与0.995。模型对微震事件坐标的预测结果稍差，明显低于与的余弦相似度，其平均相似度为0.939。模型对事件时移的预测出现了明显波动，其平均的预测余弦相似度为0.523，表明其预测值与真实值之间仍存在不错的相关性，MSNet模型挖掘出了微震序列的时间规律。模型对微震事件能量几乎没能做出正确预测，超过一半测试样本的余弦相似度徘徊在0附近，表明预测值与真实值之间几乎没有关联。综合分析，MSNet模型对微震事件的位置预测最准确，时间预测精度次之，对微震事件的能量预测精度较差。

图5展示了模型训练数据(1—10月的9 863次微震事件)的5个属性进行[0,1]标准化之后的分布箱型图。

图5 模型训练数据分布箱型图

图5中绿色菱形表示数据异常值，蓝色圆点表示数据平均值，长方形上下边缘分别代表属性的上四分位数与下四分位数，红色长线表示数据中位数。从图5中可以发现微震事件的位置数据(,,)分布较为理想，中位数与平均值均处于数据集中区域，且与平均值相差不大，异常值并未过多影响数据分布。微震事件的时移数据分布稍差，数据集中在一个较小区域，有部分异常值超过了数据内限，影响了数据分布，但其中位数与平均值仍处于数据集中区域。微震事件能量分布出现了明显的极端异常值，其平均值与中位数处于数据集中区，但是数据能级出现大范围空白，直接影响了数据分布。这是由于工程现场绝大部分微震事件是小能量事件，超过1×10J的大能量事件占比极少。MSNet模型训练输入的绝大部分能量数据位于小能级，因此模型对大能量事件预测结果带有明显的小能量偏向性(bias)，能量数据不平衡(imbalance)极大影响了模型对事件能量的预测能力。

4 基于“双驱动”模型的冲击地压智能预警平台开发

4.1 平台三维几何模型

基于矿井采掘工程平面图、矿井地貌图、地质说明图、早期勘探及测绘相关资料，利用3Dmax建模软件通过虚拟三维空间构建三维几何模型，能够对回采工作面、掘进工作面、巷道布置、围岩体等元素进行展示，准确表征断层、岩性、煤层厚度变化等特征信息。如图6所示，通过Unity 3D平台技术建立了煤矿地质三维模型，包括整个矿区地上地下的三维空间结构、各采区的相对位置、巷道的地质构造及岩性等特征信息，实现了较为准确的三维表征。

图6 采区三维几何模型

4.2 “双驱动”预警模块实现流程

平台搭建预警模块的步骤如下：第1步，数据转换。对监测系统获得的数据进行数据转换，把信息转换成机器能够接收的形式。第2步，数据分组。将数据按照不同类型进行分组，对煤岩数据、地质数据、回采数据、微震数据等进行分类。第3步，数据组织。把分好类的数据放到一个容器中，主要有数组、链表、列表等，笔者选择的是列表容器。第4步，数据计算、存储。通过C#语言对评价模型算法进行代码编写，实现对接入数据的计算，并对数据进行保存。第5步，数据显示。根据需求将容器当中需要进行界面显示的计算结果通过Unity3D平台进行展示。

对于微震事件预测模块，在训练完毕以后进行在线部署，利用实时获取的微震事件序列对未来事件的属性进行动态预测。采用Intel公司的OpenVINO深度学习套件进行模型的在线部署。首先从Keras平台导出TensorFlow冻结图模型文件(*.pb文件)，使用OpenVINO套件中的Model Optimizer工具优化导出的冻结图模型，获得描述MSNet网络拓扑结构的*.xml文件与储存模型参数的*.bin文件，调用套件里面的推理引擎(Inference Engine)完成AI的推理计算。MSNet做预测时输入的实时数据从AraWin微震数据处理软件自带数据库中获取，模型的预测结果以文本形式本地保存并被实时调用。

在数据处理和代码化的一系列工作之后，实现最终的三维预警平台需要进行如下操作：对所用到的应用程序接口、编写的代码、3D模型、3D模型材质以及应用程序需要展示界面的图片进行打包集成，通过unity导出形成一个三维透明化的应用程序。通过数据收集及处理，得到评价可用数据，将数据按动静态类别整理形成动静态指标数据库。预警软件平台调用数据库数据，利用嵌入的评价模型对冲击地压综合危险性进行等级评价并进行区域与时间的预警。图7展示了具有数据处理及评价功能的三维预警平台。三维平台界面右边主要显示微震事件及应力数据，左边显示曲线变化。

图7 三维预警平台危险性等级展示

5 工程应用

现场应用地点为河南省义马市耿村煤矿，其示范工作面煤层厚度为0～20.32 m，平均厚度为9.86 m，倾角8°～15°，平均10°，工作面斜长245 m，走向长为745 m，平均采深为630 m。近水平层理，裂隙发育，厚度10.7～17.8 m，煤层结构复杂，含夹矸3～5层，夹矸为灰黑色泥岩或细砂岩，夹矸累厚1.5 m，为较稳定煤层，煤层变化系数为12%。直接顶为泥岩，厚4.04～42.64 m，深灰色，致密，块状构造，含植物化石。煤层直接底为泥岩，厚度约2.8 m，基本底为泥岩砂岩交替层，厚度35 m以上，致密坚硬。基本底为泥岩、细～中砂岩和砾岩，厚0～35.7 m，平均13.24 m，煤层单向抗压强度8.58 MPa，煤的弹性能指数2.54。

采用井下监测系统和预警系统对工作面冲击地压灾害进行监测预警，如图8所示。在工作面前方250 m左右区域的上、下巷外帮埋设应力计，浅孔埋深12 m，深孔埋深18 m，浅孔与深孔每组应力计的间距为0.5～1.0 m，每组20～25 m。微震拾震器布署在整个工作面上下巷道，间距为250 m左右。随着工作面的推进，钻孔应力计及微震拾震器不断向前移动。通过应力数据、微震及结合地震CT形成的波速异常系数对工作面冲击危险性进行实时动态评价。

图8 工作面应力-微震在线监测设计

(1)应力监测结果。通过地震CT主动爆破的探测方式，获得工作面监测区域相对应力情况，再利用产生的微震事件作为被动震源进行三维波速反演，形成三维速度梯度。根据反演的三维数据进行三维建模，形成三维数据云图。同时结合监测区域钻孔应力的变化特征，实现对工作面采动区域应力集中情况透明化显示，如图9所示。

图9 工作面三维应力云图

(2)微震监测结果。图10展示了2019年监测的微震事件及部分微震多维指标的变化结果，其中图10(c)，(d)显示了2019-08-01—2019-08-07微震多维指标的变化曲线。从微震事件和指标对比来看，微震出现较为明显的低值异常，与微震事件能量具有负相关。()和等效能级参数指标的变化规律具有相似性，同累计差量指标具有明显高值异常。

图10 监测的微震事件及微震多维指标变化曲线

(3)冲击危险性等级预警结果。当平台接入评价指标数据库后，某个指标参数的改变都会使其区域灾害危险性的结果及时更新。评价指标参数一旦发生改变，平台则立刻对指标进行重新计算，实时更新综合等级结果。

通过预警平台对示范工作面2019年的数据进行现场测试，对比现场实际危险情况，确定预警模型的准确率。根据示范工作面的现场防冲日志报告，在2019年，共有25 d具有危险性，均属于弱危险性。根据示范工作面2019年现场实际的动力现象和工人勘测结果及管理统计来看，没有特别明显的支架压力增大、卡钻等中等及强冲击危险性现象，顶底板、支架未见明显变形，时有轻微煤爆或煤炮声等。根据《防治煤矿冲击地压细则》并结合现场专家意见，认为该示范工作面2019年的防冲日志具有参考意义和可比性。

通过平台的预警结果表明，有22 d的预警结果同矿上的评价等级及实际危险情况是相符合的，有3 d出现了不一致，其预警平台的评价结果均显示为中等危险，都高于该矿的防冲日报结果，具体日期见表4。因此，在冲击危险性预警方面，平台的准确率为88%。需要说明的是，该模型准确率验证的前提是防冲日志结果和实际现场情况是吻合的。评价结果产生偏差的原因可能是各指标结果的综合叠加计算引起的，可通过不断优化各指标的权重比和指标的隶属关系来提高准确率。

表4 冲击危险性平台预测结果与现场结果对比

(4)双驱动预警结果。结合物理驱动与数据驱动的预测输出，形成冲击地压等级、时间、区域的综合预警结果。图11为冲击地压“双驱动”综合预警结果，图中黄色和淡蓝色圆圈分别为平台确定的冲击地压危险性等级区域和预测的下一时刻潜在的危险区域，此时的危险性等级为安全。随着采掘的推进，危险性等级区域将动态化地涵盖整个回采工作面。危险预测区域确定的方式是，以预测的6个微震事件中能量最大的事件位置为圆心画圆，其边缘通过距离圆心最远的微震事件位置，且危险性沿径向递减。该平台一直保持显示最近12个连续微震事件，同时预测接下来6个微震事件。当微震系统监测到最新的微震事件时，平台右侧的数据将会动态更新，并同步获得最新的预警结果。冲击地压时间预测的方式是，首先根据本矿井冲击地压历史事故发生前监测到的微震事件最大能量值作为判断阈值，一旦预测的微震事件能量超过该值，则触发冲击地压时间预警。

图11 冲击地压”双驱动”预警三维平台

6 结论与展望

(1)建立冲击地压多源多参量动态评级体系，对扰动区域局部近场危险性进行评价。利用贝叶斯组合概率方法，实现了冲击地压动静态协同的远近场联合监测预警。通过嵌入评价模型的三维预警平台进行测试对比，在一年的防冲日报表里，具有冲击危险性的25 d中有22 d评价结果相同，另外3 d的评价结果偏高，其危险性评价模型的准确率为88%，但仍需进一步扩大示范应用地点和模型的应用时间，以提高模型适用性和精度。

(2)以微震事件定量预测为切入点，量化冲击地压的危险性区域与时间。研究构建了一个普通卷积模块、循环神经网络模块与AR自回归模型相结合的深度学习网络MSNet，利用深度学习技术首次实现了未来微震事件的时间、地点、能量预测。预测结果表明，MSNet模型对微震事件位置的预测最为准确，,,坐标的预测余弦相似度都超过了0.93，对事件发生时间的预测结果稍差，平均余弦相似度为0.523。由于训练样本中能量属性的不平衡，MSNet对事件能量的预测具有明显的小能量偏向性。

(3)开发了冲击地压综合智能预警3D可视化平台，实现了现场监测数据的实时采集、同步传输、智能预警、动态展示，直观地展示了冲击危险性等级，危险性区域及潜在的灾害时间。

该研究的核心是构建贝叶斯动静态协同评价模型进行冲击地压危险性等级预测，同时利用MSNet模型预测未来微震事件，并以此为依据判定的冲击地压危险区域及时间。MSNet的预测精度直接决定了冲击地压量化预测的准确性。MSNet虽然在位置预测方面精度很高，但在能量及时间预测上仍有提升空间。以下列出的几个方面可以成为进一步研究的重点：

(1)调整模型的损失函数结构。该研究计算MSNet模型损失时(式(26))，把独立微震事件作为一个长度为5的向量，并计算真实值与预测值的向量距离。这其中隐含了一个假设，即向量的5个维度(、、、能量、时移)是等权重的。如果工程现场更关注预测能量准确性，可以将向量维度中的时移权重调小，从而让模型损失的大部分由其他4个属性的损失构成，保证能量预测精度。

(2)考虑微震事件的全波形信息。微震事件的全波形是压电传感器电压随时间变化的曲线，其包含了震源位置、介质属性等，信息量丰富。而本研究中对波形特征只考虑了微震事件能量(通过波形的时间积分求得)与时间(通过AraWin软件拾取到时)，忽略了其他大量波形信息。未来研究中，可以把微震事件全波形(独立的超高维向量)作为模型的一个输入属性，这将大幅提升模型信息量，有望提升模型能量与时间的预测精度。

(3)冲击危险性区域与时间确定。本研究中确定危险性区域的方式并没有考虑6个微震事件在能量与时间上的梯度变化，因此确定的冲击区域并无危险性高低的进一步精准划分。另一方面，该研究给出的冲击地压时间的预测方法只考虑了历史灾害数据中能量最大的微震事件作为触发时间预警的阈值。未来研究中可以更进一步建立历史冲击地压发生前的微震事件的能量序列分布，从数据分布中提取特征值或者特征向量作为判断阈值，增强预测的鲁棒性。

动力灾害的智能监测预警是煤矿智能化的重要组成部分。笔者利用物理驱动与数据驱动相结合的方法解决冲击地压量化预测的难题，既考虑了冲击地压发生的物理机制，又充分利用了海量的灾害前兆数据。研究提出的“双驱动”智能预警架构，给出了冲击地压的危险性等级、区域及发生时间的预测方法。同时，研究给出了“双驱动”预警架构未来的研究方向，为彻底解决冲击地压的量化预测难题奠定基础。