基于聚类算法的钻孔瞬变电磁视电阻率立体成像方法

范 涛,李鸿泰,刘 磊，赵 睿,李博凡,郭建磊，李宇腾

(1. 中煤科工集团西安研究院有限公司，陕西 西安 710077; 2. 长安大学 地球物理场多参数综合模拟实验室(中国地球物理学会重点实验室)，陕西 西安 710054； 3. 四川省地质矿产勘查开发局物探队，四川 成都 610072)

0 引 言

随着中国煤炭开发逐渐向深部发展，复杂地质构造和隐蔽致灾因素对煤矿安全生产影响加重，尤其是在煤矿巷道掘进过程中，面临复杂未知的地质条件，对地质保障技术的先进性与精准性提出更高的要求。2018年颁布的《煤矿防治水细则》中就对地球物理超前探测技术提出了新的要求，其中明确提出探测点20 m范围内不得有积水和金属物体，然而目前主流的超前探测技术均不符合该项规定，急需研发新的超前探测手段[1]。

研究表明，利用掘进工作面的钻孔进行瞬变电磁探测工作，可以在掘进前开展远距离、高精度的隐伏水害超前预报。钻孔瞬变电磁法基本施工装置如图1所示。图1中，发射线圈和接收线圈共同集成在孔内探管中，通过孔口控制一次场发射与关断，激发探管附近的低阻异常体产生二次场，再接收二次场的三分量信号，最后解释得到钻孔周围的电性分布。该方法的优势在于：①发射线圈在目标体附近激发，能在防爆限制下最大限度激发煤层附近的目标体；②孔中接收线圈既避开了巷道中的铁磁性干扰，还能最大限度减少目标体二次场因距离带来的能量损耗；③通过对水平分量的处理定位孔旁地质异常体所在方位，可利用单钻孔探测资料对异常体进行立体解释。

图1 钻孔瞬变电磁法基本施工装置示意图Fig.1 Schematic View of Basic Construction Device forBorehole Transient Electromagnetic Method

早在20世纪70年代，国外学者就已开展了利用钻孔瞬变电磁三分量探测信号解释孔旁地质异常体的研究：Woods通过比例模型实验研究，总结出一套解释板体模型不同参数的特征关系曲线[2]；Macnae等阐述了导电背景中的井中响应符号变化现象和特征[3]；Buselli等模拟了导电覆盖层下多个目标体的信号响应[4]；Kozhevnikov等研究了钻孔套管对孔中瞬变电磁响应的影响[5]。中国学者自20世纪80年代引入地-井瞬变电磁装备后也开展了相关研究[6-28]：胡平等开展了基于自由空间球体和板体的地-井瞬变电磁响应理论计算，对国外已报道的结果进行了补充[6]；张杰推导了矩形回线在空间任意点处产生的一次场表达式，提出三分量数据矢量交汇技术[7]；杨毅等提出基于导电薄板等效涡流的异常反演方法[8]；孟庆鑫等通过大地介质影响下的正演模拟，确定了围岩背景场对总响应的影响结果[9]；徐正玉等采用时域有限差分法模拟研究了接触带不同埋深位置和接触面两侧不同电阻率对信号的影响[10-12]；杨海燕等研究了覆盖层影响下板状体异常响应规律[13]；武军杰等定义了电性源地-井瞬变电磁全域视电阻率[14]；陈卫营等对电性源在地下激发的6个电磁场分量的扩散、分布特性和探测能力进行了分析研究[15]。在隧/巷道内工作的钻孔瞬变电磁法近几年才被提出，相关研究资料较少。国外只有Vella曾将地-井瞬变电磁发射线圈移到金属矿巷道中来探测含金块状黄铁矿体[16]。国内王世睿研究了隧道10 m以内浅孔中的瞬变电磁响应特征，提出利用移动掌子面上发射线圈位置来定性判断异常体方位的施工技术[17]；孙怀凤等通过物理模拟实验证明了孔中瞬变电磁信号可用于判断隧道掘进工作面前方是否存在异常构造[18]；陈丁等通过在全空间一维背景上增加三维异常体的积分方程数值模拟研究了煤矿巷道垂直孔中瞬变电磁特性[19]；范涛研究了钻孔瞬变电磁的叠加超前探测方法和径向探测数据的二维拟地震反演方法[20-23]。

