煤矿采空区覆岩“三带”智能识别方法

李建文，赵 文，吴振坤，徐小兵，王庆涛，段隆臣

(1.中煤科工生态环境科技有限公司，北京 100013；2.中国地质大学(武汉) 自动化学院，湖北 武汉 430074；3.中国地质大学(武汉) 工程学院，湖北 武汉 430074)

煤矿采空区治理是践行矿山绿色发展理念的重要举措，而对采空区的覆岩裂隙发育情况进行勘探是开展采空区治理的基础条件。煤层开采后由于岩层移动、变形和破坏，在覆岩中形成采动裂隙，导致矿区的地质结构遭到破坏，极易引起地表塌陷、裂缝和下沉等地质灾害，严重影响城市居民生活与建设发展。因此，探明采空区地下裂隙分布及覆岩破坏情况，对重新利用采空区场地和满足城市建设开发需求有重要意义。

“三带”型是一种典型的采空区覆岩发育破坏类型。根据覆岩中的裂隙发育情况，自下而上可分为垮落带、断裂带和弯曲带，简称“三带”。其中，断裂带和垮落带合称导水裂隙带，在《煤矿采空区岩土工程勘察规范(2017 年版)》(简称勘察规范)中也称垮落断裂带[1]。

目前确定“三带”分布与高度的研究主要包括实测“三带”裂隙发育程度和预测导水裂隙带高度两方面。其中，常见的实测手段包括工程物探[2-4]、工程钻探[5]、煤田测井[6]、钻孔成像[7-8]、井下注水观测[9-10]等。实测法是最可靠的“三带”识别方式，但总体成本较高，并需要专门的技术人员对数据进行分析才能得到采空区覆岩裂隙发育情况。预测导水裂隙带高度的研究方法主要包括经验公式法[11-13]、相似材料模拟法[14]、应力计算法[15-18](也指数值模拟法)，以及机器学习法[19-20]等。其中，勘察规范中的经验公式应用较为广泛，可以得到断裂带和垮落带的发育高度范围，但不能得到垮落带和断裂带的具体发育高度和“三带”的具体分布深度，需要结合实测方法才能确定[21]。通过预测导水裂隙带高度确定“三带”分布与高度的主要优点是应用成本低，对实际工程有一定指导意义，但大多基于室内试验获得，因而不能完全反映实际情况[22-23]。

随着采空区治理自动化、智能化的发展，将物联网技术、数据融合和分析技术进行有机结合，从而预测断裂带高度，为钻进过程中自动获取“三带”高度提供了新的思路[24]。目前钻探现场可以安装各种类型的传感器，自动记录钻进时间、钻进深度及冲洗液漏失量等钻进参数，但缺少自动分析数据的功能。技术人员完钻后对数据进行统一分析时具有滞后性，且分析结果受技术人员主观判断影响显著。此外，施工现场对采空区的初步判别需要技术人员持续关注钻进参数的变化，而当出现采空区误判或漏判时会则影响下一步的实时决策，例如针对可能出现的孔内事故采取预防措施和确定完钻孔深等。

综上所述，本文结合实测“三带”裂隙发育程度和预测导水裂隙带高度两类研究方法的特点，基于工程钻探中实时获取的钻进数据，采用贝叶斯变化点检测(Bayesian Online Changepoint Detection，BOCD)[25]算法在钻进过程中自动分析数据，并融合勘察规范中的经验公式，实现在钻进过程中对覆岩“三带”的智能识别。该方法可充分发挥钻进数据的时效性，提高“三带”划分效率和准确性，并同时提高煤矿采空区治理的智能化水平。

1 “三带”的钻进参数响应特征

图1 为采空区覆岩“三带”的剖面。在工程钻探中，冲洗液漏失量和钻速的变化特征与“三带”有很大相关性[7,26-28]。

图1 采空区覆岩“三带”剖面[7]Fig.1 Schematic cross section of overlying strata three zones of goaf[7]

