基于聚类的煤矿井下钻孔瞬变电磁异常响应边界成像方法

2022-08-09张幼振樊依林李宇腾

煤田地质与勘探 2022年7期

范涛，李萍，张幼振，赵睿，房哲，樊依林，刘磊，王程，李宇腾

(1.中煤科工集团西安研究院有限公司，陕西西安 710077；2.陕西省重点科技创新团队(地球物理探测技术与装备创新团队)，陕西西安 710077；3.煤炭行业工程研究中心(物探技术与装备)，陕西西安 710077)

近10 年煤矿事故总量逐年下降，但事故下降幅度趋缓。较大事故总起数、死亡人数总体呈下降趋势，但从2013 年开始曲线斜率明显趋平，2016 年后呈锯齿形下降，2017 和2019 年等一些年份反弹，2020 年再度下降。重特大事故呈锯齿形下降，同时波动幅度较大，2018 年后出现反弹，2020 年重特大事故总量与2019 年持平。事故中水害较大以上事故上升幅度大，2020 年发生较大以上水害事故3 起、死亡21 人，同比增加1 起、多死亡12 人，分别上升50%和133.3%，其中，较大水害事故2 起、死亡8 人，同比起数持平、少死亡1 人，重大事故1 起、死亡13 人，同比增加1 起、多死亡13 人，占全国煤矿较大以上事故起数和死亡人数的37.5%，经济损失和社会影响非常严重[1]。而掘进工作面则是煤矿水害事故发生最多的地点[2]。因此，掘进工作面前方隐伏水害超前探测是亟待解决的技术难题。

水害隐患超前探查代表性的方法主要是瞬变电磁法、直流电法等电磁类方法[3-5]，但是，煤矿井下探测装备功率受煤矿本安防爆限制，且井下环境的电磁干扰较大等因素的存在，致使矿井电磁类方法的探测距离短、精度偏低、多解性强[6-8]。为了解决煤矿井下掘进工作面前方的探测深度与探测精度的矛盾，许多学者逐步开始利用井下掘进工作面的钻孔进行瞬变电磁探测工作，该方法可以在掘进前开展远距离、高精度的隐伏水害超前预报[9-12]。

钻孔瞬变电磁的数据处理一般采用电阻率反演成像方法，考虑到反演存在多解性，因此，反演拟合过程中为了提高精度，往往会选择较为光滑的模型拟合实测数据，使得最终反演结果是一个电阻率连续光滑渐变的成像模型，在地质体边界处的对比程度和变化情况较为模糊，很难清晰地反映地质异常体响应与背景值的差异，对异常体规模、形态的解释工作常常需要经验丰富的专家人为干预，也不利于生产中准确指导物探工作之后的钻探和掘进工作。鉴于此，有必要研究一种提高成像结果中电阻率值聚合度，进而突出电性边界的成像方法。

显然，对电阻率的分区聚合属于典型的分类问题，可以选用无监督机器学习方法进行解决。数据挖掘领域经常使用无监督机器学习算法，主要是采用它来发现大量无标签数据的分布规律，实现对数据的区分或分类。经常使用的无监督机器学习方法主要有局部线性嵌入算法、主成分分析、局部切空间排列算法、拉普拉斯特征映射、等距映射和应用最多最广的聚类算法[13]。

聚类算法是指基于一定的优化标准将一堆无标签数据对象自动划分成若干类别的方法，这个方法要求同一类别的数据具有相似的特征，不同类别的数据具有不同的特征[14]。

近年来，聚类方法已经广泛被应用于地球物理领域的数据处理和特征挖掘中。尤其在地震勘探领域的应用较多，王伟涛等[15]通过层次聚类分析对汶川大地震的余震序列中的近似地震和重复地震进行了有效辨识；张岩等[16]采用结构聚类字典学习方法进行了地震数据随机噪声压制方面的研究；S.Scitovski[17]采用基于密度的聚类方法对地震记录进行了不同地震类型的划分。在重磁资料的处理解释中，张新兵等[18]基于改进的K-Means 聚类分析方法实现了重磁局部异常的自动圈定；李斐等[19]在优化不同区域的重力观测密度方面应用了聚类分析方法；曹书锦等[20]开展了自适应模糊聚类对多异常源的精准确定工作。在电磁数据处理解释领域，杨生等[21]在大地电磁曲线类型划分中使用了聚类方法，对削弱地质推断的多解性有较好的作用；SongYuchen 等[22]提出一种应用自适应密度聚类分类电测深曲线类型的方法；M.Audebert 等[23]结合了K-Means 聚类和电阻率CT 成像，提出一种多次反演解释垃圾场渗流区的有效方法；李晋等[24]提出一种联合聚类及递归的信噪比分离方法，对MT 低频部分的数据质量有所改善。

