基于NGWarblet-WVD 的高質量時頻分析方法

郝國成 , 馮思權 王巍 凌斯奇 譚淞元

時頻分析方法采用時間和頻率域聯合函數來處理時序信號,獲取信號頻率隨時間變化的細節信息,是現代信號處理領域的重要技術手段之一.1946 年,Gabor 首次對傅里葉變換加高斯窗函數,提出了著名的Gabor 變換,由此開啟了時頻聯合分析的新思路[1].Ville 將量子力學的Wigner 分布用于信號分析與處理領域,提出了Wigner-Ville 分布(Wigner-Ville distribution,WVD)[2].WVD 具有良好的時頻聚集特性,但存在嚴重的交叉項干擾問題.Cohen對WVD 進行時頻二維卷積得到Cohen 類時頻分布[3].Cohen 類時頻分布通過構造核函數,達到消除或抑制交叉項的目的,缺點是降低了頻率聚集度.20 世紀80 年代,Mallat[4]提出了多尺度分析思想和Mallat 算法,成功地統一了各種小波函數的構造模型,使用可調節的時間和頻率窗口,有效提高了頻率聚集度.1996 年,美國地球物理學家Stockwell等[5]對短時傅里葉變換和連續小波變換的思想進行延伸與推廣,提出了S 變換.1998 年,Huang等[6]提出經驗模態分解(Empirical mode decomposition,EMD),將信號分解為有限個固有模態函數(Intrinsic mode functions,IMF)的集合.EMD 方法可應用于多種類型信號的分解,在處理非平穩、非線性信號方面,效果良好[7?8].但相對來講,S 變換和EMD的聚集度不夠理想.2012 年,Yang等[9]提出了適用范圍更廣泛的參數化時頻分析廣義Warblet變換(Generalized Warblet transform,GWT),該方法具有真實反映信號頻率分布的特點,但存在頻率泄露現象,時頻聚集度較差,需要結合其他算法加以改進.

WVD 作為一種優良的雙線性時頻分析方法具有其他方法不可替代的高銳化時頻聚集特性,但在處理多分量信號時會出現交叉項干擾問題.文獻[10?12]提出對信號進行WVD 處理前,分別采用變分模態分解 (Variational mode decomposition,VMD)、集成經驗模態分解 (Ensemble empirical mode decomposition,EEMD)、自適應匹配追蹤 (Adaptive matching pursuit,AMP)將多分量信號轉換為單分量信號,有效避免產生交叉項;文獻[13]通過對信號進行帶通濾波和相位矯正的方法去除交叉項;文獻[14]提出矩陣旋轉變換的方法,將WVD 交叉項旋轉至與頻率軸平行,再通過濾波器濾除交叉項;文獻[15]通過二重余弦信號的WVD,推導自項與交叉項位置關系以及振蕩特性,濾除交叉項;文獻[16?17]采用兩種算法結合的思想(如Gabor-WVD、BGabor-NSPWVD、SPWVD-WVD 等),對比實驗證明此方法有效抑制了WVD 交叉項.

本文在分析WVD 交叉項的產生原因基礎上,將交叉項分為新產生的交叉項分量和混入自項的交叉項分量兩種類型.利用GWT 較好還原信號真實頻率分布的特性,將GWT 矩陣與WVD 矩陣聯合處理以實現濾波效應,抑制WVD 的兩種類型交叉項.使用兩種定量評價方法將NGWT-WVD 算法與同類算法進行對比,檢驗算法有效性.最后將該算法用于處理金屬破裂樣本信號,獲取破裂期間的時頻分布圖,找出濾波器的窗口門限近似頻率,為聲發射信號監測傳感器采集卡提供門限設置依據.

1 NGWT-WVD 算法

1.1 WVD (Wigner-Ville 分布)

WVD 是一種雙線性時頻分布,設信號為z(t),其Wigner-Ville 分布表達式為

R(t,τ)=z(t+τ/2)z?(t ?τ/2)表示信號的自相關函數,WVD 可看作信號自相關函數的傅里葉變換.WVD 具有優良的時頻聚集度,但處理多分量信號時會出現交叉項,如式(2)所示