弯曲带为断裂带顶部到地表之间的地层，基本呈整体移动，受采动影响不大。在垂直剖面上，地层上部很少出现离层裂隙；相比而言，地层下部更常见离层裂隙，但这些离层裂隙仅局部充水，且不与断裂带连通。因此，在弯曲带中钻进时，冲洗液漏失量和钻速总体趋于平稳，局部存在轻微波动。

断裂带位于弯曲带下方、垮落带上方。其上部有裂隙，可以间接导水和积水，但因垂直裂隙不发育，一般不与下部裂隙沟通；其下部垂直裂隙逐渐发育增强，离层裂隙和垂向裂隙连通，导水性明显增加，并能向下渗流至采空区。因此，钻进至断裂带时，冲洗液漏失量会突然大幅增加，相较于正常地层可能增加数倍。对于垂直裂隙发育的地层，冲洗液甚至会全部漏失。与此同时，由于裂隙发育，钻速会大幅提升，并存在较大波动。

垮落带位于采空区覆岩最下部，通常为较大空洞或充填有垮落的岩块，且岩块之间存在很多空隙，连通性较强。因此，大多数“三带”发育的采空区钻孔在钻进至垮落带时，冲洗液漏失量很大，甚至全部漏失，而钻速相较于正常地层也会显著提高，且可能存在明显波动。

当钻孔进入采空区底板后，由于地层受采动影响不大、相对完整，钻速相较于垮落带会明显降低，且逐渐趋于稳定。

2 “三带”识别方法

通过“三带”钻进响应特征分析可知，从弯曲带进入断裂带时冲洗液漏失量和钻速均会明显增加，而进入采空区底板时钻速会明显降低。因此，将这两个特征的提取看作对变化点的识别，并采用BOCD 算法来识别弯曲带下限和采空区底板。

然而，仅通过变化点检测无法获得实际的“三带”界限，主要原因为：(1)识别变化点与目标变化点为多对一，需要从多个变化点中筛选出目标变化点；(2)通过钻进参数仅能确定弯曲带下限深度、采空区底板深度和垮落带下限深度，无法确定断裂带下限深度。因此，本文在BOCD 算法识别变化点的基础上，结合勘察规范中有关断裂带和垮落带高度的经验公式，以及相关钻前资料，以此自动识别准确的“三带”界限，并计算得到“三带”分布的高度值。

2.1 BOCD 算法原理

BOCD 算法基于贝叶斯框架，可以实时识别数据随时间的变化情况，目前已被广泛应用于金融分析、生物识别等领域。BOCD 算法的核心思想是：假设时间序列x1:t(x1:t为序列x1,x2,···,xt的集合，t∈N*为序列中数据点所对应的时间或序号)中任意两个相邻变化点之间的数据相互独立，且分布类型相同，则新观测数据xt+1是否为变化点的判断依据是该数据的预测条件概率值P(xt+1|x1:t)的大小。若新观测的数据与上一变化点之后的数据非常相似，则其预测条件概率值将会很高，即新的观测数据不是变化点。反之，其预测条件概率值会很低，即新的观测点为变化点。

定义行程长度rt为从上一变化点到t时刻走过了r个单位时间，表示t时刻行程长度为r时的数据集合，则P(xt+1|x1:t)的表达式为：

将风险函数f看作离散指数分布，即f=1/λ，其中λ为尺度参数，与数据的采样频率有关。P(rt,x1:t)可以通过上一时刻的概率迭代计算得到。

2.2 “三带”智能识别方案

假设待识别钻孔的设计孔深为d，垮落带、断裂带分布在设计孔深往上Hu之内，候选变化点的有效识别深度区间为[d-Hu,d]，所有“三带”界限深度均位于该区间内。Hu可按如下公式计算：

