許玉瑩， 陳 強， 常鎖亮， 張 生， 張興平

(1.太原理工大學(xué)礦業(yè)工程學(xué)院， 太原 030024； 2.山西山地物探技術(shù)有限公司， 晉中 030600)

隱伏陷落柱指未發(fā)育至目標煤層，僅發(fā)育于其底板下方巖層或灰?guī)r中(太灰或奧灰)的陷落柱[1]。由于其膠結(jié)程度低，裂隙發(fā)育，可能溝通奧灰?guī)r溶水、煤系含水層及地表水，易導(dǎo)致重大突水事故，而且具有很強的隱蔽性和突發(fā)性，嚴重影了煤炭的安全生產(chǎn)[2]。因此準確識別隱伏陷落柱對于煤礦安全高效生產(chǎn)至關(guān)重要。

由于奧灰頂界面通常為平行不整合界面，波阻抗差異較小，反射波能量弱，導(dǎo)致隱伏陷落柱的地球物理特征不明顯，反射波法識別隱伏陷落柱難度較大。針對這一問題，中外學(xué)者做了大量研究，總體來看，主要集中在以下3個方面：首先是地震采集和處理成像技術(shù)[3]，如采用寬方位采集和偏移距向量片(OVT)域地震數(shù)據(jù)處理技術(shù)等，OVT域全方位地震數(shù)據(jù)體分辨率高于常規(guī)數(shù)據(jù)體，可以提高隱伏陷落柱的識別精度[4]；其次是多屬性綜合識別技術(shù)，如采用相干、混沌、波形差異、相位、振幅等對陷落柱異常敏感度較高的屬性進行綜合識別[5-8]；三是對陷落柱成因地質(zhì)因素進行分級，采用模糊層次分析法對其賦以權(quán)重，形成模糊評判體系，從而達到預(yù)測陷落柱位置或者靶區(qū)的目的[9]。以上研究有效地提高了隱伏陷落柱的識別率，并能夠在平面中識別陷落柱位置，但在識別過程中多采用沿煤層屬性或?qū)娱g屬性，多解性強，需要對所識別的陷落柱異常逐一排除，效率較低；此外，相干、方差、瞬時相位、波形差異等屬性雖然可反映陷落柱的陷落邊界、冒落高度等垂向特征，但屬性剖面中陷落邊界不清，解釋位置擺動誤差較大，且無法判斷伴生斷裂發(fā)育情況[10]。

因此，現(xiàn)采用波形聚類法綜合層間屬性和沿層屬性對陷落柱異常識別，通過異常疊加，以實現(xiàn)對隱伏陷落柱的位置的自動識別，提高隱伏陷落柱的解釋效率；其次，將斷層概率體屬性引入隱伏陷落柱識別，并將其與導(dǎo)向濾波后的數(shù)據(jù)體進行融合，以精細刻畫陷落柱邊界、陷落高度及伴生斷裂發(fā)育特征。

1 方法原理

1.1 聚類分析

聚類分析的定義是將數(shù)據(jù)分組形成多個聚類，同一個聚類中，其對象之間具有較高的相似度，不同聚類中，對象差別較大[11]。其流程如下：①數(shù)據(jù)篩選，篩選出準確的、代表性的數(shù)據(jù)；②數(shù)據(jù)特征分析和提煉；③聚類算法選擇；④聚類效果測評；⑤分析解釋聚類結(jié)果。

研究采用分區(qū)聚類算法中的K-means聚類(K-means clustering)算法。具體計算流程為，將對象的數(shù)據(jù)集劃分為k個聚類，給定初始聚類的中心，計算各個聚類中心與其他n-k個數(shù)據(jù)點的距離，并將其歸并到距離最近的聚類，則一次迭代完成，利用式(1)重新計算每個聚類的中心：

(1)