Ra(t,τ)表 示自相關成分(自項),是有用信息;Rs(t,τ)表示互相關成分(交叉項),是干擾信息.對于同時間區間的多分量信號,交叉項具有兩個特點: 任意兩個信號分量均會產生一個交叉項;交叉項位于兩頻率分量中間.

選取三分量線性調頻信號z1(t) 為例,說明信號各成分之間的關系.其表達式如式(3)所示

信號二維、三維WVD 時頻分布如圖1 所示.由圖1(a)和圖1(b)對比可見,三分量線性調頻信號的Wigner-Ville分布引入了一些新分量,這些分量是新產生的交叉項分量;而由圖1(c)的三維圖可見,中間分量能量遠高于兩邊分量,說明有部分交叉項混入進自項成分.本文需要解決的問題有兩個: 一是抑制新產生的交叉項分量;二是抑制混入自項的交叉項分量.

圖1 三分量信號的WVD 時頻圖Fig.1 Time-frequency diagram of WVD of three-components signal

1.2 GWT (廣義Warblet 變換)

GWT 是核函數以傅里葉級數為模型定義的參數化時頻分析方法.其定義為

其中,t0∈R 表示時間窗滑動時的窗中心所在時間,wσ∈L2(R)定義了一個非負對稱的標準化實窗,通常是高斯窗,(t) 計算式為

參數化時頻分析方法選取合適的核函數對其分析效果有很大影響,采用傅里葉級數變換核的GWT能夠分析具有周期性或非周期性時頻特征的非平穩信號,以及具有強振蕩時頻特征的信號,使其適用范圍更加廣闊.

三分量線性調頻信號的GWT 時頻圖如圖2 所示.

圖2 中GWT 算法采用傅里葉級數模型作為核函數,較好地還原了信號的真實頻率分布,雖然時頻聚集度較差,但能夠保留真實頻率分量,不會產生交叉項.由此,基于WVD 的高銳化特性和GWT真實還原信號時頻分布的特性,將二者相結合,得到GWT-WVD 算法,既能抑制交叉項,又能保留高銳化時頻聚集度的性能.

圖2 三分量信號的GWT 時頻圖Fig.2 Time-frequency diagram of GWT of three-components signal

1.3 GWT-WVD 算法: 抑制新產生的交叉項分量

本文提出的GWT-WVD 算法,可有效抑制WVD 的交叉項分量.該算法通過對GWT 算法與WVD 算法得到的矩陣進行運算,在保留較好的時頻聚集度的同時,能較好地消除或抑制交叉項.其表達式為

其中,GWarx(t,f)與W x(t,f) 分別代表GWT與WVD,p(x,y) 為聯合處理函數.本文采用了3種不同的函數,得到GWT-WVD 算法的3種定義式

式(7)～(9)分別采用最小值法、二值化法、冪指數調節法對兩矩陣進行處理,其各自實現思想如下.

1)最小值法.比較WVD 矩陣及GWT 矩陣對應位置元素,篩選其中較小的元素,按比例處理后組成GWT-WVD 矩陣,使WVD 矩陣新產生的交叉項所在位置的元素被GWT 矩陣元素取代.

2)冪指數調節法.調節冪指數,增強GWT與WVD矩陣中信號數據對應的元素,削弱交叉項數據對應的元素,再將兩矩陣點乘,得到GWT-WVD矩陣.

3)二值化法.選取合適的閾值將GWT 矩陣二值化得到新矩陣,用新矩陣點乘WVD 矩陣得到GWTWVD 矩陣.WVD 矩陣有交叉項的元素位置對應的GWT 矩陣元素會小于閾值,二值化后為0,與WVD 矩陣相乘后可消除新產生的交叉項.

