黃 斌，張宏兵，王 強(qiáng)，喬 驍

（河海大學(xué)地球科學(xué)與工程學(xué)院，南京 211100）

廣義S變換與短時傅里葉變換在地震時頻分析中的對比研究

黃 斌，張宏兵，王 強(qiáng)，喬 驍

（河海大學(xué)地球科學(xué)與工程學(xué)院，南京 211100）

短時傅里葉變換（STFT）和廣義S變換（GST）都被應(yīng)用到地震時頻分析中，但對兩者在信號分析過程中的特點(diǎn)和差異的研究相對較少。通過比較兩者的理論公式、窗口函數(shù)及地震信號的實(shí)際處理效果發(fā)現(xiàn)：短時傅里葉變換在地震信號分析過程中整個時頻域具有相同的分辨率，整體性較強(qiáng)，缺少時頻聚焦能力，不能對信號重點(diǎn)觀測區(qū)域有針對性的提高時頻分辨率；廣義S變換對地震高頻信號具有較高的時間分辨率，對低頻信號具有較高的頻率分辨率，且可以通過改變參數(shù)p的值對廣義S變換窗口函數(shù)的形態(tài)做出較大的調(diào)整，也可以改變λ的值實(shí)現(xiàn)窗口形態(tài)微調(diào)，通過對窗口函數(shù)的調(diào)整，廣義S變換可以對信號特定區(qū)域進(jìn)行時頻聚焦。

短時傅里葉變換；S變換；廣義S變換；時頻窗；地震記錄

0 引言

地震勘探數(shù)據(jù)是一種非平穩(wěn)的信號，其頻率值隨時間的變化較快，用傳統(tǒng)傅里葉變換進(jìn)行分析不能將時間值和頻率值對應(yīng)起來，分析效果不佳。采用短時傅里葉變換、S變換和廣義S變換等時頻分析方法則可以很清晰的反應(yīng)出地震信號在每個時頻點(diǎn)的幅值大小。S變換由R.G.Stockwel[1]于1996年提出，是在短時傅里葉變換和小波變換的基礎(chǔ)上發(fā)展而來的一種新型時頻分析方法。廣義S變換提出后即在地震勘探數(shù)據(jù)處理領(lǐng)域中得到廣泛的應(yīng)用，劉喜武[2]等人用廣義S變換來提取地震旋回，結(jié)論表明地震數(shù)據(jù)的振幅和頻率的旋回性和沉積旋回具有很好的對應(yīng)性；高靜懷[3]將廣義S變換用于研究薄互層地震響應(yīng)特征，研究得出恰當(dāng)選用基本小波，廣義S變換在薄互層識別中將起到重要作用；陳學(xué)華[4]、王長江[5]等人用廣義S變換實(shí)現(xiàn)了地震資料的高效時頻譜分解；姬戰(zhàn)懷[6]、齊春艷[7]研究和改進(jìn)了廣義S變換理論，使得廣義S變換類型增多，應(yīng)用范圍更廣。近年來在雷達(dá)信號[8]、圖像信息處理[9]等領(lǐng)域也開始使用廣義S變換且取得了較好的效果。

雖然短時傅里葉變換和廣義S變換已經(jīng)在地震信號時頻分析上得到應(yīng)用，但是關(guān)于兩者在信號分析上的優(yōu)缺點(diǎn)并沒有系統(tǒng)的研究，本文首先由短時傅里葉變換理論推導(dǎo)到廣義S變換理論，再研究兩者窗函數(shù)的差別，進(jìn)而研究廣義S變換中兩調(diào)節(jié)參數(shù)的具體作用，最后將兩者應(yīng)用到地震記錄時頻分析中，對比兩者在地震信號分析中的特點(diǎn)，為實(shí)際數(shù)據(jù)處理提供參考。

1 廣義S變換與短時傅里葉變換理論對比

1.1 從短時傅里葉（SFFT）到廣義S變換（GST）