本文参考以上资料，考虑不同的地质体存在电阻率差异，基于统计思想提出将具有相近电阻率值的成像元素划分至一个类别，进而实现快速识别地质异常响应边界的方法。本文讨论了选用聚类方法的原则，介绍了聚类的方法原理，研究了确定最佳聚类条件的方法，最后使用三维数值模拟和井下实测数据检验了方法的有效性和实用性。

1 方法原理

1.1 聚类方法的选取

对电阻率进行聚类本质上是对一个一维数据体进行聚类，相当于对电阻率进行自适应层级划分。绝大多数聚类算法都是针对二维及以上的数据，针对一维数据的聚类算法非常少，主要有K-Means 方法和Jenks Natural Breaks 方法。

K-Means 算法，也被称为K-平均或K-均值算法，是一种广泛使用的聚类算法。K-Means 算法以距离作为数据对象间相似性度量的标准，即数据对象间的距离越小，则它们的相似性越高，则它们越有可能在同一个类簇。之所以被称为K-Means 是因为它可以发现K个不同的簇，且每个簇的中心采用簇中所含值的均值计算而成。

Jenks Natural Breaks 算法，也就是自然断点分类，分类的原则就是将方差接近的放在一起，分成若干类。通过计算每类的方差和方差和，用方差和的大小来比较分类的好坏，值最小的就是最优的分类结果(但并不唯一)。

K-Means 和Jenks Natural Breaks 在处理一维数据时完全等价。它们的目标函数一样，但是算法的步骤不完全相同。K-Means 是先设定好K个初始随机点。而Jenks Natural Breaks 则是用遍历的方法，一个点一个点地移动，直到达到最小值。显然，n个数分成k类，Jenks Natural Breaks 就要从n-1 个数中找k-1 个组合。当数据量较大时，如果分类又多，计算量会显著增加，因此，一般针对一维数组的聚类，更多选用计算量小的K-Means 算法。

1.2 K-Means 聚类算法

K-Means 算法的思想很简单，对于给定的样本集，按照样本之间的距离大小，将样本集划分为K个簇。让簇内的点尽量紧密地连在一起，而让簇间的距离尽量大。

假设数据集为：

式中：n为数据个数；xi为除质心外的其他点。

划分的簇为：

式中：K为预先设置的分类数。

此时，最小化平方误差E即为：

式中：μi是簇Ci的均值向量，也称为质心，表达式为：

计算开始时，先采用随机方式选定K个质心对所有数据进行初始分类，共分为K个初始类别，之后对所有数据计算它们与质心之间的欧氏距离，再依据这些距离的平均值更新质心和划分类别，最后反复迭代该更新过程，直到满足迭代的停止条件为止。

迭代的停止条件一般是质心变化率满足下式：

1.3 聚类中关键条件的确定方法

K-Means 聚类算法需要提前给出的条件主要有3 个：初始质心、聚类时的距离计算规则和聚类数目。

由于K-Means 聚类易陷入局部极小，不同的初始质心可能会导致不同的结果，本文根据聚类数目K的不同，选择距离尽可能远的K个点为初始质心，具体做法为：随机选择一个点作为第一个初始簇质心，然后选择距离该点最远的那个点作为第二个初始簇质心，然后再选择距离前2 个点的最近距离最大的点作为第3 个初始簇的质心，以此类推，直至选出K个初始类质心。

考核类簇之间的相似性程度主要依靠各类簇之间的距离，较为常见的距离计算方法有欧氏距离、曼哈顿距离、夹角余弦距离和相关距离等。一般来说，欧氏距离最为简单直观，也更能反映数据在数值特征上的差别，因此，本文选择使用欧氏距离。

聚类数目是K-Means 聚类算法中最重要的参数，因为不准确的聚类数目会明显导致分类效果变差，但是目前并没有依托数学原理的完美评价标准，本文采用基于组内平方误差和(Sum of Squared Error，SSE)的肘部法则来确定最佳聚类数目。

组内平方误差和(SSE)表示一个类簇内各点与该类质心的平方误差之和，可由下式计算。

ESS越小则说明各个类簇越收敛，但显然不是越小越好。因为考虑一种极端情况：将所有的样本点均视作单独类簇，此时ESS为0，而并未达到分类的目的。因此，需要在聚类数目和ESS之间寻找一个平衡点。

肘部法则就是这样一种方法。指定一个j值，即可能的最大聚类数目。然后将聚类数目从1 开始一直递增到j，计算出j个ESS。随着聚类数目增多，每一个类簇中数据点数量越来越少，距离越来越近，因此，ESS值必然随着聚类数目增多而减少。但当ESS减少得很缓慢时，可以认为进一步增大聚类数目分类效果也并不能增强，这个“肘点(拐点)”就是最佳聚类数目。

