郝國成 談帆 程卓 王巍 馮思權(quán) 張偉民

時(shí)頻分析(Time-frequency analysis,TFA)是從時(shí)間–頻率的角度來處理非平穩(wěn)信號的方法,其基本思想是建立時(shí)間和頻率的聯(lián)合函數(shù),描述信號在時(shí)間–頻率平面的能量密度或強(qiáng)度[1?2].它將信號從時(shí)間域變換到時(shí)間–頻率聯(lián)合域,不同頻率分量的時(shí)間關(guān)聯(lián)特性能夠在時(shí)–頻平面上有效地表示出來.信號在任意時(shí)刻的能量都聚集在此瞬時(shí)頻率附近,基于線性方法的逆變換則可以重構(gòu)其等效的時(shí)間域信號[3].目前時(shí)頻分析方法已經(jīng)廣泛應(yīng)用于自動化控制、信號處理、數(shù)據(jù)通信、地質(zhì)勘探、基礎(chǔ)物理、工業(yè)生產(chǎn)等各個領(lǐng)域[4?5].在合成人造金剛石加工環(huán)節(jié),自動檢測頂壓機(jī)頂錘是否破裂是該項(xiàng)生產(chǎn)過程的重要步驟之一,可有效降低人工憑經(jīng)驗(yàn)聽音辨別的誤差.該金屬破裂信號屬于典型的非平穩(wěn)信號,利用合適的時(shí)頻分析方法能夠?qū)Υ祟愋盘栠M(jìn)行有效的時(shí)頻表示,針對信號的破裂局部信息給出對應(yīng)的高幅值頻率分布,為數(shù)據(jù)采集卡的頻率判別窗口提供閾值依據(jù).本文提出基于二值化Gabor的歸一化Wigner-Ville分布(Wigner-Ville distrivution,WVD)和歸一化偽平滑WVD(BGabor-NWVD,BGabor-NSPWVD)時(shí)頻分析算法,具有良好的高銳化時(shí)頻聚集度和魯棒性,對金屬破裂樣本信號有較好的去噪效果,可以得出有參考意義的時(shí)間和頻率聯(lián)合分布.

1 時(shí)頻分析

時(shí)頻分析通過設(shè)計(jì)時(shí)間和頻率的幅度或能量密度關(guān)聯(lián)函數(shù),將一維的時(shí)間序列信號以二維的時(shí)間–頻率密度函數(shù)形式表示,旨在揭示信號中包含多少頻率分量,以及每一分量隨時(shí)間的變化規(guī)律.時(shí)頻分析方法可以分為線性型時(shí)頻表示、非線性型時(shí)頻表示和二次型時(shí)頻表示.1)典型的線性時(shí)頻表示有短時(shí)傅里葉變換(Short-time Fourier transform,STFT)[6?7]、S 變換(S transform,ST)[8]、和小波變換(Wavelet transform,WT)[9]等,線性時(shí)頻表示會受到不確定性原理的制約,存在時(shí)頻模糊等問題.2)非線性型時(shí)頻表示包括基于經(jīng)驗(yàn)?zāi)Ｐ头纸?Empirical mode decomposition,EMD)的希爾伯特–黃變換 (Hilbert-Huang transformation,HHT)[10]、集合經(jīng)驗(yàn)?zāi)Ｐ头纸?Ensemble empirical mode decomposition,EEMD)及其改進(jìn)算法[11?12]、同步壓縮變換(Synchrosqueezing transform,SST)[13]等.基于EMD和EEMD的HHT具有良好的窄帶自適應(yīng)性,但缺乏嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)支撐,且信號頻率帶寬較大時(shí),其時(shí)頻聚集度有待提高.SST方法是Daubechies于2011年提出的一種時(shí)頻重排技術(shù),其變換核可以基于小波變換或者基于STFT變換,SST支持信號重構(gòu),其時(shí)頻輸出具有高銳化聚集度的特點(diǎn)[13],但直接SST存在復(fù)雜頻率成分交叉點(diǎn)附近的描述模糊和魯棒性弱的缺點(diǎn).改善SST的魯棒性,可采用與壓縮感知稀疏方法[14]相結(jié)合的方法,此部分本文作者另有相關(guān)文章闡述.3)二次型時(shí)頻分布則是一種更為嚴(yán)格的時(shí)頻表示,能夠描述信號的瞬時(shí)功率譜密度,可從時(shí)間–頻率–能量三者聯(lián)合分布的維度來了解信號的特點(diǎn).二次型時(shí)頻分布一般包括Cohen類時(shí)頻分布[6]、WVD[15],以及對WVD方法的改進(jìn),如STFT-WVD[16]、STFT-SPWVD[17]和NSTFT-WVD[18]等.