傅里葉變換是將信號的時域和頻域聯(lián)系起來的方法，但是傅里葉變換主要用于平穩(wěn)信號的時頻轉(zhuǎn)換，是信號的一種整體變換，缺乏定位能力。為了解決這個問題，人們提出在時間域加窗的方法，即將時間域分割成不同的片段，對這些片段分別進(jìn)行傅里葉變換得到對應(yīng)的頻率值，這樣就實(shí)現(xiàn)了時間域和頻率域的局部對應(yīng)，將一維時間域轉(zhuǎn)換成二維時頻域，這也是短時傅里葉變換和廣義S變換的基本原理。

其中g(shù)(τ -t)為移動窗口函數(shù)，τ為時移因子，（1）式是由傅里葉變換公式加了一個窗口函數(shù)構(gòu)成。在時間域用g(t)窗函數(shù)去截取信號x(t)，可以認(rèn)為截取的信號為平穩(wěn)信號，對其做傅里葉變換得到t時刻附近的信號頻域值，移動窗口即可得到時間序列不同t值附近的頻域值，這些傅里葉變換的集合就是X(τ ，f)。X(τ ，f)是三維的復(fù)合函數(shù)，表示信號x(t)隨時間和頻率變化的相位和幅度值。

根據(jù)測不準(zhǔn)原理，當(dāng)窗函數(shù)為高斯類型函數(shù)時，短時傅里葉變換時頻窗口半徑最小。高斯窗可用下式表示：

觀察（2）式可知短時傅里葉變換時移因子τ一旦確定，窗函數(shù)也就確定，其窗寬和幅值不會隨著頻率值改變。為了讓窗函數(shù)和頻率值聯(lián)系起來，R.G. Stockwe[2]提出的S變換（ST）將（2）式中的σ與頻率f相關(guān)聯(lián)，令可得S變換的窗口函數(shù)為:

將（3）式帶入（1）式就可以得到S變換的公式：

廣義S變換又將S變換的窗函數(shù)加以改造，通過添加兩個調(diào)節(jié)因子λ和p來達(dá)到調(diào)節(jié)窗函數(shù)值的目的，其窗函數(shù)為：

將（5）式帶入（1）中即可得廣義S變換的公式：

通過上述公式推導(dǎo)可知S變換是在繼承短時傅里葉變換原理的基礎(chǔ)上，通過采用高斯類窗函數(shù)來提高時頻分辨率并且在窗函數(shù)中加入頻率變量，使得窗口的寬度和幅值能隨著頻率值做線性變化，具有與頻率相關(guān)的的分辨率，這就大大改善了短時傅里葉變換窗口在整個頻率軸上固定不變的缺點(diǎn)。廣義S變換是在S變換的窗口函數(shù)中引入調(diào)節(jié)參數(shù)，使得窗口函數(shù)隨頻率變化的規(guī)律可以調(diào)節(jié)，解決了S變換的時窗寬度和幅值隨頻率的變化值固定的這一缺點(diǎn)，使得窗口函數(shù)可以根據(jù)實(shí)際需要更加靈活的調(diào)控。

1.2 窗函數(shù)形態(tài)特征

圖1a、圖1b、圖1c分別為短時傅里葉變換高斯窗，S變換窗，廣義S變換窗。從圖中可以明顯看出：①短時傅里葉窗口函數(shù)窗寬和幅值沿頻率軸沒有發(fā)生改變。②S變換窗口函數(shù)幅值隨頻率的增大而增大，窗寬隨頻率的增大而降低，在高頻處有高幅值和窄窗口，可以判定其高頻處的時間分辨率很高，在低頻處具有低幅值和寬窗口，可以得出其低頻時間分辨率較低，根據(jù)測不準(zhǔn)原理又可得其高頻處的頻率分辨率很低而低頻處的頻率分辨率則較高。③廣義S變換窗函數(shù)形態(tài)大致和S變換窗函數(shù)相同，但是通過改變λ和p值可以發(fā)現(xiàn)其在高頻處的幅值較S變換窗幅值大，窗寬則變得更窄，表明廣義S變換通過參數(shù)調(diào)節(jié)使得其在高頻處的時間分辨率較S變換提高了。

為了更具體的了解廣義S變換窗口函數(shù)中調(diào)節(jié)參數(shù)λ和p值的調(diào)節(jié)特點(diǎn)，本文采用控制變量的方法，分別對單變量λ和p值做圖分析。

