韩刚， 解嘉豪， 秦喜文， 王星， 郝晓琦

（1. 中煤能源研究院有限责任公司，陕西 西安 710054；2. 中煤冲击地压与水害防治研究中心，内蒙古 鄂尔多斯 017200；3. 中煤西安设计工程有限责任公司，陕西 西安 710054）

0 引言

冲击地压是我国煤矿开采过程中面临的一种典型动力灾害[1-2]。近年来，我国煤矿多次发生重大冲击地压事故，造成严重的经济损失和人员伤亡[3]。冲击地压危险性评价是冲击地压防治的基础，《防治煤矿冲击地压细则》明确规定：开采具有冲击倾向性的煤层，必须进行冲击地压危险性评价。冲击地压矿井每年需开展大量的冲击地压危险性评价工作，编制报告多、评价周期长、人工成本消耗大，进一步加重了矿井人员的工作负担。在智能化背景下，借助先进的信息技术，实现冲击地压危险区域智能化评价，对于提高煤矿智能化水平、为防冲工作减负具有较大的现实意义。

冲击地压危险区域的准确识别是实现智能化评价的关键。已有国内学者借助大数据、云计算、图像识别等新兴技术，开展了冲击地压危险区域的评价和识别工作。窦林名等[4]、何生全等[5]基于GIS、大数据、物联网、云计算、采矿地球物理等技术，搭建了冲击矿压风险智能判识与多参量监测预警云平台；潘俊锋等[6]开发了基于冲击启动载荷分源监测的冲击地压分源权重综合监测预警平台；姜福兴等[7]运用大数据分析方法和云平台技术，开发了多参量联合监测的冲击地压监控预警平台。上述平台侧重于冲击地压危险区域的动态识别和监测预警，但尚未集成冲击地压危险区域的采前评价和静态识别功能。现有的冲击地压危险性评价方法主要有多因素耦合叠加法、地质动力区划法、模糊层次分析评价法、基于应力与围岩分类的冲击地压危险性评价方法等[8-13]。相比冲击地压危险区域的动态识别和监测预警技术，现有评价方法还存在手动计算量大、危险区域划分精度低、智能化水平不高等问题。

为进一步提高煤矿冲击地压危险性评价的智能化水平，借助现有计算机视觉技术，提出了一种基于图像识别技术的冲击地压危险区域智能化评价方法。基于OpenCV 机器视觉库和深度学习模型，实现冲击地压危险主控因素的精确识别；通过构建应力集中系数与图像灰阶的映射矩阵，实现主控因素影响范围和应力增量的叠加计算；构建了冲击地压危险区域的“无、弱、中等、强”4 级判别标准，实现冲击地压危险的分级分区评价。基于Python 程序语言，开发了冲击地压危险智能化评价软件，为煤矿采掘区域的冲击地压危险性评价提供了新的方法和工具。

1 冲击地压危险性评价的理论基础

1.1 冲击地压危险主控因素及影响范围

研究普遍认为，冲击地压的发生主要受力源因素、内在因素和结构因素3 个方面的影响[14]。其中内在因素为煤岩体固有的冲击倾向性，结构因素为能量释放的软弱结构面，但凡有采掘活动，该因素就难以避免。因此，精确识别各项力源因素是进行冲击地压危险性评价和冲击地压防治的重点。窦林名等[15]将冲击地压发生的力源因素分为静载主控因素和动载主控因素两类，认为当煤岩体中的静载与动载相叠加超过冲击地压发生的临界应力时，冲击地压显现。因此，通过识别主控因素，并对各项主控因素的影响范围和应力增量进行叠加，得到采掘区域的总应力估算云图，以此作为评价区域冲击地压危险性的依据。

在国内主要冲击地压矿井中，造成采掘工作面静载应力集中的因素主要有自重应力、构造应力和采动应力，其中，自重应力的主控因素为采深；构造应力的主控因素为褶曲、断层等地质构造；采动应力的主控因素为采空区、区段煤柱、开切眼外错、终采线外错、巷道交叉等。造成动载应力集中的主控因素有工作面初次来压、“见方”破断及矿震动载等因素。11 项冲击地压危险主控因素中，煤层埋深、矿震动载为全局影响因素，其余因素均为局部影响因素。冲击地压危险主控因素及影响范围如图1 所示。

图1 冲击地压危险主控因素及影响范围Fig. 1 Main controlling factors and influence range of rock burst hazard