常用的STFT方法由Gabor于1946年提出[7],即加窗的傅里葉變換.窗函數(shù)可以選擇矩形窗、三角窗、漢寧(Hanning)窗、海明(Hamming)窗、布拉克曼(Blackman)窗、愷撒(Kaiser)窗等,當(dāng)STFT的窗函數(shù)選擇高斯窗時(shí),加窗短時(shí)傅里葉變換稱為Gabor變換.根據(jù)Heisenberg測不準(zhǔn)原理,高斯函數(shù)窗口面積已達(dá)到測不準(zhǔn)原理下界,是時(shí)域窗口面積達(dá)到最小的函數(shù),Gabor變換是最優(yōu)的STFT.文獻(xiàn)[19]采用稀疏分析窗的方法來求解離散Gabor變換,一定程度提高了時(shí)頻的聚集度.STFT和Gabor變換使用的是大小和形狀固定的滑動窗口,變換基函數(shù)為非正交系,對于突變和非平穩(wěn)信號,不能精確分解周期比時(shí)間窗大的低頻信號,且高頻的時(shí)頻分辨率比較差,需要輔以其他高時(shí)頻聚集度方法加以改進(jìn),WVD就是較好的選擇,該方法具有理論上最高時(shí)頻分辨率和許多優(yōu)良的數(shù)學(xué)性質(zhì).

WVD是一種基礎(chǔ)并且重要的二次型時(shí)頻分布,具有較高的時(shí)頻聚集度,應(yīng)用廣泛.Wigner分布在1932年被提出,并應(yīng)用于量子力學(xué),但沒有引起重視,直到Ville在1948年將它應(yīng)用于信號分析領(lǐng)域,并稱之為Wigner-Ville分布,開始了它的快速發(fā)展時(shí)期.尤其自上世紀(jì)八十年代以來,陸續(xù)發(fā)表了相當(dāng)數(shù)量的WVD論文,對WVD的定義、性質(zhì)、時(shí)頻表示、各種改進(jìn)進(jìn)行了大量的討論,以其良好的時(shí)變特性,特別是可直接由其精確定義“瞬時(shí)頻率”而廣泛用于非平穩(wěn)信號的分析.同時(shí)需要注意的是,這種方法雖然對邊緣特性、瞬時(shí)頻率和局域刻畫等都有很好的描述,其時(shí)頻圖上的頻率成分有高銳化的聚焦,但WVD不滿足疊加原理,其變換過程中產(chǎn)生新的頻率交叉項(xiàng)分量成為該方法的瓶頸[18].交叉項(xiàng)的存在使得時(shí)頻圖上的分布混亂,額外產(chǎn)生大量的噪聲成分,充斥在真實(shí)頻率之間,干擾對信號的識別分析,極大地限制了WVD分布的實(shí)際應(yīng)用.

綜上所述,采用Gabor與WVD結(jié)合對非平穩(wěn)信號進(jìn)行時(shí)頻分析,可發(fā)揮各自的極限優(yōu)點(diǎn),去除WVD交叉項(xiàng)的同時(shí),得到較高的局部時(shí)頻分辨率,滿足對非平穩(wěn)信號進(jìn)行時(shí)頻分析的高銳化時(shí)頻聚集度和強(qiáng)魯棒性要求.本文對Gabor和WVD分別進(jìn)行改進(jìn)后再結(jié)合,提出BGabor-NWVD和BGabor-NSPWVD算法,同時(shí)兼顧線性時(shí)頻表示與二次型時(shí)頻表示的優(yōu)點(diǎn),能夠獲得較好的時(shí)頻聚集度、消除WVD交叉項(xiàng)干擾、抑制噪聲,并具有強(qiáng)魯棒性,同時(shí)能在金屬破裂樣本信號的時(shí)頻分析中有較好的應(yīng)用.

2 BGabor-NWVD和BGabor-NSPWVD

2.1 Gabor與WVD的優(yōu)缺點(diǎn)

