S變換對(duì)地震信號(hào)濾波的應(yīng)用研究

徐濤

（武夷學(xué)院機(jī)電工程學(xué)院，福建武夷山354300）

S變換對(duì)地震信號(hào)濾波的應(yīng)用研究

徐濤

（武夷學(xué)院機(jī)電工程學(xué)院，福建武夷山354300）

短時(shí)傅里葉變換的窗函數(shù)寬度固定不變，處理非平穩(wěn)信號(hào)有局限性，改變窗函數(shù)使其可調(diào)，即為S變換。S變換的窗函數(shù)可以隨信號(hào)自適應(yīng)調(diào)節(jié)，信號(hào)經(jīng)S變換得二維時(shí)頻譜，確定噪聲時(shí)頻范圍對(duì)其清零，再利用S逆變換到時(shí)間域，最終得到去噪后的有效信號(hào)。通過理論計(jì)算和地震信號(hào)仿真表明，S變換能夠得到信號(hào)的時(shí)頻信息，確定噪聲的時(shí)間范圍，并能有效濾除不同時(shí)段不同頻率的噪聲，提高信號(hào)的信噪比。

S變換；窗函數(shù)；地震信號(hào)

經(jīng)典傅里葉變換是處理穩(wěn)態(tài)信號(hào)的有效方法，可以立即得出信號(hào)的頻譜信息，但無法反映出信號(hào)頻率隨時(shí)間的變化情況，基于此，提出了窗函數(shù)的概念，使時(shí)域信號(hào)在此窗內(nèi)能夠體現(xiàn)頻率信息，這種處理方法稱為時(shí)頻分析方法[1]，其中，CWT、STFT是比較有代表性的時(shí)頻分析方法。但時(shí)頻分析法的局限性在于無法有效處理變化劇烈的非平穩(wěn)信號(hào)，且時(shí)頻分析方法的窗函數(shù)仍受到Heisenberg W不確定性原理的限制，即時(shí)頻分辨率無法共同達(dá)到最優(yōu)。

地球物理學(xué)家Stockwell R G在1996年提出了一種加時(shí)窗的時(shí)頻可逆分析方法，也就是S變換[2]。S變換是在STFT基礎(chǔ)上對(duì)窗函數(shù)進(jìn)行改進(jìn)[3-4]，采用CWT思想，以Morlet小波（具有非正交性和Gaussian調(diào)節(jié)的指數(shù)賦值小波）為基本小波，根據(jù)不同時(shí)間區(qū)域信號(hào)的劇烈程度，確定合適的時(shí)窗分辨率，從而計(jì)算出不同時(shí)刻的功率譜，確定噪聲存在的區(qū)域，對(duì)該區(qū)域進(jìn)行清零，再經(jīng)S逆變換到時(shí)間域，最終得到濾波后的有效信息。與STFT、CWT相比，S變換有其獨(dú)特的優(yōu)勢(shì)[5-6]：窗函數(shù)分辨率具有多分辨性，寬度又可以隨信號(hào)頻率自適應(yīng)地做出改變，其Morlet小波可以不受容許性條件限制等[7]。

1 S變換

1.1 由傅里葉變換導(dǎo)出

設(shè)為連續(xù)時(shí)間序列，其傅里葉變換表達(dá)式可以表示為：H（f）=由于傅里葉變換無法分析非平穩(wěn)信號(hào)的局部信息，也無法體現(xiàn)信號(hào)頻率隨時(shí)間的變化情況，針對(duì)這些缺陷，在1946年，Gabor率先提出了加窗函數(shù)的思想，即在h（t）上乘以一個(gè)窗函數(shù)，設(shè)為w（t），用w（t）按時(shí)間軸去截每段信號(hào)，再對(duì)截下來的每段信號(hào)作傅里葉變換，所得到的一系列傅里葉變換結(jié)果排開則成為二維表象，這些傅里葉變換的集合記為：H（f）=（2），可以證明，窗函數(shù)的時(shí)頻寬度的乘積存在最小值，其最小值為2，這恰恰符合了測(cè)不確定性原理，該原理說明了不可能同時(shí)獲得最佳的時(shí)頻分辨率?；诖耍瑢⒋昂瘮?shù)歸一化，定義為高斯窗，且該高斯窗可以平移和伸縮，平移量為τ、伸縮量為σ，那么，w（t-τ）=此時(shí)時(shí)間序列h（t）的時(shí)間頻率二維表象又可寫為：