首先分析λ和p值的改變對窗函數(shù)的幅值和寬的影響。固定p=1，λ分別取0.8、0.9、1.0、1.1、1.2，如圖2a所示，觀察可知隨著λ的增大，窗的寬度逐漸減小，幅值逐漸增大，但是減小及增大的幅度不大。另外λ=1時為S變換的時窗曲線，可知若要在S變換的基礎(chǔ)上減小窗寬和加大幅值可以適當(dāng)增加λ的值，反之則適當(dāng)減小λ值。

固定λ=1，p值分別取0.8、0.9、1.0、1.1、1.2，如圖2b所示，觀察可知隨著p值的增大時窗幅值急劇增大，時窗寬度變窄的幅度也比較快。雖然效果和λ值變化時大致相同，但是p值的改變導(dǎo)致的時窗變化要比λ大的多。

通過對圖1c的觀察發(fā)現(xiàn)，改變λ、p的值對時窗的脊（時窗在各個頻率值剖面最大值的連線）的形態(tài)也會發(fā)生改變。為了研究其變化特性，分別取p= 0.8，λ=0.8、1.0、1.2；p=1.0，λ=0.8、1.0、1.2；p=1.2， λ=0.8、1.0、1.2作窗脊幅值隨頻率f的變化曲線如圖3所示。

圖1 時頻窗口Figure 1 Time-frequency window

圖2 λ、ρ值的改變對窗口幅值和寬度的影響Figure 2 Impacting on window amplitude and width from variation of λ and ρ values

圖3 窗脊的變化曲線Figure 3 Window ridge variation curve

觀察圖3可知p值固定時λ值的變化對曲線彎曲情況影響不大，對窗脊的幅值大小有一定影響，即λ越大，脊線的整體幅值越大。λ固定時，p值的變化對曲線的彎曲程度有較大影響，p＞1時，曲線向幅值減小方向凹陷，表明窗脊沿頻率減小的方向幅值減小的速率加快；p=1時，窗脊曲線隨頻率的變化大致呈線關(guān)系性，即幅值沿頻率軸變化的速率不變，表明S變換的時窗函數(shù)的幅值與頻率值呈線性變化關(guān)系；p＜1時，窗脊曲線大體往幅值增大方向凸起，表明窗脊沿頻率減小的方向幅值減小的速率變慢。p值的增大也會導(dǎo)致窗脊曲線的幅值整體增大，且增幅很大，比λ增大時帶來的增幅更大。

通過以上對p，λ的分析可以得出，p值的改變對時窗的窗寬、幅值、窗脊的變化具有決定性的作用，λ對窗的形態(tài)改變遠(yuǎn)不及p值，在實(shí)際應(yīng)用過程中可以起到對窗口微調(diào)的作用。

2 廣義S變換與短時傅里葉變換在地震記錄上的應(yīng)用對比

2.1 地震記錄模型的建立

在忽略各種能量衰減和干擾波的作用下，地震反射波記錄x（t）可以看成是地震子波b（t）和地震反射系數(shù)w（t）的褶積加上噪聲n（t）的結(jié)果，本文中的噪聲選用正態(tài)分布的隨機(jī)信號，其公式：

地震反射系數(shù)是上下地層波阻抗之間的函數(shù)，若上層波阻抗大于下層波阻抗，反射系數(shù)為負(fù)值，反之，反射系數(shù)為正值。本文采用的是Ricker子波模型[10]，模型公式：

式（8）中參數(shù)g＞0是子波的主頻率[11]，本文g取頻率為60 Hz，接近實(shí)際地震資料反射波主頻。圖4是本文模擬的地震信號圖形，圖4a為設(shè)置的10個反射系數(shù)，每一個反射系數(shù)值代表一個地層界面；圖4b為ricker子波模型，圖4c為其頻譜圖。可以看出ricker子波的主頻為60 Hz，最高頻約為150 Hz。圖4d為合成的單道地震記錄，采樣間隔為2 ms，采樣點(diǎn)數(shù)為300，時間從0～0.598 s。圖4e為其頻譜圖，可知合成的地震記錄在50～65 Hz以及80～100 Hz有較大的振幅值，最大頻率約為250 Hz。