1.2 主控因素的量化表征与图像识别

图像识别技术是人工智能的一个重要领域。结合OpenCV 机器视觉库和深度学习模型的图像识别技术对于信息的处理、加工、存储速度更快，在处理细节时更加精细，可以依据图像识别需求对海量信息进行分类提取，大幅提高信息处理的效率[16-18]。冲击地压危险的局部影响因素中，断层、褶曲、采空区、区段煤柱、开切眼外错、终采线外错、巷道交叉7 项影响因素具有明确的识别主体，利用OpenCV 机器视觉库，通过对识别主体进行深度学习训练，即可准确提取识别主体；初次来压、“见方”破断2 项影响因素则需要在识别开切眼位置的基础上，通过输入来压步距、工作面长度等关键参数进行信息识别。

识别主体在不同影响范围下，需分别赋值不同的应力集中系数。各项主控因素的影响范围、应力增量在参考行业经验的基础上，主要采用数值模拟方式进行半定量化估算确定[19]。以断层、采空区、区段煤柱、终采线外错、巷道交叉等局部影响因素的量化表征方法为例进行具体说明。

1.2.1 断层因素

根据某矿31103 工作面3DHF1 断层实际条件（倾角为45°，落差为25 m），建立300 m×250 m×165 m（长×宽×高）的FLAC3D 数值模型，模型采用Mohr-Coulomb 本构，结构如图2 所示。

图2 含断层面的FLAC3D 数值模型Fig. 2 FLAC3D numerical model with fault planes

开挖数值模型，模拟得到工作面回采接近断层过程中工作面前方及断层面的煤体应力集中系数变化曲线，如图3 所示。可看出当回采工作面与断层之间距离大于50 m 时，煤体应力集中系数小于1.2；当距离为30～50 m 时，应力集中系数为1.2～1.4；当距离为10～30 m 时，应力集中系数为1.4～1.9；当距离小于10 m 时，应力集中系数可达2.05。据此划分断层前后30～50 m、10～30 m、0～10 m 这3 个影响区间，并取最大值作为该区间的应力集中系数，实现对断层影响范围和应力集中系数的量化表征。

图3 工作面回采接近断层时煤体应力集中系数变化曲线Fig. 3 Change curve of stress concentration coefficient of coal body when mining face approaches fault

1.2.2 采空区因素

建立与图2 同尺寸但不含断层面的数值模型，在模型中部进行单工作面开挖，提取模型中部煤体应力集中系数变化曲线，如图4 所示。

图4 采空区形成后煤体应力集中系数变化曲线Fig. 4 Change curve of stress concentration coefficient of coal body after the formation of goaf

由图4 可知，在采空区外缘50 m 范围内煤体应力集中系数最大，为1.43～1.58，平均为1.5。因此确定采空区因素的影响范围为50 m，应力集中系数为1.5。

1.2.3 区段煤柱因素

采空区开挖平衡后，在不含断层面的数值模型中分别留设5～50 m 的区段煤柱开挖沿空巷道，提取不同区段煤柱宽度下煤体应力集中系数变化曲线，如图5 所示。由图5 可知，区段煤柱宽度不同，煤柱内部的应力集中系数也不同。当区段煤柱宽度W≤5 m 或W≥50 m 时，煤柱内峰值应力集中系数为0.7～1.1，接近原岩应力状态；当5 m＜W≤10 m时，煤柱内应力集中系数逐渐增大，最大值为1.3；当10 m＜W≤30 m 时，煤柱内应力集中系数最大，最大值为2.2；当30 m＜W≤50 m 时，煤柱内应力集中系数逐步减小，应力集中系数最大值接近1.5。根据区段煤柱宽度划分4 个影响区间，并确定了不同区间的应力集中系数，实现对区段煤柱影响范围和应力集中系数的量化表征。

图5 不同区段煤柱宽度下煤体应力集中系数变化曲线Fig. 5 Change curve of stress concentration coefficient of coal body under different coal pillar widths

1.2.4 终采线外错因素

在不改变数值模型初始应力条件的情况下，以某矿31102、31103 工作面终采线外错条件为例进行模型开挖，工作面开挖后，煤体垂直应力分布云图如图6 所示。

图6 终采线外错区域煤体垂直应力分布Fig. 6 Vertical stress distribution of coal body in staggered area outside terminal line