STFT與Gabor變換是常用的時(shí)頻分析方法之一.其優(yōu)點(diǎn)是,對信號具有一定程度的自適應(yīng)性,能夠如實(shí)還原輸入信號的時(shí)頻分布,對信號的低頻和高頻兩端都有較為一致的局部描述,得到的非平穩(wěn)信號的時(shí)頻分布較為平滑,分析多分量信號的時(shí)候,沒有交叉項(xiàng)的虛假頻率成分干擾.STFT變換的定義如式(1),把STFT的窗函數(shù)變?yōu)楦咚购瘮?shù),就得到式(2)的Gabor變換[7].

Gabor的時(shí)頻分辨率由窗函數(shù)g(τ?t)的時(shí)域大小直接決定,一旦窗口函數(shù)選定,其時(shí)頻分辨率就已確定,不隨時(shí)間和頻率的變化而變化.為了提高時(shí)間分辨率,窗函數(shù)的時(shí)間寬度應(yīng)該盡量短,但為了提高頻率分辨率,窗函數(shù)的時(shí)間寬度則應(yīng)盡可能長,受不確定性原理的約束,時(shí)間分辨率和頻率分辨率是一對矛盾體.在處理非平穩(wěn)信號的過程中,對于高頻信息,需要用較窄的窗函數(shù)進(jìn)行分析,而對于低頻信息,則需要用較寬的窗函數(shù)進(jìn)行分析.而Gabor的窗函數(shù)確定以后,只能以一種固定分辨率進(jìn)行時(shí)頻分析,無法兼顧高頻信息和低頻信息[20?21].這種方法的缺點(diǎn)很明顯,其時(shí)頻聚集性不佳,局部的頻帶粗糙,難以精確顯示時(shí)頻特點(diǎn),如圖1所示,需要通過其他方法提高其時(shí)頻聚集度.

WVD是最基本且嚴(yán)格的時(shí)間和頻率二維聯(lián)合函數(shù),可以看作信號在時(shí)間和頻率平面上的能量密度解.雖然時(shí)頻表示的線性特性是我們所希望具備的重要特性,但因?yàn)槟芰勘旧砭褪且环N二次型表示,因此,當(dāng)反映信號的時(shí)間–能量分布時(shí),采用WVD這樣的二次型時(shí)頻分布則更加合適.WVD屬于Cohen類雙線性時(shí)頻分布的一種,可在時(shí)域和頻域同時(shí)揭示信號的能量分布,并且其物理意義明確.近年來,憑借著其優(yōu)越的時(shí)頻聚集性,WVD被廣泛應(yīng)用于信號分析和處理領(lǐng)域,尤其在非平穩(wěn)信號的瞬時(shí)頻率估計(jì)、信號的相干檢測和時(shí)變?yōu)V波等諸多領(lǐng)域,是最常用的雙線性時(shí)頻分布之一.WVD分布定義為:

其中,z(t)為目標(biāo)信號的解析信號,τ是積分變量,t是時(shí)移,f是頻率.式(3)沒有使用窗函數(shù),避免了線性時(shí)頻表示時(shí)間和頻率分辨率相互制約的矛盾.WVD的時(shí)間帶寬積達(dá)到了不確定原理給出的下界,可以呈現(xiàn)較其他時(shí)頻分布更好的時(shí)頻聚集特性.但WVD在處理復(fù)雜信號,尤其是頻率相近的加性信號的時(shí)候,由自身變換過程所引入的交叉項(xiàng)問題同樣很嚴(yán)重.設(shè)z(t)=z1(t)+z2(t),則z(t)的WVD如式(4).

圖1 信號f1的理想時(shí)頻、Gabor和WVD對比圖Fig.1 Ideal time-frequency spectrum,Gabor,WVD off1

其中,等號右邊的Wz1(t,f)與Wz2(t,f)分別為z1(t)和z2(t)的WVD,即變換所要求解的時(shí)頻信息項(xiàng),第3項(xiàng)2Re{Wz1,z2(t,f)}為交叉項(xiàng).由式(4)可知,兩個信號和的WVD并不等于它們各自WVD之和,交叉項(xiàng)的存在給WVD的應(yīng)用帶來不小的困擾.

式(5)是信號f1(t)的時(shí)域表達(dá)式,圖1為信號f1(t)理想時(shí)頻表示圖、Gabor和WVD的時(shí)頻分布圖.從圖中可以看出,Gabor變換的時(shí)頻分辨率不高,而WVD方法的分辨率盡管略高于Gabor變換,但存在嚴(yán)重的交叉項(xiàng)干擾問題.因此,單獨(dú)的Gabor和WVD都不理想,需要進(jìn)行一定的改進(jìn),以克服時(shí)頻分辨率和交叉項(xiàng)的影響.