（4），上式除了時(shí)間和頻率二維表象元素外，還多了一個(gè)伸縮變量σ，作為時(shí)頻分析工具是不切實(shí)際的?？紤]將伸縮變量與頻率產(chǎn)生關(guān)系，令，這樣就得到一個(gè)公式：S（t，f）=exp（-j2πft）dt（5），這就是連續(xù)時(shí)間序列h（t）的S變換表達(dá)式。根據(jù)傅里葉逆變換很容易得到：該式為S逆變換表達(dá)式。

1.2 離散S變換及其實(shí)現(xiàn)

由于計(jì)算機(jī)處理的是離散信號(hào)，而在離散域內(nèi)，又必須滿足周期性或有限性條件，所以必須要導(dǎo)出有限序列的離散S變換。

h（KT），K=0,1,2,…,N為連續(xù)信號(hào)h（t）的離散信息，T為采樣時(shí)間間隔，其離散傅里葉變換可以表示為：N，時(shí)間序列h（KT）的離散S變換可以表示為：（8），式中N為采樣點(diǎn)數(shù)，j為正整數(shù)且無單位，m和n為頻率，是h的離散傅里葉變換。離散S反變換為：中，n≠0。

1.3 S變換的實(shí)現(xiàn)

首先記n=n/NT,m=m/NT,k=kT和j=jT。序列h（k）的S正變換步驟如下：

（1）在時(shí)間序列h（k）上進(jìn)行采樣，設(shè)采樣點(diǎn)數(shù)為N，采樣時(shí)間間隔為T，對(duì)這N個(gè)點(diǎn)求FFT（快速傅里葉變換），得到其頻譜，記為H[m]。

（2）計(jì)算頻率為n的高斯函數(shù)值G[m,n]。

（3）將H[m]向左移動(dòng)n個(gè)頻率單位得H[m+n]，即頻移譜。

（4）將G[m,n]乘以H[m+n]所得結(jié)果記作B[n,m]。

（5）將B[n,m]做快速傅里葉逆變換，得其S變換譜表達(dá)式S[n,j]，j=1,2,…,N。

（6）返回執(zhí)行步驟（3）、（4），當(dāng)所有頻率全部計(jì)算完畢，就得到了h(K)的S變換譜。

同樣，可以根據(jù)S逆變換的離散形式得出S逆變換步驟：

（1）對(duì)譜S[n，j]按行求和，計(jì)算出頻率n的傅里葉譜，記為C[n]。

（2）循環(huán)頻率n，返回再執(zhí)行步驟（1），計(jì)算出譜H[m]。

（3）對(duì)H[m]執(zhí)行傅里葉逆變換運(yùn)算，還原得出時(shí)間序列h（k）。

2 S變換在時(shí)頻域上的濾波

2.1 濾波方法簡(jiǎn)介

短時(shí)傅里葉變換在對(duì)信號(hào)進(jìn)行濾波時(shí)，是比較常用的方法，對(duì)比傅里葉變換最大的優(yōu)點(diǎn)就是，添加了可伸縮可平移的高斯窗函數(shù)，能夠獲取信號(hào)在各個(gè)時(shí)刻對(duì)應(yīng)的頻譜信息，但是高斯窗函數(shù)不夠靈活，缺乏自適應(yīng)性，因此無法跟蹤信號(hào)的突變過程，而信號(hào)的噪聲部分又往往發(fā)生在該過程，即信號(hào)突變區(qū)，所以設(shè)置隨信號(hào)頻率變化的自適應(yīng)高斯窗函數(shù)顯得尤為重要，用自適應(yīng)的高斯窗函數(shù)對(duì)信號(hào)進(jìn)行截取，得到不同時(shí)間對(duì)應(yīng)的頻域信息，在時(shí)間-頻率這個(gè)二維平面上就確定了信號(hào)噪聲的存在的時(shí)間頻率范圍，簡(jiǎn)單來說就是將信號(hào)分為有效部分和噪聲部分，利用S變換，得到信號(hào)的時(shí)頻譜，確定噪聲區(qū)域，將噪聲區(qū)域清零，然后經(jīng)逆S變換回時(shí)間域，具體做法如下，將信號(hào)表示為originalsignal(t)=realsignal(t)+noise(t)，等式右邊為原始信號(hào)，左邊為有效信號(hào)和噪聲信號(hào)，將兩邊進(jìn)行S變換得：SToriginalsignal(t)=STrealsignal(t)+STnoise(t)。在信號(hào)S譜上確定噪聲區(qū)域，對(duì)該區(qū)域清零，也就是讓STnoise(t)為0，最后將STrealsignal(t)經(jīng)S逆變換回時(shí)間域，最終得到清除干擾后的近似有效信號(hào)realsignal(t)。