2.2 地震記錄時頻分析結(jié)果對比

對圖4d所示的合成地震記錄利用matlab編程軟件分別對其做STFT、ST、GST時頻分析，結(jié)果如圖5所示。

圖4 合成單道地震記錄Figure 4 Synthetic single channel seismic records

圖5a為短時傅里葉時頻分析的結(jié)果，上文對其窗函數(shù)研究時提到，其窗函數(shù)形態(tài)與率值的大小無關(guān)，無法對特定的區(qū)域進(jìn)行聚焦，所以其分析效果比較差，沿時間軸只能大概看出反射層的位置，對于相距很近且反射系數(shù)較小的地層如圖4a所示的前三個反射系數(shù)值所代表的地層，在圖5a中基本無法識別出來。圖4e中的50～65 Hz及80～100 Hz兩個高振幅頻率區(qū)間對應(yīng)STFT時頻圖中0.25～0.3 s，30～100 Hz的能量團(tuán)，但是沿頻率軸方向觀察，能量團(tuán)粘聯(lián)在一起，無法區(qū)分兩個高振幅頻率區(qū)間，故其頻率分辨率也不足。

由于ST法是GST法取λ=1，p=1時的特例，故也屬于GST法，結(jié)果分析時統(tǒng)一看做GST法。圖5b為ST法時頻圖，有效信號在0.1～0.4 s，30～120 Hz內(nèi)，也不能得到特別好的時頻分辨率，但是相對圖5a來說已有所提升，尤其在低頻頻率分辨率上。有效信號區(qū)間外的一些淺色能量團(tuán)為噪聲信號，特點(diǎn)是噪聲能量團(tuán)比較分散不具連續(xù)性，但信噪比很低的地震信號，噪聲會對有效信號的識別帶來較大的影響。圖5c為調(diào)節(jié)參數(shù)λ=1.5，p=1.1的GST時頻圖，由上文窗函數(shù)的討論可知，這種條件下的窗函數(shù)較ST變換的窗函數(shù)在高頻處的幅值要大，窗寬要窄，故其在信號高頻處的時間分辨率比ST法要高。圖5c中的每個條帶狀能量團(tuán)代表一個地層，各條帶在時間尺度上間隔比較明顯，可見其時間分辨率比前面兩圖都要高。為了得到較好的頻率分辨率將GST窗的調(diào)節(jié)參數(shù)適當(dāng)降低，如圖5d所示的時頻圖采用λ=0.8，p=0.8的窗函數(shù)，可知其頻率分辨率較前面3個圖的要高。前面3個圖都無法分辨的50～65 Hz及80～100 Hz兩個高振幅頻率，在圖5d中0.25～0.3 s區(qū)域沿頻率軸方向被清楚的分成了2個高能量團(tuán)，表明在這種參數(shù)條件下的GST頻率分辨率高。圖5b、圖5c、圖5d聯(lián)合起來分析就可以比較細(xì)致的刻畫出地震信號在時頻域的能量分布，為后期濾波等處理提供技術(shù)支持。

圖5 單道地震記錄時頻分析Figure 5 Single channel seismic records time-frequency analysis

3 結(jié)論

①通過對STFT及GST窗函數(shù)的對比研究可得：STFT的窗函數(shù)一旦固定，其時頻分辨率在整個時頻域上相同且無法改變；GST的窗函數(shù)在高頻處具有幅值高，窗窄，在低頻處具有幅值低，窗寬等特點(diǎn)，所以其在信號高頻處具有高頻率聚焦性能，在低頻處具有高時間聚焦性能。通過對GST窗函數(shù)調(diào)節(jié)參數(shù)的研究表明：p值的改變對時窗的窗寬、幅值、窗脊的變化具有決定性的作用，λ對窗的形態(tài)改變遠(yuǎn)不及p值，在實(shí)際應(yīng)用過程中可以起到對窗口微調(diào)的作用。

