基于Morlet小波的時頻分析方法在地震槽波精準(zhǔn)反演煤厚中的應(yīng)用

郭昌放，楊 真，武 祥，陳一鼎，顧 薌

(1.中國礦業(yè)大學(xué)人工智能研究院，江蘇省徐州市，221116；2.中國礦業(yè)大學(xué)礦業(yè)工程學(xué)院，江蘇省徐州市，221116；3.昆明煤炭科學(xué)研究所，云南省昆明市，650051；4.方圓標(biāo)志認證集團，北京市海淀區(qū)，100048)

近些年來，我國袁亮院士[1]以及澳大利亞聯(lián)邦科學(xué)與工業(yè)研究組織[2]等國內(nèi)外眾多學(xué)者與機構(gòu)曾先后提出，工作面透明化作為煤礦智能精準(zhǔn)開采實現(xiàn)的基礎(chǔ)和前提，對推動煤炭資源的安全高效開采具有重要意義。而煤礦工作面透明化建設(shè)的核心是如何精準(zhǔn)識別工作面內(nèi)部的地質(zhì)異常環(huán)境[3-4]。同時，在回采過程中，煤厚的突然變化以及工作面內(nèi)部隱伏的地質(zhì)構(gòu)造經(jīng)常會造成工作人員傷亡和巨大的經(jīng)濟損失[5-6]。因此，在工作面回采之前，煤厚及地質(zhì)構(gòu)造的精準(zhǔn)預(yù)測對煤炭安全、智能、精準(zhǔn)開采具有重要意義。

目前工作面內(nèi)煤厚及地質(zhì)構(gòu)造探測方法已在煤礦得到廣泛應(yīng)用。其中，常見的鉆探方法[7]對于單個鉆孔對應(yīng)的煤厚數(shù)據(jù)是準(zhǔn)確的，然而受成本和時間限制，探測區(qū)域鉆孔數(shù)量有限，存在探測盲區(qū)；三維地震勘探[8]是在地表全方位布設(shè)炮點和檢波點，進行異常面積觀測的一種地震勘探方法，然而由于探測距離較大而且易受上覆巖層和地形條件的影響，探測精度較低；無線電波坑透技術(shù)[9]是煤厚探測的常用物探方法，然而對于內(nèi)部起伏較大的工作面探測效果較差，而且探測結(jié)果只能用于定性判斷。

槽波地震勘探技術(shù)由于其探測距離大、精度高、抗電干擾能力強以及波形易于識別[10-11]等優(yōu)點，成為目前國內(nèi)外煤礦工作面內(nèi)煤厚及地質(zhì)構(gòu)造探測和反演的主要方法[12-13]。其中，地震槽波走時的提取對煤厚及地質(zhì)構(gòu)造反演結(jié)果的精準(zhǔn)性至關(guān)重要。筆者基于Morlet小波的高分辨性，利用時頻分析方法將地震槽波在煤層中的傳播以能量流動的形式進行分析，求取地震波信號的時間-能量譜，并拾取地震槽波的走時?；贛orlet小波的時頻分析方法在邯鄲礦區(qū)郭二莊煤礦22204工作面槽波地震勘探過程中，預(yù)測結(jié)果與回采揭露結(jié)果基本保持一致，得到了良好的驗證。

1 槽波走時分析方法

地震波在煤層中激發(fā)產(chǎn)生之后，經(jīng)過頂?shù)捉缑娴亩啻稳瓷湟约跋嗷ジ缮嫘纬刹鄄?。槽波在煤層中傳播的過程中如果遇到斷層、陷落柱以及煤厚變化等異常，槽波傳播的速度、相位等參數(shù)就會發(fā)生變化。因此，槽波地震勘探方法就是利用槽波對異常體的高敏感性從而對煤層中的地質(zhì)異常區(qū)域進行預(yù)測和解釋[14]。