2.2 Gabor與WVD的結(jié)合思想

Gabor-WVD變換的基本思想是利用Gabor和WVD各自的優(yōu)點(diǎn),通過兩者重疊運(yùn)算來增強(qiáng)Gabor和WVD信息項(xiàng),抑制由WVD變換產(chǎn)生的交叉項(xiàng)部分,以達(dá)到Gabor-WVD變換在保持良好的時(shí)頻聚集特性的同時(shí),具有消除交叉項(xiàng)的效果.Gabor-WVD變換定義了G(t,f)與W(t,f)這兩個過程變量,其任意函數(shù)表達(dá)式如式(6)所示[18].

其中,p(x,y)為任意函數(shù),例如當(dāng)p(x,y)=xayb時(shí),GW(t,f)=Ga(t,f)Wb(t,f);當(dāng)p(x,y)=x+y時(shí),GW(t,f)=G(t,f)+W(t,f).Gabor-WVD得出的結(jié)果可以很好地描述非平穩(wěn)信號的二維時(shí)頻分布,并且具有良好的聚集性,對于WVD產(chǎn)生的交叉項(xiàng)也有較好的抑制.但本文作者在文獻(xiàn)[18]中已經(jīng)明確,STFT-WVD不能真實(shí)反映信號的三維幅值,如本文第三部分的數(shù)值實(shí)驗(yàn)仿真例子所示.Gabor-WVD的三維輸出時(shí)頻表示存在較大的誤差和干擾,具體體現(xiàn)在信息項(xiàng)上的交叉項(xiàng)無法消除.針對交叉項(xiàng)和三維幅值不能正確顯示的問題,需要對Gabor-WVD進(jìn)一步改進(jìn).由于交叉項(xiàng)成因與變換過程中的兩個乘項(xiàng)表現(xiàn)為時(shí)頻軸上的頻率相鄰與中間頻段位置的量[18],因此,消除交叉項(xiàng)的思路,可采用設(shè)定閾值來清除指定頻段上的虛假分量,或是利用其他合適的數(shù)值進(jìn)行替代,或是對該項(xiàng)增添冪指數(shù)進(jìn)行消除.實(shí)現(xiàn)Gabor-WVD的若干種方法如式(7)～(9).

根據(jù)函數(shù)p(x,y)形式的不同,GW(t,f)的運(yùn)算方式也有所不同.由于式(7)只取Gabor變換與WVD后的數(shù)值中的較小值,故稱為最小值法.式(8)將c設(shè)置為交叉項(xiàng)消除閾值,對Gabor數(shù)組中的部分特定數(shù)據(jù)取0或1,故稱為二值化法.式(9)設(shè)置a、b為冪指數(shù),通過冪指數(shù)來消除交叉項(xiàng),故稱為冪調(diào)節(jié)系數(shù)法.

2.3 BGabor-NWVD和BGabor-NSPWVD

盡管Gabor-WVD方法能夠有效消除WVD產(chǎn)生的虛假分量,但并不能消除疊加在信息項(xiàng)上的分量.針對此問題,本文對Gabor和WVD同時(shí)進(jìn)行改進(jìn),提出BGabor-NWVD和BGabor-NSPWVD算法.BGabor-NWVD算法的步驟為:

步驟1.首先各自對輸入信號進(jìn)行Gabor和WVD變換,分別得到數(shù)組A和數(shù)組B.

步驟2.根據(jù)時(shí)域波形的幅度,對Gabor數(shù)組進(jìn)行二值化處理(Binarization),得到數(shù)組Gabor1,同時(shí),根據(jù)WVD絕對值數(shù)組中的最大值對WVD數(shù)組進(jìn)行歸一化處理(Normalization),得到數(shù)組WV D_1.

步驟3.兩數(shù)組進(jìn)行點(diǎn)除,得到新的臨時(shí)數(shù)組Y,根據(jù)其本身的大小,設(shè)置閾值K1、K2.

步驟4.修正點(diǎn)除后的結(jié)果數(shù)組B,記錄數(shù)組Y中大于K1的位置,在數(shù)組B中將相同位置的元素置0,得到數(shù)組B0.