由图6 可知，在终采线外错区域形成“三角板”状拐角煤柱应力集中区。提取平行于切眼方向一边的煤体应力绘图分析，如图7 所示。可看出距离外错点50 m 范围内应力集中程度较高，应力集中系数普遍大于1.5，距离外错点30 m 范围内存在应力峰值，最大应力集中系数为2.1。因此，以外错点为中心，划分终采线外错区域的影响范围为50 m，并根据与外错点距离的不同进行量化取值，实现对终采线外错因素的量化表征，如图8 所示。

图7 终采线外错区域煤体应力集中系数变化曲线Fig. 7 Change curve of stress concentration coefficient in staggered area outside terminal line

图8 终采线外错影响区域Fig. 8 Staggered area outside terminal line

1.2.5 巷道交叉因素

在不改变数值模型初始应力条件的情况下，分别开挖“T”型和“X”型交叉巷道，得到煤体垂直应力分布云图，如图9（a）、图9（b）所示，绘制沿交叉点对称的应力集中系数曲线，如图9（c）所示。

图9 “T”型和“X”型交叉巷道区域数值模拟结果Fig. 9 Numerical simulation results of "T" and "X" shaped cross-roadway

根据图9 分析结果对巷道交叉因素进行量化表征：“T”型和“X”型交叉巷道影响范围均为巷道交叉点前后30 m，其中“T”型交叉巷道影响区最大应力集中系数为1.05，“X”型交叉巷道影响区最大应力集中系数为1.1。

对于力学机制复杂、受目前研究所限无法精确模拟的主控因素，采用工程类比法确定其影响范围和应力集中系数。11 种冲击地压危险主控因素的量化表征结果见表1。

表1 冲击地压危险主控因素的量化表征结果Table 1 Quantitative characterization results of main controlling factors of rock burst hazard

在分层提取出9 项局部影响因素的识别主体后，对识别主体进行图形膨胀，划分主控因素的影响范围，得到9 项局部影响因素的影响图层。

为实现基于OpenCV 的图像处理，采用Python语言编程，利用.dxf 解析工具读取CAD 文件内各影响因素的图形信息，包括工作面位置、褶曲、断层、采空区、巷道交叉等图元信息数据库，并自动将CAD 文件转换成.png 格式的评价底图，图片为固定分辨率（DPI=500）。采用固定分辨率的方式进行图像精度控制，11 种因素影响范围的划分基于同一底图，保证绘图尺度的一致性。对于输入的PNG 底图，首先对工作面位置进行图像识别，得到特征线段，计算线段中的像素点个数M，输入工作面实际宽度N，得到图上像素距离X与实际距离D的转换公式：

在建立主控因素量化表征关系的基础上，对单一因素所产生的应力集中系数进行叠加，可得到工作面不同区域的总应力集中系数分布，基于应力集中系数对工作面冲击地压危险性进行判别评价，最终实现危险区域的分级分区划分。

以中天合创能源有限责任公司葫芦素煤矿21104工作面为例，进行冲击地压危险主控因素的量化表征与图像识别。21104 工作面为临空工作面，工作面宽度为320 m，走向长度为3 015 m，所回采2-1 煤层经鉴定具有强冲击倾向性，煤层单轴抗压强度为28.519 MPa。工作面内煤层平均埋深为646.4 m，与相邻采空区之间留设35 m 区段煤柱，面内发育有B4 背斜（倾角＞15°）及F9，F10，F43，F44，F50，DF20（落差＞10 m）断层。基于Python 语言编程读取采掘工程平面图内各影响因素的图形信息后，利用OpenCV 机器视觉库，得到21104 工作面褶曲、断层、采空区、区段煤柱、巷道交叉、初次来压、“见方”破断7 个局部影响因素的量化表征图层，如图10所示。

图10 21104 工作面冲击地压危险主控因素的量化表征与图形识别Fig. 10 Quantitative characterization and graphic identification of main controlling factors of rock burst hazard in 21104 working face

1.3 灰阶映射模型构建与危险区域划分

基于图像识别量化表征得到的单一因素影响图层以灰度图的形式存储在计算机内，为实现各项冲击地压危险影响因素叠加计算，需要进一步构建图像灰阶与应力集中系数的映射模型。

若将灰度图视为m行n列的矩阵，则灰度图中第m行、第n列像素点的灰阶可表示为lmn，第i（i=1，2，···，x，x为影响因素的数量）项影响因素对应的灰阶矩阵为