頻散是地震槽波最基本的特性，槽波傳播的速度主要取決于頻率的大小。因此，槽波同時具有相速度和群速度，然而拾取單一相速度走時比較困難，而且誤差較大。因此，通常選取槽波群速度的走時進行層析成像反演計算。在槽波走時分析過程中，通常是根據(jù)煤厚的探測范圍，利用某一固定頻率在頻散曲線上拾取槽波的群速度和對應(yīng)的走時[10,13]。然而，頻率的選擇并不具有唯一性，而且，只適用某一固定的煤厚探測范圍。因此，對于較大范圍的煤厚變化，必須選取不同的頻率重復(fù)拾取走時，工作量較大。另一方面，頻散曲線的不規(guī)整性，對槽波走時的精準(zhǔn)拾取帶來不利影響。為了解決上述問題，提高槽波走時拾取的準(zhǔn)確性，采用時頻分析方法獲得透射槽波傳播的時間-能量譜，以槽波能量包絡(luò)中心點對應(yīng)的時間作為槽波的走時，為地質(zhì)異常的精準(zhǔn)反演奠定了基礎(chǔ)。

在地震勘探中，由于受到多次反射、地層吸收以及頻散的影響，地震信號往往具有非平穩(wěn)特征[15]。常規(guī)傅立葉變換無法表征任一時刻的頻率成分，而時頻分析作為分析非平穩(wěn)時變信號的主要工具[16]，能夠描述頻譜隨時間的變化情況，并可以使得地震波信號的能量密度分布在時域和頻域中同時顯示。其中，小波變換作為一種時間尺度域分析技術(shù)，廣泛應(yīng)用于地球物理勘探數(shù)據(jù)分析中[17]。通過變換能夠?qū)r間(空間)和頻率進行局部化分析，通過伸縮平移運算對信號逐步進行多尺度細化，最終達到高頻處時間細分、低頻處頻率細分，因此能自動適應(yīng)時頻信號分析的要求，從而可聚焦到信號的任意細節(jié)，是一種多分辨率的分析方法。

小波分析得益于小波基函數(shù)的完備性、自相似性和多分辨性,母小波選擇不當(dāng)，應(yīng)用效果將會受影響[18]。由于Morlet小波在時域和頻域上具有較好的分辨率，它是一種復(fù)雜的單頻正弦信號調(diào)制高斯波:

式中：ψσ(t)——Morlet小波隨時間變化的波形振幅曲線；

t——時間；

Cσ——歸一化常數(shù)；

σ——無量綱頻率；

κσ——小波容許準(zhǔn)則。

Morlet小波ψσ(t)對應(yīng)的傅立葉變換表達式為：

(4)

w——頻率。

(ωψ-σ)2-1=(ωψ2-1)e-σωψ

(5)

其中，參數(shù)σ決定了時域和頻域之間的關(guān)系。一般而言，由于波形頻率隨時間變化較快，σ的取值往往大于5。因此κσ的取值也接近于零，而且中心頻率ωψ?σ。方程(1)也可以簡化為：

(6)

由于σ為無量綱頻率，歸一化常數(shù)Cσ=1，則方程(4)可以簡化為：

(7)

為了保證Morlet變換在每一個縮放因子s下具有可比性，需要對小波函數(shù)進行歸一化處理，使其具有相同的單位能量：

ωk——第k個頻率分量；

N——頻率索引。

2 工程應(yīng)用

2.1 地質(zhì)概況

郭二莊煤礦22204工作面煤層結(jié)構(gòu)簡單，煤層厚度0.20～2.20 m，平均1.48 m，平均傾角23°，兩巷落差近60 m，由于兩巷落差較大，因此選擇槽波地震勘探方法。在工作面兩巷掘進過程中共揭露了10條斷層，位于探測區(qū)域范圍內(nèi)的斷層共有9條，均為正斷層。22204工作面探測區(qū)域如圖1所示。

圖1 22204工作面探測區(qū)域

2.2 觀測系統(tǒng)和數(shù)據(jù)采集