2.2 濾波步驟

（1）首先對(duì)待處理信號(hào)進(jìn)行S變換，得其S譜，即時(shí)間-頻率信息，如圖1。

圖1 信號(hào)S變換時(shí)頻圖Figure 1 Time frequency diagram of signal’s S-transform

（2）確定噪聲區(qū)域，即時(shí)間范圍[t1，t2]和頻率范圍[f1，f2]，將此二維時(shí)間-頻率區(qū)域標(biāo)為D，注意，有效信息和噪聲信息的邊界是一個(gè)敏感區(qū)域（誤差產(chǎn)生區(qū)域），邊界可以通過對(duì)信號(hào)分時(shí)分頻計(jì)算來確定。

（3）對(duì)噪聲區(qū)域進(jìn)行清零。

（4）對(duì)清零后數(shù)據(jù)進(jìn)行S逆變換，最終得到濾除噪聲后的時(shí)間域信息。

3 S變換與小波變換模型試算

圖2是一幅原始的待處理的信號(hào)，由人工疊加噪聲合成，為時(shí)間域信息；圖3為該信號(hào)經(jīng)S變換后的時(shí)頻譜，得到由低到高的三段頻率信息，分別為0.05、0.1和0.15 Hz，其中，0.05 Hz的低頻信號(hào)上疊加了高斯白噪聲，噪聲功率約為0.6，此時(shí)的信噪比為7.923 5；原始信號(hào)的S變換譜中可以清楚的看到三段不同頻率的信號(hào)以及它們存在的時(shí)間區(qū)域，同時(shí)也能夠清楚的看到噪聲所在的時(shí)頻區(qū)域（頻率區(qū)域大概為0.2～ 0.6 Hz，時(shí)間區(qū)域大概為65～85 s），對(duì)這部分區(qū)域進(jìn)行清零，將噪聲清除，然后將S變換后的時(shí)頻譜反變換回時(shí)間域，得到去噪后的近似有效信號(hào)，如圖4，可以看出通過S變換后濾波并不是完全徹底，這是由于在對(duì)噪聲時(shí)頻范圍的確定時(shí)存在誤差信息，濾波后信號(hào)除開始部分、與噪聲產(chǎn)生的頻率段產(chǎn)生些許畸變外，其余部分與原信號(hào)基本一致，經(jīng)計(jì)算，濾波后信號(hào)信噪比約為12.355 6，大大提高了信噪比。圖5為信號(hào)經(jīng)過小波變換濾波后的結(jié)果，經(jīng)過計(jì)算，經(jīng)小波變換濾波后的信號(hào)信噪比為11.423 3，對(duì)比得出，S變換濾波效果要優(yōu)于小波變換濾波結(jié)果。

圖2 人工合成噪聲的原始信號(hào)Figure 2 Artificial synthetic raw noise signal

圖3 原始信號(hào)經(jīng)S變換后的時(shí)頻譜Figure 3 Time spectrum after S-transform

圖4 經(jīng)S變換濾波后的時(shí)域信息Figure 4 Time domain signal after S-transform filtering

圖5 經(jīng)小波變換濾波后的時(shí)域信息Figure 5 Time domain signal after wavelet transform filtering

4 S變換對(duì)實(shí)際地震信號(hào)去噪