1.4 异常响应边界成像步骤

应用K-Means 聚类方法实现对异常响应边界成像的具体步骤如下：

(1) 对实测数据的Z分量数据进行视电阻率计算或反演成像，获得(视)电阻率参数。

(2) 计算(视)电阻率参数不同聚类数目情况下的组内平方误差和(ESS)，应用肘部法则寻找最佳聚类数目。

(3) 按照最佳聚类数目对(视)电阻率参数进行聚类。

(4) 将同一类中的(视)电阻率值全部更新为该类的质心值，并与原空间坐标重新对应。

(5) 对新的聚类后电性参数数据文件进行成像。

2 数值模拟效果

为验证本文方法的成像效果，设计如图1 所示的三维模型，采用时域有限差分方法进行了数值模拟。在钻孔深度方向40 m 处，第一象限45°放置1 个规模为10 m×10 m×10 m 的低阻异常体，异常体中心点距离钻孔15 m，采样时间范围为1×10-6～6.136×10-4s，模型其他参数见表1。

图1 数值模型Fig.1 Schematic diagram of the numerical model

表1 模型参数Table 1 Model parameters

对26～54 m 范围内的测点Z分量数据采用晚期视电阻率计算和层厚累加法进行视电阻率成像，可以得到如图2 所示的沿钻孔方向的视电阻率剖面图，图中横坐标z为钻孔深度，纵坐标r为径向探测距离。可以较为清晰地看到在钻孔深度40 m、钻孔径向15 m 位置有较为明显的低阻异常响应(蓝、绿色区域)。

图2 数值模型视电阻率剖面Fig.2 Resistivity profile of the numerical model

如果将蓝色区域解释为异常，则异常响应明显小于实际异常体，二者面积比值约为0.624；如果将绿色区域解释为异常，则异常响应明显大于实际异常体，二者面积比值约为2.458；从图2 中无法清晰反映异常响应的边界，需要根据经验人为解释确定。

采用本文1.3 节肘部法则对模型视电阻率数据进行分析，可得到如图3 所示的ESS曲线图。由图3 可知，在聚类数目为2 时曲线出现肘点，因此，选择K=2，采用1.2 节中的K-Means 算法对模型视电阻率数据进行聚类处理，可得到图4 所示的异常响应边界成像结果。由图4 可知，在钻孔深度40 m、钻孔径向15 m 位置有明显的孤立低阻异常响应(墨蓝色区域)，异常响应边界非常清晰，受晚期视电阻率计算和图像插值算法的影响，异常响应形状表现为直径约10 m 的近似圆形，其与实际异常体(红色虚线)的面积比约为0.985，整体参数与模型设置参数吻合较好。

图3 数值模型ESS 曲线Fig.3 ESS curve of the numerical model

图4 数值模型异常响应边界成像结果Fig.4 Imaging result of anomaly response boundary of the numerical model

3 应用实例

为了验证本文异常响应边界成像方法对实际生产中煤矿积水采空区的实际解释能力，在陕北某矿开展了钻孔瞬变电磁探测试验。该矿井属于侏罗纪煤田，之前由于大量使用以掘代采的方式采煤，产生较多的小型采空区，且资料保留不完全，位置不明。矿区煤层上部有砂岩裂隙水，部分区段还有第四系松散层潜水，因此，这些采空区大部分为积水采空区，对该煤矿的安全生产造成了重要影响。

钻孔瞬变电磁探测测点范围为孔深60～128 m。对测点数据进行电阻率反演成像可以得到如图5 所示的沿钻孔方向的电阻率剖面图。由图5 可知，在钻孔深度100 m、钻孔径向20 m 位置有较为明显的低阻异常响应(蓝、绿色区域)，与数值模拟结果类似，从该成果图中无法清晰反映异常响应的边界，需要根据经验人为解释确定。

图5 实测数据电阻率剖面Fig.5 Resistivity profile of measured data

同样，采用本文1.3 节肘部法则对电阻率数据进行分析，可得到如图6 所示的ESS曲线图，由图6 可知，在聚类数目为3 时曲线出现肘点，因此，选择K=3。采用1.2 节中的K-Means 算法对电阻率数据进行聚类处理，可得到如图7 所示的异常响应边界成像结果。由图7 可知，在钻孔深度100 m、钻孔径向20 m 位置有明显的孤立低阻异常响应(蓝色区域)，异常响应边界非常清晰，异常响应形状表现为边长约20 m 的三角形。矿方后期设计了上仰钻孔对该异常进行勘查验证，在进尺98 m 附近出现掉钻和出水现象，最终推断该异常为上组煤遗留的小煤窑采空区。