②通過對合成地震記錄采用多種時頻方法對比分析表明：STFT法在地震數(shù)據(jù)頻譜分析過程中雖然具有一定整體性但缺乏時頻聚焦能力，既不能得到高的頻率分辨率也不能得到好的時間分辨率，在實(shí)際地震數(shù)據(jù)時頻分析過程中可以作為一種參照方法。GST法雖然不能在一張時頻圖中同時得到高的時間分辨率和頻率分辨率，但是可以改變參數(shù)值，分別得出地震信號的高頻率分辨率的時頻圖和高時間分辨率的時頻圖，聯(lián)合起來分析就可以較為清楚的了解地震信號在時頻域的能量分布，在對時頻分辨率要求較高的地震時頻分析過程中，可以利用其刻畫特定區(qū)域的時頻能量分布。

[1]Stockwell R G，Mansinha L，Lowe R P.Localization of the complex spectrum:the Stransform[J].IEEE on Signal Processing，1996，44（4）: 998-1001.

[2]劉喜武，劉洪，李幼銘，等.基于廣義S變換研究地震地層特征[J].地球物理學(xué)進(jìn)展，2006，21（2）:440-451.

[3]高靜懷，陳文超，李幼銘，等.廣義S變換與薄互層地震響應(yīng)分析[J].地球物理學(xué)報(bào)，2003，46（4）:526-532.

[4]陳學(xué)華，賀振華，黃德濟(jì).基于廣義S變換的地震資料高效時頻譜分解[J].石油地球物理勘探，2008，43（5）:530-534.

[5]王長江，楊培杰，羅紅梅，等.基于廣義S變換的時變分頻技術(shù)[J].石油物探，2013，52（5）:489-494.

[6]姬戰(zhàn)懷，嚴(yán)勝剛.關(guān)于“一種改進(jìn)廣義S變換”的討論[J].石油地球物理勘探，2015，50（6）:1224-1230.

[7]齊春艷，李彥鵬，彭繼新，等.一種改進(jìn)的廣義S變換[J].石油地球物理勘探，2010，45（2）:215-218.

[8]劉傳武，張智軍，畢篤彥.S變換在雷達(dá)目標(biāo)識別中的應(yīng)用[J].系統(tǒng)仿真學(xué)報(bào)，2008，20（12）:3290-3292.

[9]甄莉，彭真明.基于廣義S變換的圖像局部時頻分析[J].航空學(xué)報(bào)，2008，29（4）:1013-1019.

[10]梅金順，王潤秋.Ricker類地震子波[J].石油地球物理勘探，2012（s1）:8-14.

[11]韓海英，王志章，王宗俊，等.基于Ricker類地震子波的匹配追蹤[J].石油物探，2014，53（1）:17-25.

Contrastive Study of Generalized S Transform with Short-time Fourier Transform in Seismic Time-frequency Analysis

Huang Bin,Zhang Hongbing,Wang Qiang and Qiao Xiao
(School of Earth Science and Engineering,Hohai University,Nanjing,Jiangsu 211100)

The short-time Fourier transform(STFT)and generalized S transform(GST)are all applied to the seismic time-frequency analysis procedure,but without systematic study of difference and peculiarity between the two in signal analysis procedure.The paper through correlation of theoretical formulae,window functions and real processed effects of the two has found that the STFT during the seismic signal analysis procedure has identical resolution in whole time-frequency domain with better integrity,but lack of time-frequency focusing capability and can not enhance time-frequency resolution in key signal observation area.The GST has higher time resolution on seismic high frequency signals,higher frequency resolution on low frequency signals;and can make major adjustment on GST window functions through change of parameter p;also can change λ value to realize fine adjustment of window configuration; through adjustment of window functions,GST can carry out time-frequency focusing in signal specific areas.

STFT;S transform;GST;time-frequency window;seismic record

P631.4

A

10.3969/j.issn.1674-1803.2017.01.13

1674-1803(2017)01-0059-05

國家自然科學(xué)基金（41674113)（41374116）

黃斌（1990—），男，江西撫州人，河海大學(xué)，地質(zhì)資源與地質(zhì)工程專業(yè)，碩士，研究方向?yàn)榈卣鹂碧綌?shù)據(jù)處理。

2016-12-05

責(zé)任編輯：孫常長