许玉莹， 陈 强， 常锁亮， 张 生， 张兴平

(1.太原理工大学矿业工程学院， 太原 030024； 2.山西山地物探技术有限公司， 晋中 030600)

隐伏陷落柱指未发育至目标煤层，仅发育于其底板下方岩层或灰岩中(太灰或奥灰)的陷落柱[1]。由于其胶结程度低，裂隙发育，可能沟通奥灰岩溶水、煤系含水层及地表水，易导致重大突水事故，而且具有很强的隐蔽性和突发性，严重影了煤炭的安全生产[2]。因此准确识别隐伏陷落柱对于煤矿安全高效生产至关重要。

由于奥灰顶界面通常为平行不整合界面，波阻抗差异较小，反射波能量弱，导致隐伏陷落柱的地球物理特征不明显，反射波法识别隐伏陷落柱难度较大。针对这一问题，中外学者做了大量研究，总体来看，主要集中在以下3个方面：首先是地震采集和处理成像技术[3]，如采用宽方位采集和偏移距向量片(OVT)域地震数据处理技术等，OVT域全方位地震数据体分辨率高于常规数据体，可以提高隐伏陷落柱的识别精度[4]；其次是多属性综合识别技术，如采用相干、混沌、波形差异、相位、振幅等对陷落柱异常敏感度较高的属性进行综合识别[5-8]；三是对陷落柱成因地质因素进行分级，采用模糊层次分析法对其赋以权重，形成模糊评判体系，从而达到预测陷落柱位置或者靶区的目的[9]。以上研究有效地提高了隐伏陷落柱的识别率，并能够在平面中识别陷落柱位置，但在识别过程中多采用沿煤层属性或层间属性，多解性强，需要对所识别的陷落柱异常逐一排除，效率较低；此外，相干、方差、瞬时相位、波形差异等属性虽然可反映陷落柱的陷落边界、冒落高度等垂向特征，但属性剖面中陷落边界不清，解释位置摆动误差较大，且无法判断伴生断裂发育情况[10]。

因此，现采用波形聚类法综合层间属性和沿层属性对陷落柱异常识别，通过异常叠加，以实现对隐伏陷落柱的位置的自动识别，提高隐伏陷落柱的解释效率；其次，将断层概率体属性引入隐伏陷落柱识别，并将其与导向滤波后的数据体进行融合，以精细刻画陷落柱边界、陷落高度及伴生断裂发育特征。

1 方法原理

1.1 聚类分析

聚类分析的定义是将数据分组形成多个聚类，同一个聚类中，其对象之间具有较高的相似度，不同聚类中，对象差别较大[11]。其流程如下：①数据筛选，筛选出准确的、代表性的数据；②数据特征分析和提炼；③聚类算法选择；④聚类效果测评；⑤分析解释聚类结果。

研究采用分区聚类算法中的K-means聚类(K-means clustering)算法。具体计算流程为，将对象的数据集划分为k个聚类，给定初始聚类的中心，计算各个聚类中心与其他n-k个数据点的距离，并将其归并到距离最近的聚类，则一次迭代完成，利用式(1)重新计算每个聚类的中心：

(1)

式(1)中：Ci为聚类中心的值；xi为样本值；n为样本个数；k为聚类数；新的聚类中心的值与上一次迭代的聚类中心的值相等或方差小于某个阈值，则该算法收敛；反之则再次进行迭代分析，直到算法收敛[12-13]。

在三维地震勘探中，属性平面可矢量化为离散的数值x、y、z，通过地质地震分析给定需要识别的地质异常阈值。利用聚类分析法计算时，需要优化聚类数目k，以形成最佳聚类模型。优化原则为当聚类结果达到属性阈值，聚类数目达到最优。

1.2 断层概率体属性

断层概率体属性(fault likehood attributes)是基于相干属性(coherence attributes)的一种地震属性，计算表达式[14-15]为

F=1-S8

(2)

式(2)中:S为相干属性主要表征相邻地震道反射波形的相似性，其范围为0～1，值越大表示相似程度越小[16-17]，F为断层概率体；F取值范围为0～1,值越大表示该处为断层的概率越大[18]。相较于原相干属性值，加大了F低值与高值之间的反差，更有利于凸显断裂异常。

相干属性定义[18-19]为

(3)

式(3)中:d为三维地震数据体；s表示对括号内数据进行构造导向平滑滤波；f表示对断层沿其走向和倾向进行平滑滤波[18]。

