阮清青，張會星，王 昊，李凱瑞

(中國海洋大學 海底科學與探測技術(shù)教育部重點實驗室，山東 青島 266100)

?

修正S變換與常規(guī)時頻分析方法的對比

阮清青，張會星，王 昊，李凱瑞

(中國海洋大學 海底科學與探測技術(shù)教育部重點實驗室，山東 青島 266100)

地震信號往往是非線性、非平穩(wěn)信號，基于平穩(wěn)信號理論的常規(guī)傅里葉變換方法不能刻畫任一時刻的頻率成分，而時頻分析技術(shù)能同時展示信號在時間域和頻率域的局部化特征。修正S變換作為一種較新的時頻分析方法，是針對S變換窗函數(shù)相對固定、時頻分辨率不能調(diào)節(jié)的問題提出的。通過對修正S變換與常規(guī)時頻分析方法如短時傅里葉變換(STFT)、小波變換(WT)和S變換的對比，分析了修正S變換在時頻分析方法中的作用。根據(jù)合成信號以及實測地震記錄的時頻分析可知：修正S變換較常規(guī)的時頻分析方法具有更好的時頻分辨率和能量聚集性，更有利于對非平穩(wěn)信號的處理。

時頻分析；修正S變換；S變換；小波變換；短時傅里葉變換

0 引言

在地震勘探中，地震波的傳播會受到不同的巖性、物性和含流體地層的影響，導致地震信號具有非線性、非平穩(wěn)特征，即實際地震資料的頻譜成分及信號的各種統(tǒng)計特性是隨時間而變化的。傳統(tǒng)的傅氏變換是一種整體變換，只能了解信號在時間域或者頻率域的全局特征，無法起到時頻定位作用。時頻分析( Time-Frequency Analysis)作為一種新興的信號處理方法，它將信號從時間域變換到時間-頻率域，不同頻率分量的時間關(guān)聯(lián)特性能夠在時頻平面上有效地表示出來，而且信號在任意時刻的能量都聚集在瞬時頻率附近，其逆變換則可以得到等效的時間信號[1]。在時頻分析研究和發(fā)展的短短幾十年時間里，它已經(jīng)在地震勘探數(shù)據(jù)分析的多個方面得到了重要的應(yīng)用[2-5]。

常用的時頻分析方法有短時傅里葉變換、S變換、小波變換、Wigner-Ville分布等[6-9]。短時傅里葉變換(STFT)是通過滑動時窗來計算其頻譜，因而它的時間分辨率和頻率分辨率受Heisenberg測不準原理約束[10]。小波變換繼承了STFT局部化分析信號的思想，但選取的是窗口形狀可變的窗函數(shù)，在分析時變非平穩(wěn)信號時較STFT具有更理想的效果。1996年，美國地球物理學家Stockwell[7]綜合了短時傅里葉變換(STFT)和小波變換(WT)的特點，首次提出了一種加時窗傅里葉變換方法——S變換，以較高的時頻分辨率而備受關(guān)注[11-13]。S 變換是對短時傅里葉變換和連續(xù)小波變換的繼承與發(fā)展，它綜合了兩種變換的優(yōu)點，卻彌補了它們的不足。然而，S變換具有時間和頻率的分辨率不能同時提高的問題，而且標準S變換高斯窗口隨頻率的變化固定，導致其時頻譜的時間分辨率較低[14]。