線性調頻信號z1(t) 的3種GWT-WVD 時頻分布如圖3 所示.3種結合方法均能優化信號的時頻分布,有效抑制交叉項,但冪指數調節法的時頻聚集度比另外兩種方法差(第3 節討論時頻聚集度問題).在處理復雜信號時,二值化法閾值難以確定.綜合比較,最小值法的GWT-WVD 算法性能最佳.

圖3 3種GWT-WVD 時頻圖Fig.3 Time-frequency diagram of three types of GWT-WVD

1.4 NGWT-WVD 算法: 抑制混入自項的交叉項分量

GWT-WVD 算法能夠有效抑制新產生的交叉項分量,但無法抑制混入自項分量的交叉項.針對這個問題,本文對GWT-WVD 算法進行改進,又提出NGWT-WVD 算法,其主要實現思路如下.

1)采用廣義Warblet 變換和WVD 分別對原信號進行處理,得到GWT 矩陣和WVD 矩陣;

2) 找出G W T 矩陣中元素數值的最大值GWTmax,并記錄其所對應的位置 (i,j),將GWT矩陣中的各元素除以GWTmax,即對GWT 矩陣進行歸一化,得到矩陣GWT-1;

3) 記錄GWT-1 矩陣中元素數值的最小值GWTmin,最小值GWTmin 要求非零,并用GWTmin的數值替換掉矩陣GWT-1 中所有值為0 的元素;

4)找出WVD 矩陣中位置為 (i,j) 的元素,將其記為WVDmax,同時將WVD 矩陣中的各個元素除以WVDmax,得到矩陣WVD-1;

5)用矩陣WVD-1 點除矩陣GWT-1,得到矩陣T,選取矩陣T中大于x的元素以及小于y的元素,x設置為5,y設置為2,將所對應的元素位置置1,并記錄大于x的元素位置;

6)在WVD 矩陣中找出與上一步記錄對應的元素位置,并將此位置上元素置為0,最后用WVD矩陣點除矩陣T,輸出NGWT-WVD 矩陣.其流程如圖4 所示.

圖4 NGWT-WVD 算法流程圖Fig.4 Algorithm flowchart of NGWT-WVD

線性調頻信號的GWT-WVD 算法和NGWTWVD 算法三維時頻分布如圖5 所示.由圖5(a)可見,中間分量的能量譜明顯高于其他分量的能量譜,GWT-WVD 算法不能有效抑制混入自項的交叉項分量;而圖5(b)中3 個信號分量能量一致,混入自項成分上的交叉項被有效抑制.圖5 說明NGWTWVD 算法能有效抑制線性調頻信號的兩類交叉項.

圖5 三分量信號三維時頻圖比較Fig.5 Three-dimensional time-frequency diagrams comparison of three-component signals

NGWT-WVD 算法中的閾值參數x是為了抑制存在于自項中的干擾項,閾值y是為了抑制新產生的干擾項并且去除發散的能量.閾值對NGWTWVD 算法性能起決定性作用,這里使用后文的兩種定量評價方式對兩個閾值的敏感性進行分析,其結果如圖6 所示.

圖6(a)是CM值(時頻聚集度定量評價方法,數值越大代表聚集度越高,其計算式見第3.2 節)隨閾值y變化折線圖,此時閾值x設置為5,x1,x2,x3,x4 分別代表文中使用的4種仿真信號.圖中,閾值y小于1 時,NGWT-WVD 算法的時頻聚集度低于WVD 算法;隨著閾值y的增加,CM值逐漸增加,WVD 中部分擴散的時頻系數被濾除,時頻聚集度提高.當閾值y超過2 時,時頻聚集度CM值迅速提高,但這是以濾除部分WVD 有用信息分量為代價.因此閾值y取值為2 即可,此時新產生的交叉項基本消除干凈(由圖5(b)可見),且時頻聚集度較理想.圖6(b)是時變功率譜誤差隨閾值x變化折線圖,此時閾值y設置為2 (基本消除了新產生的交叉項分量).閾值x是為濾除混入自項的交叉項(只有仿真信號1 含有混入自項的交叉項),當閾值x大于10 時,時變功率譜誤差保持在25%左右,無法濾除混入自項的交叉項;閾值x在5～6 之間時,時變功率譜誤差接近0,此時抑制效果較好;小于5 后,由于濾除了信息項,導致誤差進一步增加,因此閾值x設置為5 較合適.實際應用信號由一系列單分量信號線性疊加而成,也可按照此參數設置相應閾值.