在计算相干属性前，首先对数据进行振幅增强处理，以增强相对较弱的反射振幅，基于增强后的数据体，计算出倾角体，在相干属性的提取过程中进行倾角导向控制，从而提高其可靠性和准确性。

基于增强后的数据体，计算出应用结构向量G，结构向量公式[17]为

(4)

式(4)中：v1、v2、v3为张量矩阵中的特征向量；λ1、λ2、λ3为张量矩阵中的特征值，该矩阵中特征值均为非负值，排序满足λ1>λ2>λ3，故λ1为最大特征值，其对应的特征向量v1表示分析单元内法线方向所指示的最大变化率。

根据特征值和特征向量可计算出每一道地震数据对应的倾角和方位角，形成倾角体。倾角和方位角计算公式[20]为

(5)

(6)

式中:P(x,y,t)为倾角；Q(x,y,t)为方位角；v1x(x,y,t)、v1y(x,y,t)、v1t(x,y,t)分别为v1(x,y,t)在x、y、t方向的分量。

断层概率体属性在剖面中表现为连续的细线，该细线反映同相轴不连续点，在地质上可解释为断裂边界。断层概率体与时间剖面叠合，可直观反映陷落柱的陷落边界、冒落高度和周边伴生断层发育情况，如图1所示，图1为地震剖面和断层概率体的融合图，其中图1(b)中黑色细线即断层概率体属性，图1中1号、2号异常为断层，3号和4号异常均为陷落柱，其冒落高度特征明显，且顶部均伴生小断裂。

1.3 识别方法

基于以上原理，对隐伏陷落柱进行识别，首先进行构造导向滤波，主要目的为保护和突出不连续点，以使得隐伏陷落柱的属性特征更加明显。在滤波的基础上进行相干属性和断层概率体属性计算。

在属性体的基础上，识别流程主要分为两方面，一方面为平面识别，另一方面为剖面分析。平面识别中主要采用波形聚类方法对沿层和层间(奥灰与上组煤层之间)两类相干属性平面进行陷落柱识别，并将结果叠加，以达到快速筛选隐伏陷落柱的目的。聚类数目是波形聚类法的主要影响参数，不同的聚类数识别效果差异较大，需对比分析，选取合适的聚类数目。在剖面分析中，将断层概率体属性和地震滤波后地震数据体进行融合，对隐伏陷落柱进行垂向展布及伴生断裂等三维空间特征刻画。在断层概率体属性计算过程中，采用不同阈值所识别的陷落柱尺度不同，对于中等尺度陷落柱采用主断裂(阈值为0.96)，研究中所识别的主要为奥灰隐伏陷落柱，故采用了保留小断层及裂缝的概率体(阈值为1)。识别流程如图2所示。

1为断层；2为断层；3为隐落柱；4为隐伏陷落柱图1 地震原始剖面、断层概率体属性与地震数据体融合剖面对比Fig.1 Comparison of original seismic section and co-rendering of seismic section and fault likehood attribute

图2 隐伏陷落柱识别流程图Fig.2 Identification of hidden collapse column

2 实例应用

2.1 研究区概况

研究区位于沁水盆地西部、太岳山隆起的东南部，主要煤层为山西组1号煤层和太原组9+10号煤层，埋深500～900 m。开采煤层为太原组9+10号煤层，奥灰顶界与9+10号煤层层间距约50 m。研究表明，该区存在众多陷落柱，且发育特征明显，呈一走向北东的条带。其中部分陷落柱为奥灰隐伏陷落柱，因与所开采的太原组9+10号煤层间距很小，存在极高的导水隐患，故查明奥灰隐伏陷落柱的分布及冒落高度对研究区突水治理至关重要。

2.2 隐伏陷落柱异常识别

为突出陷落异常，在进行数据体属性计算前，先对数据体进行构造导向滤波，如图3所示，图3(a)为原始剖面，图3(b)为构造导向滤波后剖面，滤波后地震数据更加平滑，并且保存了信号的不连续性，使地震剖面不连续点更加容易识别，为识别陷落柱边界提供了基础。

图3 原始地震时间剖面与滤波后时间剖面对比Fig.3 Comparison of original seismic section and filtered section

图4 聚类数对比Fig.4 Comparison of cluster numbers