从前人研究成果可以看出，水平分量形态组合和幅值差异对孔旁地质异常体位置敏感，结合垂直分量视电阻率成像结果可对异常体进行立体解释。但是，根据水平分量异常形态组合确定异常体所在象限需人工进行识别和判断，效率较低，尤其当测点较多时，人工逐点逐道识别异常曲线形态更是不现实的工作。因此，钻孔瞬变电磁法当前的立体解释还处于定性水平，水平分量异常曲线形态自动识别是阻碍钻孔瞬变电磁法实现立体成像的关键因素，本文通过机器学习中的K-means聚类算法实现对大数据水平分量异常曲线形态的自动分类。

1 方法原理

1.1 垂直分量数据处理解释方法

本文使用的钻孔瞬变电磁法本质上仍属于中心回线装置类型，其垂直(z)分量实测数据曲线形态与矿井瞬变电磁探测数据曲线形态基本相同(图2)，仅是发射线圈尺寸与匝数不同导致电感影响较大，因此，数据处理方法可参考矿井瞬变电磁，采用范涛等提出的预处理技术[29]对电感影响进行校正，对校正后的数据可采用晚期视电阻率计算和层厚累加法进行处理解释。

B为磁感应强度图2 矿井与钻孔瞬变电磁探测数据曲线对比Fig.2 Comparison of Transient Electromagnetic DetectionData Curves Between Mine and Borehole

晚期视电阻率计算公式为

(1)

式中：ρ为晚期视电阻率；μ0为真空磁导率，取值为4π×10-7H·m-1；ST为发射线圈面积；SR为接收线圈面积；t为测道时间，V(t)/I为归一化感应电动势，I为电流；C为视电阻率修正系数。

深度计算公式为

(2)

(3)

式中：Hi为第i道深度；K为深度修正系数；hi为第i道双程传播深度；ρi为第i道晚期视电阻率；ti为第i道采样时间。

1.2 水平分量与异常体所在象限的规律

设计如图3所示模型，规定水平x与y分量正方向之间区域为第一象限，逆时针依次定义为第二、三、四象限，在钻孔深度方向50 m处，依次放置16个规模为20 m×20 m×6 m的板状异常体，异常体中心点距离钻孔15 m，模型其他参数见表1。绘制16个模型的水平分量异常响应多测道图(图4～11)。

表1 模型参数

图(a)中数字为异常体编号图3 模型示意图Fig.3 Schematic Views of Model

图4～11分别为4个象限水平分量异常响应。从图4～11可以看出，以水平钻孔为参考系，所有水平分量异常响应均呈“S”或“反S”形态，且当某一水平分量与该水平分量坐标轴夹角为0°时，该水平分量的响应幅值达到最小。两组水平分量形态组合与异常体所在象限之间的关系如图12所示。

1.3 异常体工具面角计算方法

钻孔瞬变电磁法径向探测时，异常体引起的二次场是矢量场。根据水平涡流场的空间分布特征，假设异常体为板状体，在钻孔中观测到两个磁场水平分量Vx、Vy的矢量和Vxy，其方向一定处于由钻孔指向异常体的等效涡流中心上，那么只需要求出Vxy的方向，就知道异常体中心的具体方位[6]。

如图13所示，设Vxy与x轴夹角为θ，则

(4)

其中，Vx、Vy均为已知值，求反正弦即得到θ的表达式为

(5)

最后，根据异常体所在象限，可分别求出对应的工具面角α：①异常体在第一象限，α=θ；②异常体在第二象限，α=π-θ；③异常体在第三象限，α=π+θ；④异常体在第四象限，α=2π-θ。

1.4 水平分量异常曲线类型自动识别方法

欲对水平分量异常场曲线进行自动分类，应以水平分量每一个测点为插值窗口中心，采用Hermit插值求得所有窗口对应的异常场。插值窗口大小的选择可根据垂直分量中主要异常区统计平均值确定。之后对提取出的异常数据进行正规化，将曲线横坐标范围统一为从编号1开始的相同点号区间，然后对区间正规化后的所有异常数据进行特征提取，提取数据中极值对应的正规化点号，最后以极大值点号为x轴，极小值点号为y轴，形成特征值分布图(图14)。

