S變換對地震信號濾波的應用研究

徐濤

（武夷學院機電工程學院，福建武夷山354300）

S變換對地震信號濾波的應用研究

徐濤

（武夷學院機電工程學院，福建武夷山354300）

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

S變換；窗函數；地震信號

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

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

1 S變換

1.1 由傅里葉變換導出

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

（4），上式除了時間和頻率二維表象元素外，還多了一個伸縮變量σ，作為時頻分析工具是不切實際的。考慮將伸縮變量與頻率產生關系，令，這樣就得到一個公式：S（t，f）=exp（-j2πft）dt（5），這就是連續時間序列h（t）的S變換表達式。根據傅里葉逆變換很容易得到：該式為S逆變換表達式。

1.2 離散S變換及其實現

由于計算機處理的是離散信號，而在離散域內，又必須滿足周期性或有限性條件，所以必須要導出有限序列的離散S變換。

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

1.3 S變換的實現

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

（1）在時間序列h（k）上進行采樣，設采樣點數為N，采樣時間間隔為T，對這N個點求FFT（快速傅里葉變換），得到其頻譜，記為H[m]。

（2）計算頻率為n的高斯函數值G[m,n]。

（3）將H[m]向左移動n個頻率單位得H[m+n]，即頻移譜。

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

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

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

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

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

（2）循環頻率n，返回再執行步驟（1），計算出譜H[m]。

（3）對H[m]執行傅里葉逆變換運算，還原得出時間序列h（k）。

2 S變換在時頻域上的濾波

2.1 濾波方法簡介

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

2.2 濾波步驟

（1）首先對待處理信號進行S變換，得其S譜，即時間-頻率信息，如圖1。

圖1 信號S變換時頻圖Figure 1 Time frequency diagram of signal’s S-transform

（2）確定噪聲區域，即時間范圍[t1，t2]和頻率范圍[f1，f2]，將此二維時間-頻率區域標為D，注意，有效信息和噪聲信息的邊界是一個敏感區域（誤差產生區域），邊界可以通過對信號分時分頻計算來確定。

（3）對噪聲區域進行清零。

（4）對清零后數據進行S逆變換，最終得到濾除噪聲后的時間域信息。

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

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

圖2 人工合成噪聲的原始信號Figure 2 Artificial synthetic raw noise signal

圖3 原始信號經S變換后的時頻譜Figure 3 Time spectrum after S-transform

圖4 經S變換濾波后的時域信息Figure 4 Time domain signal after S-transform filtering

圖5 經小波變換濾波后的時域信息Figure 5 Time domain signal after wavelet transform filtering

4 S變換對實際地震信號去噪

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

圖6 原始地震信號Figure 6 Original seismic signal

圖7 去噪后的地震信號Figure 7 Seismic signal after denoising

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

5 結論

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

[1]焦敘明.時頻分析及其在地震資料處理分析中的應用[D].青島：中國海洋大學,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]楊志強,單娜林,劉占興,等.S變換在巖溶區地震映像資料處理中的應用[J].工程地球物理學報,2012（2）：227-230.

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

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

[6]劉霞,徐濤,段玉波,等.基于廣義S變換的時頻濾波技術研究[J].自動化技術與應用,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.

（責任編輯：葉麗娜）

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-)男,漢族,助教，主要從事故障診斷與容錯控制的研究。