圖6 閾值敏感性測試圖Fig.6 The test chart of threshold sensitivity

NGWT-WVD 算法通過對GWT 矩陣進行歸一化且做去零處理得到的GWT-1 矩陣,此矩陣為對照矩陣,通過設置兩個閾值,實現了濾波效應,剔除了WVD 中發散能量,且抑制了混入自項的交叉項和新產生的交叉項,得到了更加理想的時頻分布結果.

2 數值實驗

為了評價NGWT-WVD 算法的處理效果,本文另外構造了3種具有代表性的信號,并通過與Garbor-WVD、VMD-WVD 兩種較有代表性的方法進行對比,驗證該算法的有效性.

2.1 分段信號

函數表達式為

分段信號z2(t) 的WVD、Gabor-WVD (文獻[14]中采用的處理方法)、VMD-WVD (選取模數為4)、NGWT-WVD 算法時頻分布如圖7 所示.圖7(a)中WVD 在端點處出現嚴重的交叉項干擾;圖7(b)中Gabor-WVD 算法處理信號時仍然會出現少許交叉項;圖7(c) 中,VMD 的參數設置準確的情況下,VMD-WVD 算法交叉項抑制效果較好;圖7(d)中NGWT-WVD 算法得到的時頻分布也基本沒有交叉項.

圖7 分段信號時頻圖Fig.7 Time-frequency diagram of segmented signal

2.2 交叉型信號

函數表達式為

交叉信號z3(t) 的WVD、Gabor-WVD、VMDWVD (選取模數為2)以及NGWT-WVD 算法的時頻分布如圖8 所示.此信號為兩個線性調頻信號交叉于一點,交叉點的頻率分量混合在一起,分離困難.圖8(b)中Gabor-WVD 算法在處理該信號時,不能十分有效地抑制交叉項,且時頻聚集度較低.圖8(c)中VMD-WVD 算法交叉項抑制效果較好,時頻聚集度略差;圖8(d)中NGWT-WVD 算法對此信號的頻率分量刻畫準確,時頻聚集度高.

圖8 交叉型信號時頻圖Fig.8 Time-frequency diagram of cross-type signal

2.3 兩分量調頻信號

函數表達式為

信號z4(t) 是由一個拋物線調頻信號和一個正弦調頻信號組成,并存在交叉點.該信號的WVD、Gabor-WVD、VMD-WVD (選取模數為2)以及NGWT-WVD 算法的時頻分布如圖9 所示.由圖可見,圖9(b)中Garbor-WVD 不能完全抑制,圖9(c)中VMD-WVD 在處理這一類復雜調頻信號時,即使參數設置準確,也無法有效將單分量信號完全分開,交叉項干擾較為嚴重;而NGWT-WVD 算法對交叉項具有較好的抑制效果.由此,以上3 個仿真實驗都能說明NGWT-WVD 算法具有較好的抑制信號交叉項效果,同時保留了高時頻聚集度,提高了時頻分析質量.

圖9 兩分量調頻信號時頻圖Fig.9 Time-frequency diagram of two-component frequency modulated signal

3 算法性能評價

為定量說明NGWT?WVD 算法性能,本文選取交叉項抑制效果評價和時頻聚集度評價兩個指標評判該算法的時頻分析效果.

3.1 交叉項抑制效果評價

交叉項抑制效果是WVD 改進算法性能評價的重要指標,以往均使用定性分析,即通過觀察時頻分布結果得出判斷.本文在分析WVD 交叉項出現規律的基礎上,提出一種定量的交叉項抑制效果評價方法.由式(2)知,信號經WVD 處理后,分量包括自項成分與交叉項成分,本文計算信號在時間?頻率平面的時變功率譜Pi,并與信號標準時變功率譜IPi(不含交叉項的信號功率譜)作對比,計算出時變功率譜誤差,以此評價交叉項抑制效果.選取信號z1(t),z2(t),z3(t),z4(t) 作為研究對象,計算各算法的功率譜與標準功率譜的平均相對誤差,其計算式為

