基于廣義S變換的透射槽波埃里相識別

陳波，朱國維，武延輝，楊振強，周俊杰

(1.中國礦業(yè)大學(xué)(北京) 地球科學(xué)與測繪工程學(xué)院，北京 100083； 2.中國礦業(yè)大學(xué)(北京) 煤炭資源與安全開采國家重點實驗室，北京 100083)

0 引言

我國煤炭資源豐富，分布廣[1]，但煤層賦存地質(zhì)條件差，致災(zāi)因素較多，而地面地震勘探精度低，無法滿足煤礦企業(yè)對安全生產(chǎn)的需求，因此在煤礦工作面內(nèi)開展透射地震勘探為精細(xì)查明工作面內(nèi)的地質(zhì)情況提供了可能[2-3]。相對于圍巖，煤層速度低，在煤層內(nèi)激發(fā)地震波會產(chǎn)生沿煤層傳播的槽波[4-5]。由于在煤層內(nèi)傳播的槽波攜帶了煤層的地質(zhì)信息，因此通過對透射槽波的分析可以查明工作面內(nèi)的地質(zhì)異常。頻散是透射槽波的重要特征，埃里相分布在頻散曲線上速度最低能量強的位置，其特征明顯，通過對頻散曲線埃里相的分析可以獲得工作面內(nèi)煤層的速度及厚度[6-8]，從而查明工作面內(nèi)的地質(zhì)情況。因此，準(zhǔn)確地從透射槽波頻散曲線中識別埃里相顯得尤為關(guān)鍵。

目前，頻散曲線的提取方法主要有基于單道分析的多次濾波法、移動時窗法、時頻分析法及基于多道分析的F-K域變換和拉東變換等[6]。在透射槽波的分析中，需要對單道槽波提取頻散曲線，研究高精度的單道頻散曲線的提取方法對透射槽波分析具有重要意義。Feng和Zhang[9]提出了時窗可變的S變換提取透射槽波頻散曲線，取得了較好效果，但未討論如何準(zhǔn)確識別埃里相。筆者受Li等[10]的啟發(fā)，將廣義S變換引入到透射槽波的頻散曲線提取中，并對參數(shù)的取值做出分析，通過選取合理參數(shù)突出埃里相特征，確保準(zhǔn)確地識別出頻散曲線中的埃里相。

1 廣義S變換

Stockwell等[11]提出了S變換，其定義式如下：

(1)

其逆變換為：

(2)

Li等[10]構(gòu)造了一種形式更為簡潔的兩參數(shù)廣義S變換并展示了其良好的應(yīng)用效果，定義式如下[10]：

圖1 廣義S變換窗函數(shù)時間半高寬隨頻率變化Fig.1 Temporal FWHM of generalized S-transform window versus frequency

(4)

(5)

式中：f1、f2為頻率，tFWHM1、tFWHM2為窗函數(shù)時間半高寬。當(dāng)頻率為f1時，窗函數(shù)的時間半高寬為tFWHM1；當(dāng)頻率為f2時，窗函數(shù)的時間半高寬為tFWHM2。附錄A給出了詳細(xì)推導(dǎo)，可以通過給定兩個頻率點的時間半高寬快速確定A和B。

為檢驗廣義S變換在提高信號時頻分辨率方面的性能，利用合成信號進(jìn)行了初步測試，圖2a為頻率呈階梯變化的非平穩(wěn)信號[22]，其中前部分為低頻并包含一個高頻信號，中間部分為中頻，后部分為高頻，信號的具體表達(dá)式如式(6)所示：

(6)

圖2b、2c分別為圖2a所示信號的S變換和廣義S變換結(jié)果，從變換的結(jié)果可以看到，廣義S變換信號的頻率分辨率顯著提高，在時頻域中，階梯狀信號的能量聚集性更好，即頻率分辨率提高，由海森堡不確定性原理，其時間分辨率降低，圖2c中持續(xù)時間較短的高頻信號能量相對發(fā)散，表明了其時間分辨率降低。圖2d展示了圖2b和圖2c虛線處的線譜，相比S變換，廣義S變換信號能量更集中，頻率分辨率得到大幅提高。圖2e為50 Hz時的S變換和廣義S變換窗函數(shù)比較，廣義S變換窗函數(shù)呈現(xiàn)出“矮”、“寬”特征，圖2f為窗函數(shù)時間半高寬對比圖，廣義S變換窗函數(shù)整體比S變換寬，因此，其頻率分辨率高于S變換。任取兩點，由式(4)和式(5)即可確定A、B值，本例中A=0.4239、B=0.471。