Hli与采空区顶板倾角α、采空区顶板岩层的碎胀系数k等有关[1,7]，该公式给出了Hli的最大和最小值，为了扩大候选变化点的有效深度区间，尽量覆盖更多的候选点，这里取Hli的最大值；与Hli的取值类似，这里取最大值，取最小值；M为已知值，Hn为设计值。

由“三带”钻进中的冲洗液漏失量和钻速变化特征可知，区间[d-Hu,d]内的第一个变化点所对应的深度即为弯曲带下限深度d1；最后一个变化点所对应的深度即为采空区底板深度d4，则垮落带下限深度d3为d4-M。

上述界限深度确定后，接下来需要确定断裂带下限深度d2。由于断裂带下限深度附近可能存在多个候选变化点，应用垮落断裂带和垮落带高度的经验公式[1]，结合弯曲带下限深度d1和垮落带下限深度d3，由下式计算断裂带下限深度d2的识别深度范围[d2]：

式(5)中，括号内的值为估算的Hm变化比例范围，正负号分别代表范围的上下限；寻找距离[d2]区间中心最近的变化点，即为断裂带下限深度d2；若[d2]区间内无变化点，则取该区间的中心点为断裂带下限深度d2。

根据上述“三带”界限识别结果，可得断裂带实际高度为d2-d1、垮落带实际高度为d3-d2。至此，“三带”的界限深度和高度已经全部识别完成，完整的识别流程如图2 所示。

图2 “三带”识别流程Fig.2 Three zones identification flow chart

按照图2 流程，弯曲带下限深度、采空区底板深度可以实时确定。在通过采空区底板确定垮落带下限后，结合勘察规范中经验公式的结果可以确定断裂带下限深度。当采空区底板深度确定后，即可确定完钻深度。在实际工程中，该流程中涉及到的钻前地质资料等参数的值(如M、d、k、α、Hn等)可以在钻进前通过工控机或物联网设备输入给定，钻进时无需其他辅助即可自动识别“三带”结果。需要指出的是，上述“三带”识别方法是基于回转钻进的钻进参数特征实现的，如果钻进方式发生明显改变，则其钻进过程中的“三带”识别还需要进一步开展相关适配性研究。

上述BOCD 算法和“三带”识别过程可以通过Python 代码实现，并将其作为Web 应用部署在服务器上，从而通过无线网络分别与数据感知端和用户应用端相连，实现“三带”识别的物联应用。“三带”识别的Web 应用如图3 所示，数据感知端包含不同的钻井监测站点，每个站点由各类传感器(如钻速传感器、流量传感器等)组成，用于获取基础数据。Web 服务器是主要功能的实现载体，如数据存储功能、“三带”识别功能等，其通过无线网络接收数据感知端的数据和将识别结果等数据传输到应用终端。应用终端直接面向用户，包括移动终端、大屏等，用于显示识别结果或输入数据。

图3 “三带”识别的Web 应用Fig.3 Illustrating the application of three zones identification for goaf in Web applications

3 工程应用

3.1 项目概况

为了验证所提方法的应用效果，选择山东某矿区一口勘察井的监测数据进行测试。矿区采用立井分水平开拓(主井、副井、风井)，中央并列式通风，条带短壁式综合机械化一次采全高采煤工艺，顶板自由垮落法管理[7]。测试钻孔的采空区顶板为砂质泥岩、泥岩，单轴饱和抗压强度平均值分别约为40.77 和31.66 MPa，属于较硬岩。煤层平均倾角α=6°，顶板岩层的碎胀系数为k=1.35，平均采高为M=2.8 m。钻孔为垂直钻孔，设计孔深d=570.00 m。

该项目所在勘察区使用回转钻进工艺，并采用PDC 钻头作为碎岩工具。测试钻孔只统计纯钻进时间内的冲洗液漏失量和钻速，其中冲洗液漏失量通过流量计测量进浆量和出浆量得到，深度和钻速通过无线射频位移传感器得到。从钻进孔深为266.00 m 时开始，每钻进约0.5 m 记录一次数据，钻进至孔深565.87 m 时停止记录，共657 组数据。