灰阶矩阵Li与其代表的应力集中系数矩阵 λi的映射关系为

式中：k为灰阶与应力集中系数的映射系数；λmn为第m行、第n列像素点对应的应力集中系数。

在得到各项影响因素对应的灰度图后，通过灰阶运算，可实现应力集中系数的叠加。如果应力增量采用线性累加方式进行叠加，全部影响因素的总灰阶矩阵L可通过下式计算。

如果应力增量采用非线性累加方式进行叠加，全部影响因素的总灰阶矩阵L可通过下式计算。

为进一步明确各区域的冲击地压危险程度，定义冲击发生的最小应力集中系数λmin及临界应力集中系数λlim，将总灰阶矩阵进一步划分为“无、弱、中等、强”4 个等级。根据工程经验和冲击地压发生的应力准则[20]，当煤体中的垂直应力超过煤层单轴抗压强度的1.5 倍时，煤体就具备了发生冲击地压的应力条件；当煤体中的垂直应力为煤层单轴抗压强度的2.5 倍时，冲击地压发生的危险性急剧升高。据此确定λmin及λlim取值。

式中：σz为煤层单轴抗压强度；γ为上覆岩层容重；h为采深。

采用min-max 标准化法，参照行业内冲击地压危险区域的划分方式，对λmin到λlim范围内的应力集中系数进行标准化，得到“无、弱、中等、强”4 种冲击地压危险等级下应力集中系数与灰阶的取值，见表2。

表2 冲击地压危险区域的划分标准Table 2 Classification standard of rock burst hazard area

为实现冲击地压危险区域叠加结果的可视化，通过Python 语言编程定义色相，将“无、弱、中等、强”冲击地压危险区域分别用从冷到暖的渐变色表示，并通过3 条危险等级划分线进行色彩区分，3 条危险等级划分线的灰阶分别为kλmin，k(λmax-λmin)/3，2k(λmax-λmin)/3。冲击地压危险区域可视化评价结果如图11 所示。

图11 冲击地压危险区域可视化评价结果Fig. 11 Visual evaluation results of rock burst hazard area

在最终输出的评价云图中，映射系数k控制云图绘制的精细程度。k越大，云图显示的冲击地压危险区域越精细，但可显示的应力集中系数范围变窄；k越小，可显示的应力集中系数范围越大，但云图的精度降低。通过调整映射系数k，可实现冲击地压危险评价结果的动态调整。

2 冲击地压危险智能化评价软件开发与应用实例

2.1 软件开发

冲击地压危险智能化评价软件采用生命周期的软件开发方法，以及自顶向下、逐步求精的设计方法。软件开发语言为Python，HTML，CSS，JavaScript。软件主要由CAD 文件解析模块cad2pic、危险因素识别叠加模块fc_plus、后端管理模块Django、前端操作页面模块html+Layui 及MySQL 数据库组成，如图12 所示。

图12 软件组成Fig. 12 Software composition

软件核心模块由CAD 文件解析模块cad2pic 和危险因素识别叠加模块fc_plus 组成。其中，CAD 文件解析模块基于ezdxf 库开发，实现对.dxf 文件格式的采掘工程平面图进行解析，并将图层转换为.png格式图片。危险区域绘制模块主要功能：① 识别并定位图片中的冲击地压危险因素位置。② 逐个计算冲击地压危险因素的影响范围和程度，计算单一影响因素的灰阶矩阵。③ 叠加指定因素，生成总灰阶矩阵。④ 生成冲击地压危险区域评价云图，划分危险区域。通过前端操作页面模块完成参数输入、文件上传、叠加因素选择等操作，后端管理模块及MySQL 数据库进行数据的存储、管理及后台计算。

2.2 应用实例

冲击地压危险智能化评价软件通过流程化设计实现冲击地压危险性的智能自动评价，以葫芦素煤矿21104 工作面为例进行说明。