a—頻率階梯變化的非平穩(wěn)信號;b—圖2a的S變換結(jié)果;c—圖2a的廣義S變換結(jié)果;d—圖2b、圖2c虛線處的線譜;e—50 Hz時的S變換和廣義S變換窗函數(shù)振幅隨時間變化;f—S變換和廣義S變換窗函數(shù)時間半高寬對比a—non-stationary signal with frequency step change;b—S-transform results of Fig. 2a;c—generalized S-transform results of Fig. 2a;d—line spectra at dotted lines in Fig. 2b and Fig. 2c;e—S-transform and generalized S-transform window amplitude versus time at 50 Hz;f—comparison of the temporal FWHM between S-transform and generalized S-transform圖2 頻率階梯變化的非平穩(wěn)信號的S變換和廣義S變換(A=0.4239，B=0.471)對比分析Fig.2 Comparison of S-transform and generalized S-transform (A=0.4239，B=0.471) for non-stationary signal with frequency step change

2 合成數(shù)據(jù)應(yīng)用

利用文獻(xiàn)[23]給出的Love型槽波計算方法合成槽波，計算參數(shù)如表1所示。

表1 合成Love型槽波物性參數(shù)

合成的Love型槽波如圖3a所示，道間距為5 m，共41道。取第30道進(jìn)行群速度頻散曲線提取，提取過程簡述如下：

首先將地震道變換到時頻域S(t,f)，然后用炮檢距x除以t，即得v=x/t，最后對S(v,f)重新排序并插值得到群速度頻散曲線圖。

圖3b為第30道Love型槽波記錄，圖3c為圖3b數(shù)據(jù)S變換提取的群速度頻散曲線圖，其中星號線為Love型槽波群速度理論頻散曲線[24-25](下同)，由圖可見，提取的頻散曲線與理論頻散曲線一致性較好，說明了提取的頻散曲線的正確性。在本例中，廣義S變換的目的是適度提高提取的頻散曲線頻率分辨率，突出埃里相特征。圖3d為圖3b數(shù)據(jù)廣義S變換結(jié)果，從圖可知，頻散曲線埃里相(圖3c、3d黑色圓圈所示區(qū)域，下同)能量聚集性明顯比S變換更好，更利于準(zhǔn)確快速拾取埃里相。本例中，廣義S變換A、B值與圖2一致，廣義S變換的目的都是適度提高頻率分辨率。

a—合成Love型槽波記錄;b—圖3a中第30道槽波記錄;c—圖3b數(shù)據(jù)S變換提取的頻散曲線;d—圖3b數(shù)據(jù)廣義S變換提取的頻散曲線a—synthetic Love channel waves;b—the 30th channel wave in Fig. 3a;c—dispersion curve of Fig. 3b extracted by S-transform;d—dispersion curve of Fig. 3b extracted by generalized S-transform圖3 合成Love型槽波數(shù)據(jù)S變換和廣義S變換(A=0.4239，B=0.471)提取頻散曲線對比分析Fig.3 Comparison of the dispersion curve extracted by S-transform and generalized S-transform (A=0.4239，B=0.471) for synthetic Love channel waves

為進(jìn)一步說明廣義S變換A、B取值對頻散曲線埃里相特征的影響，與圖3相反，我們?nèi)∫唤M值A(chǔ)=2.1195、B=2.355，使廣義S變換窗函數(shù)時間半高寬整體比S變換窄，廣義S變換窗函數(shù)呈現(xiàn)出“高”、“窄”特征，如圖4a所示，圖4b為窗函數(shù)時間半高寬對比，廣義S變換窗函數(shù)整體比S變換窄，其時間分辨率高于S變換，由海森堡不確定性原理，頻率分辨率低于S變換。利用此廣義S變換對圖3b數(shù)據(jù)進(jìn)行處理，提取的群速度頻散曲線如圖4c所示(埃里相所在區(qū)域由黑色圓圈所示)，可以看到，相比S變換提取的群速度頻散曲線圖3c，圖4c中埃里相所在區(qū)域黑色圓圈所示完全無法識別，因此，通過圖3、圖4的分析可知，利用廣義S變換適度提高頻率分辨率更有利于埃里相的識別。