步驟5.對數(shù)組Y進(jìn)行修正,將數(shù)組Y中大于K1的元素置1,將數(shù)組Y中小于K2的元素也置1,得到新數(shù)組Y0.

步驟6.數(shù)組B0點(diǎn)除Y0,輸出BGabor-NWVD數(shù)組,算法結(jié)束.

數(shù)值實(shí)驗(yàn)結(jié)果表明,這一改進(jìn)算法對信息項(xiàng)上交叉項(xiàng)的疊加有較好的抑制效果.SPWVD是對WVD的一種加窗平滑改進(jìn),SPWVD本身即具有消除WVD交叉項(xiàng)的作用,將BGabor與SPWVD相結(jié)合,也取得了較好的效果.需要特別說明的是,與SPWVD結(jié)合后,盡管犧牲了一定的二維時(shí)頻聚集度,但大幅度提高了非平穩(wěn)信號三維時(shí)頻表示的準(zhǔn)確度.BGabor-NSPWVD算法流程與BGabor-NWVD算法類似,將信號的WVD換成SPWVD即可,BGabor-NWVD變換算法流程框圖如圖2所示.

3 數(shù)值仿真實(shí)驗(yàn)

3.1 時(shí)頻分布對比

本節(jié)采用2個構(gòu)造函數(shù)來進(jìn)行多個時(shí)頻分析方法的測試實(shí)驗(yàn),重點(diǎn)討論每種方法的時(shí)間和頻率聚集度,以及每個頻率分量的幅值分布是否正確.構(gòu)造的2個函數(shù)屬于易產(chǎn)生交叉項(xiàng)的多頻分量信號,如式(10)和式(11),分別為較為復(fù)雜的四分量線性調(diào)頻信號和三分量正弦頻率信號,數(shù)值實(shí)驗(yàn)的方法采用STFT、WVD、SPWVD、Gabor-WVD、Gabor-SPWVD、BGabor-NWVD和BGabor-NSPWVD.圖3為三分量和四分量構(gòu)造函數(shù)的二維時(shí)頻圖,采用本文改進(jìn)后的BGabor-NWVD和BGabor-NSPWVD,可以看出,這兩種方法都能準(zhǔn)確地顯示函數(shù)的時(shí)頻分布,沒有交叉項(xiàng)的干擾,且時(shí)間–頻率軸都有較好的聚集度.在進(jìn)行BGabor-NWVD和BGabor-NSPWVD處理時(shí),同時(shí)比較使用最小值法、二值化法和冪指數(shù)調(diào)節(jié)法實(shí)現(xiàn)Gabor-WVD和Gabor-SPWVD的結(jié)果,討論各方法的三維幅值分布是否準(zhǔn)確.線性調(diào)頻四頻率分量信號表示如式(10)所示,m為一般參數(shù),可以根據(jù)需要調(diào)整,這里選取m=0.2.

正弦三頻率分量信號如式(11)所示.

圖2 BGabor-NWVD算法流程圖Fig.2 BGabor-NWVD algorithm flow chart

圖3中,BGabor-NWVD 和BGabor-NSPWVD方法都能較好地實(shí)現(xiàn)復(fù)雜信號的二維時(shí)頻分布,消除了交叉項(xiàng)的干擾,BGabor-NWVD較BGabor-NSPWVD的二維時(shí)頻聚集度略好,采用NSPWVD方法犧牲了少許的平面聚集性,但這并不影響時(shí)頻表現(xiàn)和分布的判斷.下面討論這些方法的時(shí)頻三維分布情況,首先看基本的Gabor、WVD及WVD的改進(jìn)SPWVD,如圖4所示.

圖3 四分量f2和三分量f3的二維時(shí)頻圖Fig.3 Two-dimensional time-frequency diagram of four components signalf2and three components signalf3

圖4為函數(shù)f2的Gabor、WVD和SPWVD方法得到三維時(shí)頻分布圖,其中,圖4(a)中的Gabor方法沒有交叉項(xiàng),但是底端呈擴(kuò)散狀分布,頻率聚集性差;圖4(b)中WVD的交叉項(xiàng)干擾嚴(yán)重;圖4(c)中SPWVD沒有交叉項(xiàng),但是頻率項(xiàng)上的幅度有尖端失真,需要予以消除.

3.2 BGabor-NWVD