软件共设计文件上传、图层选择、叠加因素选择、图片编辑及导出4 个界面，如图13 所示。在文件上传界面，需输入工作面长度、煤层埋深、初次来压步距、单轴抗压强度、矿震能级5 个参数，并上传.dxf 格式的采掘工程平面图。在图层选择界面，软件将自动读取上传的CAD 文件，并将图层以列表形式展示，勾选巷道图层、断层、褶曲等图层，可生成CAD 底图，用于冲击地压危险影响因素的识别与危险区域划分云图的叠加。在叠加因素选择界面，可手动选择自重应力等8 项静载影响因素、初次来压等3 项动载影响因素，并提供了“叠加”与“叠乘”2 种应力增量的叠加方式，选择后软件自动计算并生成冲击地压危险区域评价云图，并将冲击地压危险区域划分为“无、弱、中等、强”4 个等级。在图片编辑及导出界面，可手动添加文字、工作面回采方向等简单标记，并对云图进行缩放、裁剪、旋转等基本操作，最后导出调整完毕的高清图片。

图13 冲击地压危险智能化评价软件界面Fig. 13 Interface of intelligent evaluation software of rock burst hazard

与传统的冲击地压危险性评价方法相比，本软件将采掘、地质、危险区域划分等冲击地压要素集成于一张图显示，将传统仅针对巷道的一维线性危险区域划分方法改进为针对整个采掘空间的二维平面划分方法，评价结果所含信息更多、精细度更高。同时，所有前后处理工作可在10 min 之内完成，大幅提升了评价效率，降低了人工成本。

3 应用效果检验

冲击地压危险智能化评价软件已应用于多个矿井采掘工作面的冲击地压危险性评价。以葫芦素煤矿21104 工作面为例，对软件应用效果进行检验。

微震是检验区域冲击地压危险性最常用、最有效的手段之一[21-22]。采掘区域内微震活动性越强、能量释放越高，表明区域内冲击地压危险程度越高。为定量评价采掘区域内微震能量集中程度，首先定义微震能量密度指标U：

式中：Ej为落入统计区域内第j个矿震的能量；S为统计区域面积。

对21104 工作面回采前1 000 m 范围内103J 以上大能量微震事件进行平面定位，结果如图14（a）所示。将工作面划分为10 m×10 m 的统计网格，绘制21104 工作面全部微震事件的能量密度云图，如图14（b）所示。冲击地压危险智能化评价软件评价结果如图14（c）所示。

图14 冲击地压危险性评价结果与微震能量密度云图对比Fig. 14 Comparison between rock burst hazard assessment results and cloud picture of microseismic energy density

由图14 可知，大能量微震事件出现的区域对应高能量密度区，高微震能量密度区与划分的中等以上冲击地压危险区域对应关系较好。主要表现在：① 21104 回风巷一侧微震事件整体多于运输巷一侧，与冲击地压危险区域划分结果一致。② 21104工作面在F9，F10 断层及B4 背斜与区段煤柱相交区域的微震事件能量释放远高于其他区域，表明该区域具有较高的冲击地压危险性，与冲击地压危险区域划分结果一致。③ 21104 工作面回采至450～500 m时，出现沿工作面倾向的初次“见方”影响区，与冲击地压危险区域划分结果相吻合，但“见方”位置略超前于评价结果。

为进一步对比软件评价结果与现场实际矿压显现之间的关联性，统计了21104 工作面现场出现中等及以上冲击地压危险预警，以及工作面大煤炮、煤体弹射、巷道变形量大等矿压异常区域的分布情况，如图15 所示。

图15 21104 工作面冲击地压危险预警及矿压异常区域分布Fig. 15 Distribution of rock burst hazard early warning area and abnormal area of mine pressure in 21104 working face

由图15 可知，21104 工作面在回采至310～990 m时频繁发生预警。据统计，现场矿压异常区主要集中于回风巷一次“见方”附近100 m、回风巷与B4 背斜、F9 断层交叉区域200 m，以及二次“见方”附近160 m 等区域，有明显大煤炮，局部存在支护失效、煤体弹射等动力现象，工作面具有强冲击地压危险性。软件划分的临空巷中等冲击地压危险区及局部的强冲击地压危险区与现场一致性较好。

4 结论

1） 通过半定量化估算确定了11 项冲击地压危险动静载主控因素的图形量化表征方法，并基于OpenCV 机器视觉库和深度学习模型，实现了对单一主控因素的图像识别。

2） 构建了图像灰阶与应力集中系数的映射矩阵，实现了单一影响因素的线性与非线性叠加，得到评价区域的应力集中系数矩阵，明确了“无、弱、中等、强”冲击地压危险区的划分标准。

3） 开发了冲击地压危险智能化评价软件，采用微震等手段对评价结果进行了检验，评价结果与现场实测矿压规律一致性较高。该软件提高了评价效率与冲击地压危险区域划分精度。