图14是一组数值模拟数据的水平分量异常场特征值分布。从图14可以看出其具有明显的二分类特性，与水平分量异常场的“S”或“反S”形态存在显著相关性。基于该特性，无需提前进行标签样本的监督训练，可直接选用无监督机器学习算法对数据进行分类。

无监督机器学习算法常常被应用在数据挖掘领域，用于在大量无标签数据中发现规律。它的训练数据是无标签的，训练目标是能对观察值进行分类或区分等。常用的无监督学习算法主要有主成分分析方法、等距映射方法、局部线性嵌入方法、拉普拉斯特征映射方法、黑塞局部线性嵌入方法、局部切空间排列方法和最常用的聚类算法。

图4 第一象限水平x分量异常响应Fig.4 Horizontal x Component Anomaly Responses ofthe First Quadrant

图5 第一象限水平y分量异常响应Fig.5 Horizontal y Component Anomaly Responses ofthe First Quadrant

图6 第二象限水平x分量异常响应Fig.6 Horizontal x Component Anomaly Responses ofthe Second Quadrant

图7 第二象限水平y分量异常响应Fig.7 Horizontal y Component Anomaly Responses ofthe Second Quadrant

图8 第三象限水平x分量异常响应Fig.8 Horizontal x Component Anomaly Responses ofthe Third Quadrant

图9 第三象限水平y分量异常响应Fig.9 Horizontal y Component Anomaly Responses ofthe Third Quadrant

图10 第四象限水平x分量异常响应Fig.10 Horizontal x Component Anomaly Responses ofthe Fourth Quadrant

图11 第四象限水平y分量异常响应Fig.11 Horizontal y Component Anomaly Responses ofthe Fourth Quadrant

图12 异常体位于不同象限时水平分量响应形态Fig.12 Horizontal Component Response Forms ofAnomalous Body in Different Quadrants

图13 异常偏转角示意图Fig.13 Schematic View of Anomalous Deflection Angle

图14 特征值分布Fig.14 Distribution of Eigenvalues

聚类算法是指将一堆没有标签的数据自动划分成几类的方法，这个方法要保证同一类数据有相似的特征。考虑到水平分量异常场特征值本身就有明显的二分类特性，无需采用复杂的聚类算法，因此，本文选择应用广泛、速度快、鲁棒性强的K-means聚类算法。K-means聚类算法使用最大期望算法求解的高斯混合模型在正态分布的协方差为单位矩阵，且是在隐变量的后验分布为一组狄拉克δ函数时所得到的特例，其假设相同类别中数据之间的距离应该都很近，即数据之间的相似度与它们之间的欧式距离成反比。

本文需要将N个水平分量纯异常数据{xn}聚为2类，令经过聚类之后每个数据所属的类别为{tn}，而这2个聚类的中心为{μm}，则损失函数(L)计算公式为

(6)

实际计算时，先随机设置2个质心把所有数据粗略分成2个初始类，计算所有数据与质心的欧式距离，再根据平均值重新计算质心和类别，对以上过程反复迭代，直至达到终止条件。终止条件可设置为簇中心点变化率(RLabel,n)。其表达式为

(7)

将分类好的数据类别分别同“S”与“反S”形态进行对应，根据图12就可以准确确定钻孔瞬变电磁法观测数据中任意测点任意测道反映的视电阻率信息所在的象限，再按照异常体工具面角计算方法，将每一测点每一测道对应的视电阻率信息视为异常体带入计算，即可获取对应的工具面角。

1.5 视电阻率立体成像算法

应用1.1节方法可由钻孔瞬变电磁法探测的垂直分量计算并绘制以钻孔深度为横坐标、以钻孔径向探测距离为纵坐标的视电阻率剖面成像图。在已知图中每一测点每一测道视电阻率对应的工具面角信息时，假设钻孔为直钻孔，可通过三角函数关系将其投影在xOy平面上，将剖面图中每一测点的一维视电阻率信息曲线转换为二维视电阻率信息平面，多个测点的二维成像结果组合即可实现钻孔瞬变电磁视电阻率立体成像。

单个测点视电阻率信息一维坐标扩展至二维坐标的计算公式为