由于S變換的基本小波是固定不變的，使得S變換的應(yīng)用受到限制，許多學者通過改進基本小波或窗函數(shù)得到了改進的S變換，并取得了良好的應(yīng)用效果[15]。Mansinha 等[16]用(f/r)代替f，得到調(diào)諧的高斯函數(shù)，允許使用者自定S變換在時頻面上時間和頻率的分辨率。Pinnegar 等[17]提出用非對稱的雙曲窗代替高斯窗，用于地震波的P波首波時間的判定。高靜懷等[18]用4個待定參數(shù)的調(diào)幅簡諧波來代替S變換中的基本小波，用于薄層及薄互層的地震探測。陳學華等[19]引入調(diào)節(jié)參數(shù)λ和p對S變換的窗函數(shù)進行改造，用于信號的提取和抑噪。Duan Li和John Castagna[20]在2013年提出一種修正的S變換新方法，對標準S變換中的高斯窗函數(shù)進行改進，利用一個線性頻率方程代替高斯窗的頻率，調(diào)節(jié)時窗寬度隨頻率呈反比變化的速度，提高了S變換在具體應(yīng)用中的實用性和靈活性。本文將修正S變換與常規(guī)的時頻分析方法進行對比分析，通過對幾種合成信號的頻譜分解，研究修正S變換在時間分辨率、頻率分辨率和能量聚集性方面的特點。

1 時頻分析算法

1.1 短時傅里葉變換(STFT)

Dennis Gabor[6]針對傅里葉變換在時頻分析的不足，提出了短時傅里葉變換。其基本思想是：引入一個時間局部化的窗函數(shù)，假定信號在窗函數(shù)內(nèi)是平穩(wěn)的，采用傅里葉變換分析窗內(nèi)的信號，以確定窗內(nèi)存在的頻率成分，在信號時間方向移動窗函數(shù)，得到頻率隨時間的變化關(guān)系，即所需要的時頻分布。短時傅里葉變換的定義為：

(1)

式中：x(t)為輸入信號，f是頻率，g(t)為窗函數(shù)，不同的窗函數(shù)得到的時頻分布也不相同。

1.2 小波變換(WT)

小波變換繼承和發(fā)展了短時傅里葉變換的局部化思想，同時又克服了窗口大小不隨頻率變化的缺點，是處理時變非平穩(wěn)信號比較理想的工具，該方法分別采用了窄時窗及寬時窗對高頻及低頻信號進行處理。信號f(t)的小波變換定義[8]為：

〈f(t),Ψa,b(t)〉

(2)

1.3 S變換(ST)

Stockwell[7]在前人研究時頻分析方法的基礎(chǔ)上，提出了S變換。S變換組合了短時傅里葉變換和小波變換的優(yōu)點,同時也克服了它們的不足,采用與頻率有關(guān)的可變高斯窗函數(shù),而基本小波不必滿足容許性條件。時域連續(xù)信號h(t)的S變換定義如下：

(3)

式中：h(t)為信號時間函數(shù)，f為頻率，t為時間，τ控制高斯函數(shù)在時間軸上的位置。e-i2πft是傅里葉變換為S變換的振蕩部分提供的指數(shù)函數(shù)，當τ隨著時間移動時它不隨規(guī)范化高斯窗口變化。

在S變換中，高斯窗函數(shù)和基本小波分別定義為：

(5)

由式(4)可知，S變換采用寬度可變的高斯窗函數(shù),其時窗寬度隨頻率f呈反比變化。低頻時，時窗較寬，從而獲得較高的頻率分辨率；高頻時，時窗較窄，故可獲得較高的時間分辨率。

1.4 修正S變換(MST)

對某一非平穩(wěn)信號進行時頻分解，由于S變換的窗函數(shù)是以固定的趨勢隨頻率變化的，所以獲得的時頻譜是唯一的，其時頻分辨率有限，因此需要對高斯窗進行改進[21-23]。為提高S變換對信號的時頻分辨率，在高斯窗中加入一個調(diào)節(jié)參數(shù)δ，δ為頻率的一階線性方程：

δ(f)=mf+n

(6)

其中，m、n為常數(shù)，共同控制高斯窗的寬度，m、n兩個參數(shù)均為正數(shù)，即調(diào)節(jié)參數(shù)δ隨著頻率f的增大而線性增大。

高斯窗函數(shù)則變?yōu)椋?/p>