針對圖4出現(xiàn)的問題,采用二值化改進(jìn)的Gabor與WVD、SPWVD相結(jié)合的算法,時(shí)頻分析三維效果如圖5所示.

將Gabor與WVD直接相結(jié)合,按照式(7)～(9)分別進(jìn)行Gabor-WVD最小值法、二值化法和冪系數(shù)調(diào)節(jié)法的三維時(shí)頻分布畫圖,如圖5(a)、5(b)和5(c)所示,3個圖中函數(shù)f2和f3的頻率分量幅度出現(xiàn)了較大失真.由此需要對Gabor-WVD算法進(jìn)一步改進(jìn),先將Gabor二值化(BGabor),再結(jié)合歸一化的WVD(NWVD),得到BGabor-NWVD算法.實(shí)驗(yàn)仿真函數(shù)f2和f3的三維時(shí)頻分布如圖5(d),較好地克服了Gabor-WVD方法的幅值失真問題,其幅值頂端依舊存在模糊現(xiàn)象,還需要進(jìn)一步加以改進(jìn).

圖4 四分量f2的三維時(shí)頻圖Fig.4 Three-dimensional time-frequency diagram of four components signalf2

圖5 基于Gabor和WVD的四分量f2(上)和三分量f3(下)的三維時(shí)頻比較圖Fig.5 Three-dimensional time-frequency diagram of four-componentsf2(upper)and three-componentsf3(bottom)based on Gabor and WVD

3.3 BGabor-NSPWVD

SPWVD具有較好地平滑效果,考慮將Gabor與SPWVD相結(jié)合,可以進(jìn)一步優(yōu)化幅值模糊現(xiàn)象,如圖6所示.

Gabor與SPWVD相結(jié)合,可以很大程度上平滑信號的幅度數(shù)值,使其時(shí)頻表現(xiàn)更接近真實(shí)的各頻率分量幅度.但圖6中,圖6(a)、6(b)和6(c)的幅度仍然存在一些干擾和噪聲分量,其中,圖6(a)幅度頂端較好,但是幅度的底端范圍較寬,出現(xiàn)能量泄漏現(xiàn)象.圖6(d)的BGabor-NSPWVD算法效果最佳,幅度刻畫良好,仿真函數(shù)f2和f3各分量的頂端和低端沒有出現(xiàn)明顯失真,時(shí)間–頻率–幅度分布清晰準(zhǔn)確.BGabor-NSPWVD和Gabor-NWVD算法都能去除信息項(xiàng)上的疊加殘余,對交叉項(xiàng)有較好地抑制效果.兩種算法相比較,BGabor-NSPWVD對信號分量在時(shí)頻平面的幅度還原效果更優(yōu)于Gabor-NWVD,其時(shí)間–頻率–幅度的三維分布具有良好的表現(xiàn)能力.

4 魯棒性實(shí)驗(yàn)

4.1 二維時(shí)頻分布對比

BGabor-NWVD和BGabor-NSPWVD算法不僅對復(fù)雜信號的交叉項(xiàng)有較理想的消除作用,而且對附著在信號上的高斯白噪聲同樣有較強(qiáng)的抑制效果,兩種算法在時(shí)頻分析處理效果上具有強(qiáng)魯棒性的特點(diǎn),利于處理實(shí)際的含噪信號.構(gòu)造具有多頻分量的實(shí)驗(yàn)仿真信號,加入?10dB至20dB的高斯白噪聲,對比Gabor、WVD、SPWVD、Gabor-WVD、Gabor-SPWVD、BGabor-NWVD和BGabor-NSPWVD方法的二維時(shí)頻分析效果.由于篇幅限制,本文僅列出加入2dB噪聲時(shí)函數(shù)f4的二維時(shí)頻分布,如圖7所示.

圖6 基于Gabor和SPWVD的四分量f2(上)和三分量f3(下)的三維時(shí)頻比較圖Fig.6 Three-dimensional time-frequency diagram of four-componentsf2(upper)and three-componentsf3(bottom)based on Gabor and SPWVD

圖7 含噪信號f4的二維時(shí)頻分布比較(SNR=2dB)Fig.7 The two-dimensional time-frequency distribution of the noisy signalf4(SNR=2dB)