图4 为现场获取的冲洗液漏失量和钻速数据所绘制的曲线，两条曲线均可以划分为两段。在第一段曲线中，冲洗液漏失量整体数值在110 L/m 左右，部分数据的值较高，主要是由于监测误差和这一段地层中存在的破碎带造成；钻速曲线整体数值在11 mm/min 左右，部分数据的值较低，考虑主要是由于监测误差造成。在第二段曲线中，冲洗液漏失量和钻速值都明显高于第一段曲线，这主要与钻头进入垮落带和断裂带有关。

图4 原始钻进监测数据Fig.4 Raw drilling monitoring data

3.2 BOCD 变化点识别

利用钻进数据识别“三带”时，同一采样频率的λ大致相同，且采样频率越大，λ也越大。实际应用中可通过多次测试确定某一采样频率对应的λ值，本文中的λ取值为50。给定上述条件，即可以识别钻进数据的变化点。

在识别前，为了提高识别效果，选择中值滤波的方法降低异常数据对变化点识别结果的干扰。中值滤波窗口大小设置为5。图5 为BOCD 算法的运行结果，每条垂直虚线都与一个行程的开始对应，即对应一个变化点，其中图5a 和图5b 分别为冲洗液漏失量和钻速滤波后的曲线及对应的变化点识别结果。与原始数据相比，滤波后冲洗液漏失量和钻速的第一段曲线更加平稳，基本没有异常值，第二段曲线同样波动减小，并且曲线的变化特点更加明显，表明滤波可以有效地去除异常值。如图5 所示，冲洗液漏失量曲线识别出来的变化点对应深度为494.22、511.36、515.36、543.50 m，钻速曲线识别出来的变化点对应深度为269.72、505.68、510.42、515.71、544.26、562.87 m，所有这些深度值都是用来确定“三带”界限的候选变化点。

图5 滤波后钻进监测数据变化点识别结果Fig.5 Identification results of Bayesian Online Changepoint Detection from the drilling monitoring data after filtering

3.3 “三带”识别结果与分析

试验钻孔所在地区已探明“三带”的钻孔中，已知实际垮落断裂带高度的最大值为87.30 m，设计钻孔要求钻入采空区底板深度Hn不少于10 m，采高M为2.8 m。

由给定的钻前地质资料可知，采空区顶板属于较硬岩，则根据勘察规范附录L 可以确定试验钻孔垮落断裂带高度Hli、垮落带高度Hm的经验计算公式如下：

取 max(Hli)为61.14 m。本实例中，区域内只有厚度为2.8m 的单一煤层，因此根据规范计算获得的Hli和相同。取min为43.34m，由式(4)可得Hu为135.95 m，则候选变化点的有效深度区间为444.04～570.00 m。这一范围内的第一个变化点所对应深度为494.22 m，则弯曲带下限深度d1为494.22 m。这一范围内的最后一个变化点对应深度为562.87 m，则采空区底板深度d4为562.87 m。进而确定垮落带下限深度d3为560.07 m。

由式(5)可得d2的识别深度范围为[550.17，552.65]，该界限内没有变化点，则断裂带下限深度d2取识别区间中心点所在位置，即551.41 m。

至此，已得到“三带”的下限深度自上而下分别为494.22、551.41、560.07 m。为评价本文识别方法的准确性，选择与本文依托相同项目和测试钻孔的文献[7]中的“三带”结果作为对比。该文献综合利用工程地质钻探、孔内电视与煤田测井3 种方法，得到测试钻孔的“三带”界限深度自上而下分别为493.55、551.10、559.55 m。图6 为智能识别得到的“三带”结果(灰色垂直虚线)与文献[7]所得结果(黑色垂直虚线)的对比情况，可以看出两者较吻合。“三带”界限深度的智能识别结果与实际结果对比见表1，误差均小于1 m。“三带”高度的智能识别结果和经验公式计算结果与其实际高度的对比见表2。其中，智能识别结果的高度误差均小于3%，而经验公式计算得到的断裂带和垮落高度误差分别为9.23%和4.85%。结果表明，“三带”智能识别方法显著优于经验计算方式，且具有较高的识别精度。