其中,N為信號采樣點數.標準時變功率譜計算式為

其中,n表示信號含有單分量的個數,Pfj(i) 代表第j個單分量的時變功率譜.將各單分量分別求時頻分布后相加,得到沒有交叉項的標準時變功率譜.四種仿真信號平均時變功率譜誤差對應的柱形圖如圖10 所示,數據如表1 所示,z1,z2,z3,z4 分別代表信號z1(t),z2(t),z3(t),z4(t).

圖10 各算法的時變功率譜誤差柱形圖Fig.10 Time-varying power spectrum error column chart of each algorithm

表1 各算法的時變功率譜誤差比較Table 1 Time-varying power spectrum error comparison of each algorithm

表1 中不同信號的NGWT-WVD 算法的時變功率譜誤差均為最小,反映出算法抑制交叉項性能最佳.圖11 為表1 的平均時變功率譜誤差折線圖形式.

圖11 各算法的時變功率譜誤差折線圖Fig.11 Time-varying power spectrum error line chart of each algorithm

圖11 中WVD 算法的誤差值最大,說明時頻分布中含有大量的交叉項,Gabor-WVD、GWTWVD 算法都能一定程度上抑制交叉項,但是無法有效抑制混入自項成分中的交叉項,其時變功率譜誤差仍很大;VMD-WVD 算法會先分解信號,在處理恒頻信號時,分解效果好,交叉項抑制效果也好,在處理復雜信號時,分解效果不佳,時變功率譜誤差遠高于NGWT-WVD 算法;而NGWT-WVD 算法可以較好地抑制兩種交叉項,時變功率譜誤差較小,且不需要設置初始參數,算法的適應性強.

3.2 時頻聚集度評價

Shafi等[18]于2009 年提出了關于時頻聚集度量化評定的方法,本文選取CM值作為時頻聚集度量化標準,評價信號的時頻圖的聚集度,其計算式為

其中,n為時間窗長度,ω為信號在某點的頻率,Q表示信號的時頻分布.選取上述4種仿真信號作為實驗對象,計算CM值.圖12 為4種信號各時頻分析方法的CM值柱形圖,其數據如表2 所示.

圖12 各算法的CM 值柱形圖Fig.12 CM value column chart of each algorithm

表2 中NGWT-WVD 算法在評價不同信號的時頻聚集度時,其CM值均為最大,反映出算法時頻聚集度性能最佳,銳化程度最高.圖13 時頻聚集度CM值折線圖.

表2 各算法的CM 值比較(×10?3)Table 2 CM value comparison of each algorithm(×10?3)

圖13 中,GWT 算法CM值最小,時頻聚集度最差.Gabor-WVD 算法CM值低于GWT-WVD算法,時頻聚集度較GWT-WVD 算法低,其去除交叉項的同時,降低了時頻聚集度.信號z1 因為是一種線性調頻信號,其VMD-WVD、NGWT-WVD算法的CM值與WVD的CM值相差不大,但對于另外3種信號,VMD-WVD、NGWT-WVD 算法的時頻聚集度CM值明顯高于WVD 算法,且NGWT-WVD 算法略優于VMD-WVD 算法.通過比較4種信號各時頻分析方法的時頻聚集度CM值大小,驗證了NGWT-WVD 算法具有高銳化的時頻聚集度.

圖13 各算法的CM 值折線圖(×10?3)Fig.13 CM value line chart of each algorithm (×10?3)

4 金屬破裂檢測信號的時頻分析應用