圖7(b)為改進(jìn)前的Gabor方法,即高斯變換核STFT,時(shí)頻聚集度較差且噪聲分布干擾嚴(yán)重.圖7(c)的WVD方法存在嚴(yán)重的交叉項(xiàng)干擾,需要去除存在于真實(shí)頻率之間的虛假頻率分量.圖7(d)的SPWVD較好地去除了交叉項(xiàng)和噪聲干擾,但是由于采用偽平滑的手段,使得該方法的時(shí)頻聚集度較差.圖7(e)、7(f)和7(g)是結(jié)合后的Gabor-WVD方法,其中冪系數(shù)調(diào)節(jié)法去除交叉項(xiàng)的效果略好,但是這三種方法對噪聲的抑制和去除交叉項(xiàng)都不滿足要求.圖7(h)采用進(jìn)一步結(jié)合后的BGabor-NWVD方法,較好地去除了交叉項(xiàng)和噪聲干擾,具有最佳的二維時(shí)頻分布表現(xiàn).

圖7(i)、7(j)和 7(k)是 Gabor與 SPWVD相結(jié)合的三種方法,對交叉項(xiàng)的抑制尚可,但對噪聲干擾的去除效果還需提高.圖7(l)為BGabor-NSPWVD方法,對交叉項(xiàng)的抑制和噪聲的去除都較為理想,時(shí)頻聚集度略遜于圖7(h)的BGabor-NWVD方法,但是對交叉項(xiàng)的抑制和魯棒性而言,BGabor-NSPWVD和BGabor-NWVD都是可以采用的方法.

4.2 時(shí)頻聚集度評價(jià)

文獻(xiàn)[20]給出了評價(jià)時(shí)頻聚集度的量化公式,如式(12),n為時(shí)間窗長度,ω為頻率,式(12)通過調(diào)整窗口參數(shù),求得最大的時(shí)頻分布的第四冪范數(shù)與第二冪范數(shù)之商來評價(jià)時(shí)頻聚集度.向?qū)嶒?yàn)函數(shù)f4添加?10dB至20dB的高斯白噪聲,求解函數(shù)f4的EJP數(shù)值,用以衡量各個方法的聚集度分布,如表1所示.將表1的數(shù)值用圖8的折線趨勢來表示,能夠較為直觀地比較每種方法的聚集度效果.

時(shí)頻聚集度評價(jià)圖8中,信噪比(SNR)大于0dB時(shí),BGabor-NWVD和BGabor-NSPWVD方法的聚集度EJP開始平穩(wěn),且其值高于其他方法,隨著SNR的變化,這兩種改進(jìn)算法皆具有較高的時(shí)頻聚集度和強(qiáng)魯棒性,BGabor-NWVD的聚集度優(yōu)于BGabor-NSPWVD算法,與圖7的二維時(shí)頻分布效果圖一致.

5 硬質(zhì)合金頂錘破裂信號監(jiān)測

圖8 時(shí)頻聚集度參數(shù)EJP評價(jià)比較Fig.8 Comparison of time-frequency aggregation degree evaluation onEJP

表1 仿真函數(shù)f4在不同噪聲條件下各方法的聚集度EJP數(shù)值比較Table 1 TheEJPnumerical comparison of experimental functionf4in different noise conditions

圖9 硬質(zhì)合金頂錘的現(xiàn)場實(shí)物圖與三維模擬圖Fig.9 Carbide anvil physical site map and 3D simulation figure

在人造金剛石合成加工過程中,硬質(zhì)合金頂錘在交變的外載應(yīng)力及熱應(yīng)力作用下,易發(fā)生彈性形變和塑性形變.如由于疲勞損壞而產(chǎn)生微觀裂紋,這些裂紋不斷發(fā)展貫穿許多晶粒成為宏觀裂紋,使頂錘斷面進(jìn)而發(fā)生橫向或縱向的破裂或壓潰[22].作為高脆性材料,受材料自身特點(diǎn)的正常性破壞和結(jié)構(gòu)、工藝、人為操作等非正常性破壞等因素影響,硬質(zhì)合金頂錘對于微觀裂紋極為敏感,當(dāng)產(chǎn)生微觀裂紋時(shí),由裂紋成核、裂紋擴(kuò)展延伸并迅速發(fā)展成宏觀裂紋,產(chǎn)生強(qiáng)烈的聲裂發(fā)射頻率段信號[23].大量實(shí)驗(yàn)研究表明,聲發(fā)射源主要有塑性形變(滑移和孿生)、斷裂(裂紋的形成和擴(kuò)展、第二相質(zhì)點(diǎn)或夾雜物)、相變(馬氏體相變、共晶反映等)、磁效應(yīng)和表面效應(yīng)等[22].圖9(a)為硬質(zhì)合金頂錘的工作現(xiàn)場圖,圖9(b)為硬質(zhì)合金頂錘的三維模擬效果圖.