表1 识别“三带”界限与实际“三带”界限[7]对比Table 1 Comparison of the intelligent identification results and actual results[7] boundary depths in the three zones of goaf

表2 “三带”高度结果对比Table 2 Comparison of heights in the three zones of goaf

图6 “三带”识别结果与实际人工划分结果[7]对比Fig.6 Comparison of the intelligent identification results and actual results[7] in the three zones of goal

应用结果表明，在变化点识别的基础上，结合勘察规范中垮落断裂带和垮落带的经验公式，可将断裂带下限的范围由弯曲带下限和垮落带下限之间缩到更小的范围，使得识别结果更有针对性。因此，融合了BOCD算法的采空区“三带”智能识别方法比传统方法更接近实际情况。

4 结论

a.提出的基于冲洗液漏失量和钻速等钻进数据的变化点检测方法，可以实现在钻进过程中的“三带”智能识别，无需等完钻后再进行数据分析，充分发挥了钻进数据的时效性，有效避免了技术人员的主观判断和经验对结果产生的影响。

b.与实际多种勘察手段综合划分的结果对比表明，智能识别的“三带”深度误差均小于1 m，高度误差均小于3%，具有较高的识别精度。将实测方法和预测方法相结合，将冲洗液漏失量和钻速这两个与“三带”有密切关系的数据作为输入参数，提高了识别的精度。提出的“三带”智能识别方法可以代替多种手段综合划分的方式，有助于提高“三带”划分效率，降低采空区勘察成本。

c.提出的“三带”智能识别方法适用于开采方式、地质条件、煤层厚度、底板位置等相关信息比较明确的采空区勘探。现有监测参数仅考虑了钻速和冲洗液漏失量，且其响应特征主要针对以冲洗液为循环介质的回转钻进，未来可进一步研究其他钻进参数及其在不同钻进方式中与“三带”之间的响应关系，以提高智能识别方法的适用性。

符号注释：

d为设计孔深，m；d0为地表深度，m；d1为弯曲带下限深度，m；d2为断裂带下限深度，m；d3为垮落带下限深度，m；d4为采空区底板深度，m；f为风险函数；Hn为待识别钻孔规定钻进采空区底板的深度，m；Hli为按勘察规范经验公式计算得到的待识别钻孔垮落断裂带高度，m；为区域内已探明的实际垮落断裂带高度，m；为按勘察规范经验公式计算得到的已完成钻孔垮落断裂带高度，m；Hm为按勘察规范中经验公式计算得到的待识别钻孔垮落带高度，m；Hu为待识别钻孔弯曲带下限与设计孔深的最大距离，m；k为采空区顶板岩层的碎胀系数；max(Hli)、min(Hli)、ave(Hli)分别为按经验公式计算得到的待识别钻孔垮落断裂带的最大高度、最小高度和平均高度，m；为实际垮落断裂带高度的最大值，m；min()为按勘察规范中经验公式计算得到的已完成钻孔垮落断裂带高度的最小值，m；M为采高，m；N*为正整数集；P(rt|rt-1)为变化点先验概率；rt-1、rt为行程长度，分别表示从上一变化点到t-1、t时刻走过的单位时间；t为观测数据点所对应的时间或序号，从1 开始；x1:t、x1:t-1分别为从第一个时间单位到第t、t-1 个时间单位的观测数据集；xt+1、xt分别为t+1、t时刻的观测数据；为t时刻对应行程长度为r内的数据集合；α为采空区顶板倾角，(°)；λ为离散指数分布的尺度参数。