(7)

將式(7)代入式(3)得到修正后的S變換的表達式[20]：

(8)

信號經(jīng)修正S變換得到的是一個時頻二維矩陣，行代表時間，列代表頻率。根據(jù)Heisenberg測不準原理，變換無法實現(xiàn)時頻分辨率同時提高，可以通過優(yōu)化選取參數(shù)m、n的數(shù)值，獲得較高的時頻分辨率。

2 時頻分析方法對比

為了從時間分辨率、 頻率分辨率和能量聚集等方面比較各種時頻分析方法的應(yīng)用效果，利用短時傅里葉變換、小波變換、S變換以及修正S變換方法分別對合成信號進行時頻分析。

圖1是不同頻率的雷克子波合成的理論地震記錄，主頻從20～70Hz，得到的時頻譜如圖2。從圖2a～圖2d可以看出， STFT的分辨率較為單一，而且成像模糊，能量較為分散。雖然小波變換和S變換都具有多分辨率分析的特點，但對比兩者發(fā)現(xiàn)小波變換的時間分辨率不如S變換。選取參數(shù)m=1,n=12得到信號的修正S變換時頻譜，和其他三種方法進行對比，發(fā)現(xiàn)修正S變換較其他三種變換在時間和頻率分辨率上皆有提高，高低頻主頻明顯分辨，且在低頻(10～20Hz)處，MST的時間分辨率提高顯著，能量也更為集中，消除了相鄰頻率造成的干擾，效果更好。

圖3是人工合成的地震信號，該信號由四種頻率成分組成：15Hz、20Hz、30Hz、50Hz，采樣點數(shù)為1024，采樣時間間隔是0.002s。對比4種方法的時頻譜圖(圖4a～圖4d)可以看出：短時傅里葉變換雖然可以給出信號的聯(lián)合時頻特征， 但由于時間窗函數(shù)固定，直接導致其時間和頻率分辨率是固定不變的，整體上呈現(xiàn)分辨率較低現(xiàn)象，這對低頻和高頻信號的檢測是很不利的；而小波變換結(jié)果則呈現(xiàn)多分辨率分析的特點，但頻散較為嚴重，與S變換相比存在一定的差異；由于S變換采用的高斯窗函數(shù)可以根據(jù)信號的頻率不同做不同調(diào)整，即高頻段采用高斯窗函數(shù)較窄，低頻段采用高斯窗函數(shù)較寬，因此在時頻譜上與短時傅里葉變換和小波變換相比，其時間分辨率和頻率分辨率都有明顯的改善；而修正S變換作為S變換的一種改進，通過調(diào)整參數(shù)m、n的值來控制高斯窗口的寬度，使得其時頻譜在低頻段的時間分辨率得到與S變換相比有明顯提高(如在0.4～0.5s，兩個能量團分開)，相鄰頻率的干擾有一定的抑制作用。因此，修正S變換在時間分辨率和頻率分辨率上較短時傅里葉變換、S變換和小波變換等線性時頻分析方法具有更佳的時頻聚集性，可以更好地適用于非平穩(wěn)地震信號的分析。

3 實際資料計算對比

圖5是實際采集的一道地震記錄，該地震記錄有5000個采樣點，采樣時間間隔為0.4ms。分別用STFT、WT、ST和MST時頻分析方法對該記錄進行了時頻特性分析(圖6a～圖6d)。從圖6a～圖6d中可以看出，短時傅里葉變換和小波變換雖然容易得到時頻分析譜，但是它們的時頻分析精度較低，短時傅里葉的分辨率單一。S變換采用的是隨頻率可變的窗函數(shù)，分別利用寬時窗和窄時窗在低高頻時分析信號，得到的時頻分析譜(圖6c)較好地改善了時頻局部化精度，較STFT和WT具有更高的時頻分辨率。 MST時頻分析方法基于ST的優(yōu)良性質(zhì)， 并修正了其時窗的缺陷， 使得MST在時頻譜(圖6d)的分辨率高于STFT、WT和ST。對比圖6a～圖6c與圖6d，在圖6～圖6c中，位于1.6s處與1.7s附近的能量團黏連在一起區(qū)分不開；而通過優(yōu)化參數(shù)(m=4/3,n=10)后的MST時頻譜時間分辨率明顯提高，黏在一起的能量團不僅區(qū)分開來(1.6～1.7s處)，而且整體具有更好的時頻聚集性，這更有利于地震資料的處理和解釋。