滤波后对数据体分别进行相干属性和断层概率体属性计算，提取沿层属性和层间属性，采用波形聚类法对以上两种属性平面进行陷落柱识别，而聚类数影响聚类效果，图4为不同聚类数识别结果对比，识别图中相干值大于0.5的异常。图4(a)为相干属性平面，黑色箭头所在位置异常小于0.5，图4(b)～图4(d)均识别出该异常，当聚类数为12时，识别结果可识别出阈值范围内异常。

在确定聚类数的前提下，对在相干属性体上提取层属性，图5(a)为10号煤层沿层属性平面，图5(b)为奥灰顶界面与9+10号煤层层间属性平面，采用聚类分析法分别对这两个属性平面中的陷落柱进行识别。聚类数选取优化的聚类数12，如图5(c)、图5(d)所示，属性平面中所有陷落柱均识别出，其中包括部分线性构造如断层。图5(e)将沿层和层间识别出的陷落柱进行叠合，其中蓝色为9+10号煤层沿层陷落柱识别结果，红色为奥灰顶界面至9+10号煤层的层间陷落柱识别结果，奥灰隐伏陷落柱(未陷落至9+10号煤层的陷落柱)表现为红色且未被蓝色覆盖的陷落柱，从图5(e)中可确定8个疑似隐伏陷落柱，通过剖面对比，确定7个为奥灰隐伏陷落柱，如图5(f)所示。隐伏陷落在平面中呈北东向条带状，与陷落柱整体发育规律一致。如采用编程进行自动识别可大大提高隐伏陷落柱的识别效率。

图5 属性平面图及隐伏陷落柱图像识别Fig.5 Attribute plan and image recognition of hidden collapse column

在确定隐伏陷落柱平面位置的基础上，结合剖面分析其垂向发育情况，图6为过陷落柱X3、X4、X5的任意时间剖面和属性剖面对比。图6中，T1为1号煤层反射波，T10为9+10号煤层反射波，To2为奥灰顶界面反射波。图6(a)为滤波后的时间剖面，陷落柱剖面特征明显，但发育形态和边界不明显；图6(b)为相干属性剖面，可直观显示波形不连续点，但陷落柱发育形态和冒落高度不明显；图6(c)为波形差异属性剖面，可反映异常位置，但无法直观判断陷落柱发育形态和陷落边界；图6(d)为瞬时属性剖面，从剖面中可大致判断出陷落柱陷落边界；图6(e)为断层概率体属性剖面，剖面中3个隐伏陷落柱发育形态和冒落高度清晰，X3、X4、X5发育形态均为锥形，冒落高度均至9+10号煤层底板下部，其中X3和X5顶部伴生小断裂，贯穿9+10号煤层，X4顶部无伴生断裂。图6(f)为断层概率体属性和时间剖面融合显示，可直观显示陷落柱垂向边界形态、陷落高度及伴生断裂特征陷落发育形态，在剖面解释过程中可提高解释效率和精度。通过对比认为断层概率体属性与时间剖面融合，可在解释过程中直观刻画陷落柱的发育形态，冒落高度及伴生断裂发育情况。

图6 地震原始剖面、多种属性剖面及时间剖面与断层概率体属性叠合Fig.6 Seismic section and variety of different attribute section and co-rendering of seismic section and fault likehood attribute

结合平面和剖面共识别隐伏陷落柱7个，其中6个顶部存在伴生断裂，存在导水性可能，为矿方安全、高效生产提供地质依据。

由此可见，通过平面快速筛选出隐伏陷落柱，结合融合后的剖面分析隐伏陷落柱发育情况的隐伏陷落解释手段，在隐伏陷落柱发育区可大大提高隐伏陷落柱解释的精度和效率。

3 结论

通过理论分析与实例应用得出如下结论。

(1)采用波形聚类法分别识别沿层和层间相干属性平面中的陷落柱，将结果进行叠加可快速筛选出隐伏陷落柱平面分布规律，直接形成隐伏陷落柱平面图，有利于陷落柱分布的宏观认识，指导合理布置工作面。

(2)采用构造导向滤波处理以及数据体融合技术将断层概率体属性和地震数据体融合，在剖面中陷落边界及伴生断裂显示为纤细的线，可在剖面中精细刻画隐伏陷落柱垂向边界形态、陷落高度及伴生断裂特征，同时可进行隐伏陷落柱的自动解释，减少摆动误差。

(3)采用平面快速筛选，结合剖面精确分析的方法识别隐伏陷落柱，可较大程度提高解释精度和效率。这两个过程均可实现自动识别，为隐伏陷落柱的识别提供了新的途径，具有广阔的应用前景。