式(1)中：Ci為聚類中心的值；xi為樣本值；n為樣本個數(shù)；k為聚類數(shù)；新的聚類中心的值與上一次迭代的聚類中心的值相等或方差小于某個閾值，則該算法收斂；反之則再次進行迭代分析，直到算法收斂[12-13]。

在三維地震勘探中，屬性平面可矢量化為離散的數(shù)值x、y、z，通過地質(zhì)地震分析給定需要識別的地質(zhì)異常閾值。利用聚類分析法計算時，需要優(yōu)化聚類數(shù)目k，以形成最佳聚類模型。優(yōu)化原則為當聚類結(jié)果達到屬性閾值，聚類數(shù)目達到最優(yōu)。

1.2 斷層概率體屬性

斷層概率體屬性(fault likehood attributes)是基于相干屬性(coherence attributes)的一種地震屬性，計算表達式[14-15]為

F=1-S8

(2)

式(2)中:S為相干屬性主要表征相鄰地震道反射波形的相似性，其范圍為0～1，值越大表示相似程度越小[16-17]，F(xiàn)為斷層概率體；F取值范圍為0～1,值越大表示該處為斷層的概率越大[18]。相較于原相干屬性值，加大了F低值與高值之間的反差，更有利于凸顯斷裂異常。

相干屬性定義[18-19]為

(3)

式(3)中:d為三維地震數(shù)據(jù)體；s表示對括號內(nèi)數(shù)據(jù)進行構(gòu)造導(dǎo)向平滑濾波；f表示對斷層沿其走向和傾向進行平滑濾波[18]。

在計算相干屬性前，首先對數(shù)據(jù)進行振幅增強處理，以增強相對較弱的反射振幅，基于增強后的數(shù)據(jù)體，計算出傾角體，在相干屬性的提取過程中進行傾角導(dǎo)向控制，從而提高其可靠性和準確性。

基于增強后的數(shù)據(jù)體，計算出應(yīng)用結(jié)構(gòu)向量G，結(jié)構(gòu)向量公式[17]為

(4)

式(4)中：v1、v2、v3為張量矩陣中的特征向量；λ1、λ2、λ3為張量矩陣中的特征值，該矩陣中特征值均為非負值，排序滿足λ1>λ2>λ3，故λ1為最大特征值，其對應(yīng)的特征向量v1表示分析單元內(nèi)法線方向所指示的最大變化率。