通過檢測硬質(zhì)合金頂錘工作時(shí)的疑似破裂信號,對比正常狀態(tài)下信號頻率特性及裂紋產(chǎn)生時(shí)的信號頻率特征,判別頂錘破裂的發(fā)生與否.目前存在的難點(diǎn)是微小破裂過程中釋放的應(yīng)力能太小,聲發(fā)射信號相當(dāng)微弱,在頂錘工作惡劣環(huán)境下,受復(fù)雜噪聲信號干擾,難以有效地設(shè)置采集頻率窗口和區(qū)分破裂頻率成分.為有效地接收到破裂信號,需要傳感器設(shè)置合適的頻率檢測閾值,可采用時(shí)頻分析的手段來分析疑似破裂信號,為判斷是否破裂提供參考依據(jù).本文利用BGabor-NWVD和BGabor-NSPWVD算法對金屬破裂疑似樣本信號進(jìn)行分析,獲取其時(shí)頻聯(lián)合分布特點(diǎn),找到疑似金屬破裂信號段的時(shí)頻分布表示.

圖10中,圖10(a)、10(d)、10(g)為疑似破裂樣本信號的時(shí)域波形,圖10(b)、10(e)、10(h)為采用BGabor-NWVD算法得到的時(shí)頻分析結(jié)果,圖10(c)、10(f)、10(i)為采用BGabor-NSPWVD 算法得到的時(shí)頻分析結(jié)果.圖10(c)、10(f)、10(i)的時(shí)頻集中區(qū)域較圖10(b)、10(e)、10(h)更為明顯,在疑似破裂發(fā)生時(shí)段,3個樣本信號的頻率分布集中范圍,大致分布在100kHz這個較為明顯的頻率區(qū)域.根據(jù)BGabor-NSPWVD算法得到的金屬破裂頻率窗口,將其設(shè)置為傳感器的破裂判斷閾值,用以自動在線監(jiān)測硬質(zhì)合金頂錘工作時(shí)是否發(fā)生疑似破裂,提高了判別幾率,取得了較好的效果.

6 結(jié)論

針對STFT、Gabor和WVD出現(xiàn)的時(shí)頻分辨率模糊和存在交叉項(xiàng)等缺點(diǎn),以及一些結(jié)合算法如STFT-WVD和Gabor-WVD出現(xiàn)的三維幅度失真,抗噪性能及魯棒性能還不理想的問題,本文提出BGabor-NWVD和BGabor-NSPWVD算法.通過對復(fù)雜線性調(diào)頻信號和多分量的正弦頻率分量信號進(jìn)行數(shù)值仿真實(shí)驗(yàn),BGabor-NWVD和BGabor-NSPWVD算法在抑制了交叉項(xiàng)的同時(shí),具有較高銳化時(shí)頻分辨率,兩種算法的抗噪性能和魯棒性也較為理想.由四分量線性調(diào)頻信號和三分量正弦信號的數(shù)值仿真實(shí)驗(yàn)可知,BGabor-NWVD的二維時(shí)頻表示優(yōu)于BGabor-NSPWVD,BGabor-NSPWVD的三維時(shí)頻表示優(yōu)于BGabor-NWVD.BGabor-NWVD和BGabor-NSPWVD算法綜合了STFT、Gabor、WVD和SPWVD各自的頻率自適應(yīng)性和良好的時(shí)頻表示,其時(shí)頻聚集度評價(jià)參數(shù)EJP高于其他方法,具有高銳化頻率聚集度優(yōu)點(diǎn)的同時(shí),能夠真實(shí)還原信號頻率分量的幅度.通過硬質(zhì)合金頂錘工作時(shí)產(chǎn)生的疑似破裂樣本信號進(jìn)行時(shí)頻分析,本文方法可以較為準(zhǔn)確地尋找傳感器的頻率判別窗口,為金屬破裂監(jiān)測設(shè)備數(shù)據(jù)采集卡提供有效的閾值參考.

圖10 疑似金屬破裂樣本的時(shí)頻分析Fig.10 Time-frequency analysis of suspected metal rupture samples