圖1 雷克子波合成信號Figure 1 Ricker wavelet signal

(a)stft時頻圖

(b)小波時頻圖

(c)S變換時頻圖

(d)MST時頻圖圖2 雷克子波合成信號時頻譜Figure 2 Time-frequency spectrum of Ricker wavelet synthetic signal

圖3 人工合成地震信號Figure 3 Synthetic signals

圖4 人工合成地震信號時頻譜

圖5 實際單道地震記錄Figure 5 Actual single channel seismogram

(a)stft時頻圖

(b)小波時頻圖

(c)S變換時頻圖

(d)MST時頻圖圖6 實際地震記錄時頻譜Figure 6 Time-frequency spectrum of actual seismic record

4 結(jié)論

本文研究了修正S變換與常規(guī)時頻方法在時頻分析中的效果對比。修正S變換利用一個頻率相關(guān)的一階線性方程δ(f)代替S變換中的歸一化高斯窗口的頻率f，對傳統(tǒng)S變換高斯窗函數(shù)相對固定的缺陷進行了修改。通過優(yōu)化選取兩個實數(shù)參數(shù)m、n的值，靈活地調(diào)控高斯窗函數(shù)寬度的變化速度。文中對雷克子波模擬信號、合成地震信號以及實際地震資料進行了多種時頻分析，結(jié)果表明:相比較常規(guī)的時頻分析方法，修正S變換能夠?qū)π盘栠M行多分辨率分析，可以根據(jù)待分析信號的特點選擇參數(shù)，具有更好的靈活性和更高的時頻分辨率，改善了時頻分布的聚集性，大大提高了地震信號時頻譜的可讀性和適應(yīng)性，能更好地應(yīng)用到不同需求的實際情況中。

[1]Boashash B,White L B. Instantaneous frequency estimation and automatic time-varying filtering [J]. Acoustics, Speech, and Signal Processing, 1990. ICASSP-90, 1990 International Conference on. IEEE, 1990: 1221-1224.

[2]劉葵, 劉招君. 時頻分析在石油地球物理勘探中的應(yīng)用[J]. 世界地質(zhì), 2000, 19(3): 282-285.

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

[4]鄧攻. S變換在地震時頻分析中的對比研究[J]. 科技創(chuàng)新與生產(chǎn)力, 2011 (4): 73-76.

[5]黃忠來, 張建中. 同步擠壓S變換[J]. 中國科學: 信息科學, 2016, 5: 007.

[6]Gabor D. Theory of Communication[J]. Journal of the IEEE.1946, 93:429-497.

[7]Stockwell R G, Mansinha L, Lowe R P. Localization of the complex spectrum: the S transform[J]. IEEE Transactions on signal processing, 1996, 44(4): 998-1001.

[8]Daubechies I. The wavelet transform, time-frequency localization and signal analysis[J]. IEEE Transactions on information theory, 1990, 36(5): 961-1005.

[9]Ville J. Théorie et applications de la notion de signal analytique[J]. Cables et transmission, 1948, 2(1): 61-74.

[10]胡明順, 潘冬明, 徐紅利, 等. 幾種時頻分析方法對比及在煤田地震勘探中的應(yīng)用[J]. 物探與化探, 2009, 33(6): 691-695+ 709.