(8)

(9)

式中：j是测点数；i是测道数；x是钻孔瞬变电磁坐标系中视电阻率信息对应的x方向坐标；y是钻孔瞬变电磁坐标系中视电阻率信息对应的y方向坐标；r是钻孔探测半径。

2 数值模拟

为验证基于水平分量异常场特征聚类的钻孔瞬变电磁立体成像方法探测效果，设计如图15所示的三维模型，使用孙怀凤等提出的三维时域有限差分正演程序[30]进行数值模拟。在钻孔深度方向40 m、第一象限45°处放置1个规模为5 m×5 m×5 m

图15 数值模型示意图Fig.15 Schematic Views of Numerical Model

的低阻异常体，异常体中心点距离钻孔15 m，模型其他参数如表1所示。

对26～54 m的测点垂直分量数据采用晚期视电阻率计算和层厚累加法进行视电阻率成像，视电阻率修正系数取1.58，深度修正系数取1.00，可以得到如图16所示的沿钻孔方向的视电阻率剖面。从图16可以较为清晰地看到，在钻孔深度40 m、钻孔径向15 m位置有较为明显的低阻异常响应，但从该成果图中无法反映异常的空间方位。

图16 视电阻率剖面Fig.16 Profile of Apparent Resistivity

采用1.3、1.4节方法对模型水平x、y分量数据进行处理，可得到如图17所示的单测点视电阻率展开平面图。从图17可以看出，在孔深40 m的平面图中，xOy平面第一象限有明显的低阻异常响应，与坐标系原点O(钻孔)之间的距离为15 m，形状规模与模型参数基本一致，而在孔深30 m和50 m的平面图中则没有明显低阻异常显示。

图17 单测点视电阻率展开平面图Fig.17 Expanded Plan Views for Apparent Resistivity ofSingle Measuring Point

将每一测点的x、y坐标与钻孔孔深z坐标组合，通过Voxler软件进行立体成像如图18所示。图18中，由蓝绿色表示的低阻异常体与模型设置参数一致，说明基于水平分量异常场的聚类算法对钻孔径向视电阻率立体成像有效，准确性较高。

图18 视电阻率立体成像Fig.18 Stereo Imaging of Apparent Resistivity

3 结 语

(1)以水平钻孔钻进方向为z轴正方向，以孔口所在平面右向为x轴正方向，下向为y轴正方向，钻孔瞬变电磁法所有水平分量异常响应均呈“S”或“反S”形态，且当某一水平分量与该水平分量坐标轴夹角为0°时，该水平分量的响应幅值达到最小，因此，通过水平x、y分量异常形态组合可判定异常体所在象限。

(2)由两组水平分量的幅值基于三角函数关系可计算得到异常体中心在异常所在象限内的偏转角度，据此可得出异常体工具面角。

(3)将测量数据每一测点每一测道对应的视电阻率视为一个独立的异常体，采用K-means聚类算法对相应的两组水平分量异常曲线中的特征值进行二分类，可以自动确定任意测点任意测道视电阻率的分布象限，再由水平分量异常场幅值计算出工具面角，即可结合垂直分量成像结果实现钻孔径向视电阻率的立体成像。

(4)采用三维时域有限差分数值模拟数据对基于水平分量异常场的聚类算法进行了检验，成像结果与模型吻合度较高，说明该方法准确、有效地提高了钻孔瞬变电磁法的解释水平，可以推广应用至煤矿实际生产中，为掘进工作面隐伏水害精准超前探测提供技术保障。

谨以此文庆祝母校七十周年华诞，祝愿母校光辉历程更辉煌，人才辈出代代强，桃李满天扬四海，硕果累累振中华！2001年9月到2008年6月，我在母校度过了七年的求学时光。刚入学的时候，罗马广场还是食堂，八号学生公寓还没竣工，运动操场全是泥地，军训时经常一阵风刮过就卷起一片沙尘；之后短短几年时间，母校就迎来了高速发展，渭水校区基本建成，交通科技大厦、地学科技大厦也都拔地而起，呈现一片欣欣向荣的景象。回想起曾经那个“土土”的校园和无忧无虑、一起学习进步的同学，总让人想穿越回去再次感受那份善良与纯真、真诚与热情！