人造金剛石合成過程中,六面頂壓機頂錘的破裂損壞是經常發生的生產事故.頂錘采用的鎢鈷類硬質合金,在高溫高壓環境中,長期處于超臨界應力狀態,會出現疲勞損傷,進而而導致頂錘破裂[19].若金剛石生產加工過程中未發現頂錘破裂前的異常情況,將會使六面頂壓機頂錘出現不可逆性損壞,造成嚴重生產損失.六面頂壓機和頂錘結構如圖14所示.

圖14 六面頂壓機和硬質合金頂錘Fig.14 Cubic press and carbide anvil

聲發射技術是一種有效的探傷檢測手段.構件在外力或應變力的作用下會激發一定頻譜的聲發射信號,通過判斷接收到的信號頻譜強度來預判構件的缺陷嚴重程度[20].頂錘破裂大致分為: 裂紋成核、裂紋拓展、斷裂3 個過程,在這3 個過程中應變能以彈性應力波的形式釋放出來,會產生劇烈的聲發射信號[20].

為了避免生產事故,需要在初步檢測到鎢鈷合金出現破裂時立即終止加工過程.而采集到的聲發射信號振幅微弱(如圖15(a)所示),難于甄別判斷.本文使用的疑似金屬破裂信號數據采樣頻率為40 000 Hz,選取2 000 個采樣點數據,由于篇幅限制,僅繪制時域圖、WVD 算法時頻分布圖以及NGWT-WVD算法時頻分布圖,結果如圖15 所示.兩種時頻分析方法的CM值對比如表3 所示,NGWT-WVD 的時頻聚集度CM值最大,時頻聚集度最好.

圖15 疑似金屬破裂樣本時頻分析Fig.15 Time-frequency analysis of suspected metal rupture samples

表3 六種算法的CM 值比較(×10?5)Table 3 CM value comparison of six algorithms(×10?5)

圖15(a) 為某一疑似金屬破裂樣本信號時域圖,該信號振幅微弱,不利于監測傳感器報警閾值的設置;圖15(b)和圖15(c)分別為該樣本信號WVD算法以及NGWT-WVD 算法處理得到的時頻分布圖.為便于更精確的分析,截取了頻率發生劇烈波動的片段(采樣點數為900～1 300 之間數據)進行分析,結果如圖15(d)和圖15(e)所示.圖15(d)中WVD算法得到的時頻分布圖交叉項分量與信號分量混雜,整個時頻分布圖雜亂不清,難以確定出現金屬破裂的精確時間節點,無法有效示警.圖15(e)中可較為明顯地看出抑制了圖15(d) 中的交叉項(尤其是圖15(d)方框中的主要交叉項),時頻分布銳化聚集度有明顯的提高,CM值驗證了這一結果.圖15(e)中橢圓框標記了金屬破裂過程中頻率較集中的頻率分量區域,可作為監測濾波器組的通帶上下限的選取范圍.其中濾波器通帶2 作為主要預警通帶,濾波器通帶1 由于頻率較低,較容易受到外界噪聲干擾,通帶3 則由于頻率閾值過高,容易遺漏低頻預警信號,因此濾波器通帶1、3 通常作為輔助預警通帶,此外還可以在各通帶間再設置部分濾波器通帶,提高識別幾率.實驗結果表明,NGWTWVD 算法能夠較精確地顯示出各信號出現破裂的具體時間和頻率窗口值,可為信號監測傳感器和濾波器組提供可操作的判斷閾值,提高設備的監測成功率.

5 結束語

本文分析了WVD 產生交叉項的原理,針對交叉項干擾和時頻模糊問題,提出了NGWT-WVD算法.該算法不僅能夠有效抑制新產生的交叉項分量,而且解決了Gabor-WVD 等算法無法消除混入自項成分的交叉項分量問題,在交叉項抑制效果評價和時頻聚集度評價中表現良好.仿真結果表明,NGWT-WVD 算法能夠實現保持高銳化聚集度的同時,有效抑制交叉項干擾,是一種高質量的時頻分析方法.將該算法用于處理金屬破裂樣本信號,能夠得到較為精確的信號時間和頻率窗口值,為監測傳感器報警閾值的設置和數據采集濾波器組的設計提供有效依據.