圖5為某一地區(qū)的實(shí)際地震記錄資料，選取60道地震信號(hào)進(jìn)行濾波分析，其中每一道的道長(zhǎng)為3 s，信息的采樣點(diǎn)數(shù)為1 500，采樣間隔為2 ms，從圖中可以看出每一道信號(hào)在低頻的基礎(chǔ)上疊加了高頻信息，分別取出每一道的地震信號(hào)，進(jìn)行S變換，確定噪聲存在時(shí)頻區(qū)域，對(duì)噪聲進(jìn)行清零，然后返回時(shí)間域，直到處理完所有的60道信號(hào)，最終得到去噪后的60道地震時(shí)域信息如圖6。

圖6 原始地震信號(hào)Figure 6 Original seismic signal

圖7 去噪后的地震信號(hào)Figure 7 Seismic signal after denoising

與圖5比較可以明顯看出這60道信號(hào)的高頻部分干擾已經(jīng)基本完全濾除，信號(hào)相對(duì)平穩(wěn)與清晰。將60道信號(hào)的干擾信息完全濾除后，從宏觀角度來看，這60道信息的干擾分量已經(jīng)基本完全濾除，地震信息呈現(xiàn)低頻狀態(tài)且比較平穩(wěn)，由此來看，基于S變換的去噪方法能夠很好地對(duì)實(shí)際地震信息進(jìn)行濾波。

5 結(jié)論

對(duì)比短時(shí)傅里葉變換S變換的優(yōu)勢(shì)在于，其窗函數(shù)可以根據(jù)信號(hào)特征調(diào)整其高度和寬度；對(duì)比小波變換其優(yōu)勢(shì)在于，信號(hào)的S變換結(jié)果更加直觀，能夠更加細(xì)致地對(duì)高頻區(qū)進(jìn)行分析。通過模型試算對(duì)比，S變換濾波方法要優(yōu)于小波變換濾波方法。雖然S變換的窗函數(shù)可調(diào)，但是其形態(tài)還是保持不變，僅限于平移與伸縮，這樣也就制約了S變換對(duì)信號(hào)分析的靈活度，另外，窗函數(shù)時(shí)頻分辨率依然受測(cè)不確定性原理的約束。所以，基于S變換去噪方法仍有待于完善。

[1]焦敘明.時(shí)頻分析及其在地震資料處理分析中的應(yīng)用[D].青島：中國(guó)海洋大學(xué),2007.

[2]Stockwell R G,Mansinha L,Lowe R P.Location of the complex spectrum：the stransform[J].IEEE Transactions on Signal 1996;44(4):998-1001.

[3]楊志強(qiáng),單娜林,劉占興,等.S變換在巖溶區(qū)地震映像資料處理中的應(yīng)用[J].工程地球物理學(xué)報(bào),2012（2）：227-230.

[4]遲華山,王紅星,趙培紅,等.基于S變換的線性調(diào)頻信號(hào)時(shí)頻濾波[J].無線電通信技術(shù),2012（1）：21-24.

[5]趙淑紅,王璇.S變換時(shí)頻濾波與其它濾波方法的比較[J].地震工程學(xué)報(bào),2007（3）：224-229.

[6]劉霞,徐濤,段玉波,等.基于廣義S變換的時(shí)頻濾波技術(shù)研究[J].自動(dòng)化技術(shù)與應(yīng)用,2012（2）：15-19.

[7]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(3):357-364.

（責(zé)任編輯：葉麗娜）

Application Study on S-transform to Seismic Signal Filteration

XU Tao
(School of Mechanical and Electrical Engineering,Wuyi University,Wuyishan,Fujian 354300)

The width of the short time Fourier transform window function is invariable.There are limitations in dealing with non-stationary signals.Change window function and make it adjustable,this is S-transform.The window function of S-transform can be adjusted with signals.The signal is transformed to be two dimensional time-frequency spectrum by S-transform.Determine the time-frequency range of noise and clear the noise.it is transformed to be time domain by S-inverse-transform.Finally,obtain the effective and filtered information. Through the theoretical calculation and simulation of seismic signal,S-transform can get the time frequency information of signal,determine time scope of the noise,improve the SNR of signal.

S transform;window function;seismic signal

TN911.72

A

1674-2109（2016）12-0059-04

2015-09-12

徐濤(1984-)男,漢族,助教，主要從事故障診斷與容錯(cuò)控制的研究。