本次勘探以炸藥爆破作為震源，通過單分量檢波器采集地震波，炮孔深度為3 m，藥量為300 g，發(fā)射點和接收點盡量布置在煤層中部。探測區(qū)域測點布置如圖2所示。所有的發(fā)射點(S1～S31)布置在運輸巷，間距10 m，要求炮孔直徑為42 mm，盡量保持與煤層平行。所有的接收點(R1～R69)都固定在錨桿上，間距10 m。數(shù)據(jù)采集時，炮點從S31到S1逐一爆破，產(chǎn)生地震波并被檢波器接收。在數(shù)據(jù)采集過程中，由于檢波器R1、R2、R3和R69出現(xiàn)故障，因此共采集到2015道震波數(shù)據(jù)，圖2中給出了3個炮點對應(yīng)接收點的射線路徑。

圖2 探測觀測系統(tǒng)圖

在22204工作面透射槽波勘探中，采用以Morlet小波為母小波的小波變化時頻分析方法來拾取槽波走時。震源點為S5時，在接收點R7處采集到的地震波數(shù)據(jù)時頻分析結(jié)果如圖3所示。從圖中能夠清楚看到槽波信號同時在時域和頻域中的能量分布情況，并以槽波能量包絡(luò)中心點對應(yīng)的時間作為槽波的走時，對采集的2015道地震波數(shù)據(jù)進行時頻分析，共拾取了1750個有效槽波走時數(shù)據(jù)，為進一步的層析成像反演奠定了基礎(chǔ)。

圖3 震源為S5時接收點R7采集到的地震波數(shù)據(jù)時頻分析結(jié)果

2.3 層析成像反演

槽波層析成像(CT)技術(shù)已廣泛應(yīng)用于工作面內(nèi)部的透射投影，利用槽波通過非均勻介質(zhì)傳播的時間，對其內(nèi)部慢度分布的結(jié)果進行重構(gòu)[19]。運用CT技術(shù)進行反演計算的常用算法有：最小平方法(LSQR)、代數(shù)重建法(ART)以及瞬時迭代法(SIRT)，由于SIRT對于稀疏、不規(guī)則以及低信噪比的探測數(shù)據(jù)具有較好的適用性，而且比較穩(wěn)定[20]，因此選用SIRT算法用于層析成像反演。

(10)

對探測區(qū)域進行網(wǎng)格離散化，如圖4所示，劃分為N個單元，每個單元內(nèi)的慢度值分別為S1,S2,S3,…,SN。公式(10)可轉(zhuǎn)變?yōu)椋?/p>

(11)

式中：tk——第k條射線走時；

dki——第k條射線經(jīng)過第i面元的長度；

Si——第i面元的慢度。

若射線的條數(shù)為M，則所有射線路徑走時方程可轉(zhuǎn)化為：

T=D·S

(12)

式中：T——M維時間矩陣；

S——N維慢度矩陣；

D——(M*N)矩陣算子。

利用瞬時迭代重構(gòu)法(SIRT)求解式(11)，從而可以得到慢度反演結(jié)果。

基于采集的地震槽波走時，通過SIRT算法進行層析成像迭代反演，由于網(wǎng)格的尺寸和數(shù)目都會影響層析成像的結(jié)果，因此在反演之前，需要對這些參數(shù)進行多次測試。最終網(wǎng)格尺寸定義為10.3 m × 10.3 m，共23行、61列。通過多次迭代計算，最終得到了探測區(qū)域槽波傳播慢度分布結(jié)果，如圖5所示。

圖4 探測區(qū)域網(wǎng)格離散化分析

圖5 探測區(qū)域慢度分布結(jié)果

2.4 慢度與煤厚之間的關(guān)系

槽波傳播的慢速與傳播的介質(zhì)往往具有較好的對應(yīng)關(guān)系，而且槽波在煤層中傳播的慢度比在巖層中傳播的值要高。圖5顯示了探測區(qū)域內(nèi)部的槽波傳播慢度分布結(jié)果，其中藍色低慢度區(qū)域主要分布在400～680 m之間，初步判斷該區(qū)域的煤厚值比周圍的要低。然而，通過慢度分布只能對探測區(qū)域內(nèi)煤厚的變化情況進行定性判斷。