[11]Livanos G,Ranganathan N,Jiang J. Heart sound analysis using the S-transform [J]. IEEE Computers in Cardiology,2000,27:587-590.

[12]Mansinha L, Stockwell R G, Lowe R P. Pattern analysis with two-dimensional spectral localisation: Applicationsof two-dimensional S transforms[J]. Physica A: Statistical Mechanics and its Applications, 1997, 239(1): 286-295.

[13]崔偉雄, 李德春, 徐榮華. 煤田巖溶塌陷正演模擬及屬性應(yīng)用[J]. 物探與化探, 2011, 35(5): 648-651.

[14]張繁昌, 李燦燦, 徐旺林. 改進 Stockwell 變換法識別碳酸鹽巖溶洞——以鄂爾多斯盆地奧陶系溶洞發(fā)育段為例[J]. 石油學報, 2015, 36(2): 182-187.

[15]刁瑞, 單聯(lián)瑜, 尚新民, 等. 微地震監(jiān)測數(shù)據(jù)時頻域去噪方法[J]. 物探與化探, 2015, 39(1): 112-117.

[16]Mansinha L, Stockwell R G, Lowe R P, et al. Local S-spectrum analysis of 1-D and 2-D data[J]. Physics of the Earth and Planetary Interiors, 1997, 103(3): 329-336.

[17]Pinnegar C R, Mansinha L. Time-local spectral analysis for non-stationary time series: the S-transform for noisy signals[J]. Fluctuation and Noise Letters, 2003, 3(03): L357-L364.

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

[19]陳學華, 賀振華, 黃德濟. 基于廣義 S 變換的信號提取與抑噪[J]. 成都理工大學學報 (自然科學版), 2006, 33(4): 331-335.

[20]Li D, Castagna J. Modified S-transform in time-frequency analysis of seismic data[C]//2013 SEG Annual Meeting. Society of Exploration Geophysicists, 2013.

[21]趙淑紅,朱光明. S 變換時頻濾波去噪方法[J]. 石油地球物理勘探, 2007, 42(4): 402-406.

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

[23]陳學華, 賀振華. 改進的 S 變換及在地震信號處理中的應(yīng)用[J]. 數(shù)據(jù)采集與處理, 2005, 20(4): 449-453.

Comparison of Modified S-transform (MST) and Traditional Time-Frequency Analysis Methods

Ruan Qingqing, Zhang Huixing, Wang Hao and Li Kairui

(Key Laboratory of Submarine Geosciences and Prospecting Techniques, Ministry of Education,Ocean University of China, Qingdao 266100)

Seismic signals are usually non-linear and non-stationary. The traditional Fourier transform based on stationary signal processing theory is difficult to depict frequency components at any moment. However, the time-frequency analysis technique has the capability of describing the partial features of signal both in time and frequency domains. The modified S-transform is a new method of time-frequency analysis, which was proposed due to the problems of the window function invariant and time-frequency resolution unadjustable in the S-transform. In this paper, we study and analyze the modified S-transform by comparing with the traditional time-frequency decomposition methods, such as short time Fourier transform (STFT), wavelet transform (WT) and S-transform. The time-frequency results of synthetic signals and field seismic data show that: the modified S transform has higher time and frequency resolution and better energy aggregation than the traditional time-frequency analysis methods. Therefore, it is more conducive to the processing of non-stationary signals.

time-frequency analysis; modified S-transform; S-transform; wavelet transform; short time Fourier transform

10.3969/j.issn.1674-1803.2017.04.14

1674-1803(2017)04-0066-07

國家科技重大專項(2016ZX05027-002)；國家自然科學基金(41204089)；大學生創(chuàng)新訓練項目(201510423056X)。

阮清青(1991—)，男，中國海洋大學2014級碩士研究生，主要從事時頻分析方法研究。

張會星，博士，副教授，Email:zhhuixing@sina.com

2017-03-10

文獻標識碼：A

責任編輯：孫常長