a—50 Hz時的S變換和廣義S變換窗函數(shù)振幅隨時間變化;b—S變換和廣義S變換窗函數(shù)時間半高寬對比;c—圖3b所示數(shù)據(jù)廣義S變換提取的頻散曲線a—S-transform and generalized S-transform window amplitude versus time at 50 Hz;b—comparison of temporal FWHM between S-transform and generalized S-transform;c—dispersion curve of Fig. 3b extracted by generalized S-transform圖4 合成Love槽波數(shù)據(jù)廣義S變換(A=2.1195，B=2.355)提取頻散曲線Fig.4 Dispersion curve extracted by generalized S-transform (A=2.1195，B=2.355) for synthetic Love channel waves

為測試廣義S變換在含噪情況下對數(shù)據(jù)的處理能力，在圖3a所示的數(shù)據(jù)中加入隨機(jī)噪聲，得到如圖5a所示數(shù)據(jù)。圖5a中第30道槽波記錄如圖5b所示，對圖5b所示數(shù)據(jù)進(jìn)行S變換和廣義S變換分析，廣義S變換A、B值與圖3取值一致，圖5c為S變換提取的頻散曲線，圖5d為廣義S變換提取的頻散曲線，兩種方法提取的頻散曲線與理論頻散曲線一致性均較好，但廣義S變換提取的頻散曲線埃里相更易于識別，進(jìn)一步說明廣義S變換在數(shù)據(jù)含噪情況下也能突出埃里相特征。

3 實際資料應(yīng)用

圖6a為某工作面第12炮第10道透射槽波記錄，圖6b、6c分別為圖6a數(shù)據(jù)S變換和廣義S變換(A=0.4239，B=0.471)提取的頻散曲線，從兩圖對比看到，相比圖6b，圖6c埃里相能量更集中，群速度頻散曲線埃里相特征更明顯，可以迅速確定埃里相位于圖6c中黑色圓圈所示區(qū)域，而從圖6b所示頻散曲線圖中識別埃里相則相對困難。因此，廣義S變換提取的頻散曲線更有利于準(zhǔn)確拾取透射槽波埃里相，使解釋人員能較快速地識別埃里相頻率及速度，從而為下一步工作面層析成像研究奠定堅實基礎(chǔ)。

4 結(jié)論

1)討論了廣義S變換中參數(shù)的取值與時頻分辨率的關(guān)系，并推導(dǎo)了其參數(shù)與窗函數(shù)時間半高寬的具體表達(dá)式，為定量表征時頻分辨率奠定基礎(chǔ)。

2)廣義S變換的數(shù)值試驗表明，通過合理選擇參數(shù)，廣義S變換能夠有效提高時間或頻率分辨率，根據(jù)海森堡不確定性原理，在提高頻率分辨率的同時，時間分辨率會有所下降，反之亦然，因此，需根據(jù)實際情況合理選擇參數(shù)，提高時頻分辨率。

3)將廣義S變換引入到透射槽波頻散曲線提取中，研究認(rèn)為廣義S變換能有效改善頻散曲線的時頻分辨率，通過選取合適的參數(shù)能進(jìn)一步有效突出頻散曲線中埃里相特征，使解釋人員能準(zhǔn)確快捷地識別透射槽波埃里相。

a—圖3a加入隨機(jī)噪聲后的Love型槽波記錄;b—圖5a中第30道槽波記錄;c—圖5b數(shù)據(jù)S變換提取的頻散曲線;d—圖5b數(shù)據(jù)廣義S變換提取的頻散曲線a—synthetic Love channel waves with adding random noise to Fig. 3a;b—the 30th channel wave in Fig. 5a;c—dispersion curve of Fig. 5b extracted by S-transform;d—dispersion curve of Fig. 5b extracted by generalized S-transform圖5 合成Love槽波數(shù)據(jù)中加入隨機(jī)噪聲后S變換和廣義S變換(A=0.4239，B=0.471)提取頻散曲線對比分析Fig.5 Comparison of the dispersion curve extracted by S-transform and generalized S-transform (A=0.4239，B=0.471) for synthetic Love channel waves with adding random noise

a—某工作面第12炮第10道透射槽波記錄;b—圖6a數(shù)據(jù)S變換提取的頻散曲線;c—圖6a數(shù)據(jù)廣義S變換提取的頻散曲線a—the 10th transmitted channel wave of the 12th shot in a coal mining face;b—dispersion curve of Fig. 6a extracted by S-transform;c—dispersion curve of Fig. 6a extracted by generalized S-transform圖6 實際槽波數(shù)據(jù)的S變換和廣義S變換(A=0.4239，B=0.471)提取頻散曲線對比分析Fig.6 Comparison of the dispersion curve extracted by S-transform and generalized S-transform(A=0.4239，B=0.471) for real data