為了更清楚地了解探測區(qū)域煤厚的定量分布，針對22204工作面回風(fēng)巷200～500 m以及運輸巷400～750 m區(qū)段揭露的煤厚結(jié)果(0.2～2.2 m)和對應(yīng)的慢度值進行數(shù)值擬合，由于試驗數(shù)據(jù)近似于線性分布，因此選擇一階多項式擬合煤厚-慢度關(guān)系，如圖6所示。

圖6 煤厚-慢度擬合結(jié)果

3 煤厚分布結(jié)果

利用一階多項式擬合的煤厚-慢度線性關(guān)系(圖6)，將慢度分布結(jié)果(圖5)轉(zhuǎn)換成煤厚分布結(jié)果，如圖7所示。22204工作面探測區(qū)域在700～900 m之間，由于受巷道布置限制，射線覆蓋比較稀疏，會對結(jié)果的準(zhǔn)確性造成影響。

由圖7中可看出，探測區(qū)域煤厚中間薄、兩端厚。其中兩端平均煤厚都在1.5 m以上，而中部煤厚變化在0.5～1.5 m之間，由于該區(qū)域共揭露了6條正斷層，判斷由于受斷層影響，頂?shù)装鍘r層在張力的作用下延伸至煤層導(dǎo)致該區(qū)域煤層變薄。而且該薄煤區(qū)域延伸范圍較大，對工作面回采影響較大。

22204工作面沿著回采方向生產(chǎn)過程中，在開采到706 m的位置時，由于煤厚太薄，無法繼續(xù)回采，最終在470 m位置新開切眼，自右向左繼續(xù)完成回采工作。因此經(jīng)過生產(chǎn)驗證該探測結(jié)果與實際相符，能夠為工作面合理布置提供參考依據(jù)。

圖7 探測區(qū)域煤厚分布結(jié)果

該探測區(qū)域回采完成之后，通過揭露的煤厚與預(yù)測煤厚之間的相關(guān)系數(shù)r進行對比驗證，相關(guān)系數(shù)r定義為：

(12)

式中：D——揭露的煤厚結(jié)果，m；

P——預(yù)測的煤厚結(jié)果，m；

Cov(D,P)——揭露煤厚與預(yù)測煤厚之間的協(xié)方差；

Var[D]——揭露煤厚的方差；

Var[P]——預(yù)測煤厚的方差。

22204工作面回采后揭露的結(jié)果與預(yù)測結(jié)果進行對比分析，其結(jié)果如圖8所示。揭露的煤厚與預(yù)測的煤厚之間的相關(guān)系數(shù)r為0.88，具有較好的吻合效果。

4 結(jié)論

基于Morlet小波的時頻分析方法在郭二莊煤礦22204工作面煤厚探測工程應(yīng)用中的預(yù)測結(jié)果與實際揭露結(jié)果基本保持一致，研究成果能夠為工作面內(nèi)地震槽波走時的精準(zhǔn)提取提供借鑒意義，具體結(jié)論如下所述。

圖8 預(yù)測煤厚與揭露煤厚對比結(jié)果

(1)選擇基于Morlet小波的小波變換時頻分析方法，進行槽波走時拾取時，不受固定頻率和煤厚探測范圍的限制，能夠獲得更加客觀精準(zhǔn)的槽波走時，能夠有效提高煤厚及地質(zhì)異常反演的精準(zhǔn)性。

(2)由于22204工作面傾角較大，探測效果仍然比較理想，因此該方法對于大傾角工作面具有較好的適用性。

(3)由于地球物理方法往往存在多解性，而且探測結(jié)果的準(zhǔn)確性還依賴于地球物理條件，因此應(yīng)該結(jié)合鉆探等其他的勘探手段，進一步提高預(yù)測結(jié)果的準(zhǔn)確性。