根據(jù)特征值和特征向量可計算出每一道地震數(shù)據(jù)對應(yīng)的傾角和方位角，形成傾角體。傾角和方位角計算公式[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方向的分量。

斷層概率體屬性在剖面中表現(xiàn)為連續(xù)的細線，該細線反映同相軸不連續(xù)點，在地質(zhì)上可解釋為斷裂邊界。斷層概率體與時間剖面疊合，可直觀反映陷落柱的陷落邊界、冒落高度和周邊伴生斷層發(fā)育情況，如圖1所示，圖1為地震剖面和斷層概率體的融合圖，其中圖1(b)中黑色細線即斷層概率體屬性，圖1中1號、2號異常為斷層，3號和4號異常均為陷落柱，其冒落高度特征明顯，且頂部均伴生小斷裂。

1.3 識別方法

基于以上原理，對隱伏陷落柱進行識別，首先進行構(gòu)造導(dǎo)向濾波，主要目的為保護和突出不連續(xù)點，以使得隱伏陷落柱的屬性特征更加明顯。在濾波的基礎(chǔ)上進行相干屬性和斷層概率體屬性計算。

在屬性體的基礎(chǔ)上，識別流程主要分為兩方面，一方面為平面識別，另一方面為剖面分析。平面識別中主要采用波形聚類方法對沿層和層間(奧灰與上組煤層之間)兩類相干屬性平面進行陷落柱識別，并將結(jié)果疊加，以達到快速篩選隱伏陷落柱的目的。聚類數(shù)目是波形聚類法的主要影響參數(shù)，不同的聚類數(shù)識別效果差異較大，需對比分析，選取合適的聚類數(shù)目。在剖面分析中，將斷層概率體屬性和地震濾波后地震數(shù)據(jù)體進行融合，對隱伏陷落柱進行垂向展布及伴生斷裂等三維空間特征刻畫。在斷層概率體屬性計算過程中，采用不同閾值所識別的陷落柱尺度不同，對于中等尺度陷落柱采用主斷裂(閾值為0.96)，研究中所識別的主要為奧灰隱伏陷落柱，故采用了保留小斷層及裂縫的概率體(閾值為1)。識別流程如圖2所示。

1為斷層；2為斷層；3為隱落柱；4為隱伏陷落柱圖1 地震原始剖面、斷層概率體屬性與地震數(shù)據(jù)體融合剖面對比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 實例應(yīng)用

2.1 研究區(qū)概況

研究區(qū)位于沁水盆地西部、太岳山隆起的東南部，主要煤層為山西組1號煤層和太原組9+10號煤層，埋深500～900 m。開采煤層為太原組9+10號煤層，奧灰頂界與9+10號煤層層間距約50 m。研究表明，該區(qū)存在眾多陷落柱，且發(fā)育特征明顯，呈一走向北東的條帶。其中部分陷落柱為奧灰隱伏陷落柱，因與所開采的太原組9+10號煤層間距很小，存在極高的導(dǎo)水隱患，故查明奧灰隱伏陷落柱的分布及冒落高度對研究區(qū)突水治理至關(guān)重要。

2.2 隱伏陷落柱異常識別

為突出陷落異常，在進行數(shù)據(jù)體屬性計算前，先對數(shù)據(jù)體進行構(gòu)造導(dǎo)向濾波，如圖3所示，圖3(a)為原始剖面，圖3(b)為構(gòu)造導(dǎo)向濾波后剖面，濾波后地震數(shù)據(jù)更加平滑，并且保存了信號的不連續(xù)性，使地震剖面不連續(xù)點更加容易識別，為識別陷落柱邊界提供了基礎(chǔ)。

圖3 原始地震時間剖面與濾波后時間剖面對比Fig.3 Comparison of original seismic section and filtered section

圖4 聚類數(shù)對比Fig.4 Comparison of cluster numbers

濾波后對數(shù)據(jù)體分別進行相干屬性和斷層概率體屬性計算，提取沿層屬性和層間屬性，采用波形聚類法對以上兩種屬性平面進行陷落柱識別，而聚類數(shù)影響聚類效果，圖4為不同聚類數(shù)識別結(jié)果對比，識別圖中相干值大于0.5的異常。圖4(a)為相干屬性平面，黑色箭頭所在位置異常小于0.5，圖4(b)～圖4(d)均識別出該異常，當聚類數(shù)為12時，識別結(jié)果可識別出閾值范圍內(nèi)異常。

在確定聚類數(shù)的前提下，對在相干屬性體上提取層屬性，圖5(a)為10號煤層沿層屬性平面，圖5(b)為奧灰頂界面與9+10號煤層層間屬性平面，采用聚類分析法分別對這兩個屬性平面中的陷落柱進行識別。聚類數(shù)選取優(yōu)化的聚類數(shù)12，如圖5(c)、圖5(d)所示，屬性平面中所有陷落柱均識別出，其中包括部分線性構(gòu)造如斷層。圖5(e)將沿層和層間識別出的陷落柱進行疊合，其中藍色為9+10號煤層沿層陷落柱識別結(jié)果，紅色為奧灰頂界面至9+10號煤層的層間陷落柱識別結(jié)果，奧灰隱伏陷落柱(未陷落至9+10號煤層的陷落柱)表現(xiàn)為紅色且未被藍色覆蓋的陷落柱，從圖5(e)中可確定8個疑似隱伏陷落柱，通過剖面對比，確定7個為奧灰隱伏陷落柱，如圖5(f)所示。隱伏陷落在平面中呈北東向條帶狀，與陷落柱整體發(fā)育規(guī)律一致。如采用編程進行自動識別可大大提高隱伏陷落柱的識別效率。

圖5 屬性平面圖及隱伏陷落柱圖像識別Fig.5 Attribute plan and image recognition of hidden collapse column

在確定隱伏陷落柱平面位置的基礎(chǔ)上，結(jié)合剖面分析其垂向發(fā)育情況，圖6為過陷落柱X3、X4、X5的任意時間剖面和屬性剖面對比。圖6中，T1為1號煤層反射波，T10為9+10號煤層反射波，To2為奧灰頂界面反射波。圖6(a)為濾波后的時間剖面，陷落柱剖面特征明顯，但發(fā)育形態(tài)和邊界不明顯；圖6(b)為相干屬性剖面，可直觀顯示波形不連續(xù)點，但陷落柱發(fā)育形態(tài)和冒落高度不明顯；圖6(c)為波形差異屬性剖面，可反映異常位置，但無法直觀判斷陷落柱發(fā)育形態(tài)和陷落邊界；圖6(d)為瞬時屬性剖面，從剖面中可大致判斷出陷落柱陷落邊界；圖6(e)為斷層概率體屬性剖面，剖面中3個隱伏陷落柱發(fā)育形態(tài)和冒落高度清晰，X3、X4、X5發(fā)育形態(tài)均為錐形，冒落高度均至9+10號煤層底板下部，其中X3和X5頂部伴生小斷裂，貫穿9+10號煤層，X4頂部無伴生斷裂。圖6(f)為斷層概率體屬性和時間剖面融合顯示，可直觀顯示陷落柱垂向邊界形態(tài)、陷落高度及伴生斷裂特征陷落發(fā)育形態(tài)，在剖面解釋過程中可提高解釋效率和精度。通過對比認為斷層概率體屬性與時間剖面融合，可在解釋過程中直觀刻畫陷落柱的發(fā)育形態(tài)，冒落高度及伴生斷裂發(fā)育情況。

圖6 地震原始剖面、多種屬性剖面及時間剖面與斷層概率體屬性疊合Fig.6 Seismic section and variety of different attribute section and co-rendering of seismic section and fault likehood attribute

結(jié)合平面和剖面共識別隱伏陷落柱7個，其中6個頂部存在伴生斷裂，存在導(dǎo)水性可能，為礦方安全、高效生產(chǎn)提供地質(zhì)依據(jù)。

由此可見，通過平面快速篩選出隱伏陷落柱，結(jié)合融合后的剖面分析隱伏陷落柱發(fā)育情況的隱伏陷落解釋手段，在隱伏陷落柱發(fā)育區(qū)可大大提高隱伏陷落柱解釋的精度和效率。

3 結(jié)論

通過理論分析與實例應(yīng)用得出如下結(jié)論。

(1)采用波形聚類法分別識別沿層和層間相干屬性平面中的陷落柱，將結(jié)果進行疊加可快速篩選出隱伏陷落柱平面分布規(guī)律，直接形成隱伏陷落柱平面圖，有利于陷落柱分布的宏觀認識，指導(dǎo)合理布置工作面。

(2)采用構(gòu)造導(dǎo)向濾波處理以及數(shù)據(jù)體融合技術(shù)將斷層概率體屬性和地震數(shù)據(jù)體融合，在剖面中陷落邊界及伴生斷裂顯示為纖細的線，可在剖面中精細刻畫隱伏陷落柱垂向邊界形態(tài)、陷落高度及伴生斷裂特征，同時可進行隱伏陷落柱的自動解釋，減少擺動誤差。

(3)采用平面快速篩選，結(jié)合剖面精確分析的方法識別隱伏陷落柱，可較大程度提高解釋精度和效率。這兩個過程均可實現(xiàn)自動識別，為隱伏陷落柱的識別提供了新的途徑，具有廣闊的應(yīng)用前景。