金 陽，郝志勇，鄭 旭，景國璽

(1.浙江大學(xué) 能源工程學(xué)系，杭州 310027;2.湖北汽車工業(yè)學(xué)院 汽車工程系，十堰 442002)

往復(fù)運動的內(nèi)燃機其確定性激勵與其轉(zhuǎn)速間存在固有的聯(lián)系，激勵頻率與轉(zhuǎn)頻成倍數(shù)關(guān)系。內(nèi)燃機在工作轉(zhuǎn)速范圍內(nèi)加減速類似于給它自身施加了一種天然的掃頻激勵，只不過每個時刻的激勵頻率不是一個，而是多個成倍數(shù)關(guān)系的頻率。與此過程同步測量的振動或噪聲信號中必然包含了系統(tǒng)的動力學(xué)特性和聲輻射特征，并且還是結(jié)構(gòu)設(shè)計者所關(guān)心的工作頻率內(nèi)的特性。因而對變速階段信號的特征提取具有工程應(yīng)用意義［1］。

變速階段的振聲信號是非穩(wěn)態(tài)的，在旋轉(zhuǎn)機械領(lǐng)域，由于其信號中具有上述的以轉(zhuǎn)速為基頻的諧波性質(zhì)，每個諧波稱為階次，因而對這類特定信號的分析也被稱為階次分析。已發(fā)展出數(shù)種相關(guān)的處理方法［2］，各有利裨及其適用場合。本文研究的是兩種不需要角度域重采樣的階次分析技術(shù)-STFT與Gabor階次跟蹤［2-5］-兩個方面的問題:

(1)STFT幅值譜及Gabor系數(shù)譜對階次分量的分隔能力與什么有關(guān)?能否找出一個工程上可用的定量的評判標(biāo)準(zhǔn)?

(2)這兩種技術(shù)在階次分量瞬時幅值估計上的誤差水平如何?

1 高斯窗STFT幅值譜的對階次分量的分隔能力與階次切片的幅值誤差

1.1 離散STFT幅值譜

與FT諧波分析一樣，對于離散信號s［i］，當(dāng)采用的窗函數(shù)為 γ［i］時的離散 STFT 幅值{Am，n}m，n∈Z也采用類似的幅值校正算法:

1.2 高斯窗

為了考查STFT幅值譜對階次分量的分隔能力與窗參數(shù)的定量關(guān)系，除了選用時頻集中性最好的高斯窗外，還將連續(xù)時間域高斯窗的時間標(biāo)準(zhǔn)差σi(s)(以下將簡稱為“窗寬”，它是時間信號持續(xù)時間的一種度量。)做為輸入?yún)?shù)來生成STFT數(shù)值計算中的離散窗。這樣做的好處是，我們研究的信號原本就是屬于連續(xù)時間域上的，在同一域內(nèi)便于發(fā)現(xiàn)變換結(jié)果與窗之間的聯(lián)系。

設(shè)采樣頻率為fs(Hz)，離散點數(shù)為L(偶數(shù))，窗序列的編號為i，則連續(xù)時間域中時間標(biāo)準(zhǔn)差為σt(s)的高斯窗函數(shù)離散后成為:

其中:

在數(shù)值計算中，L的選取應(yīng)保證σN足夠小。本研究的數(shù)值計算中，一般均保證σN≤0.1，這意味著在高斯窗的兩端點的值已降至低于窗的峰值的0.2%。連續(xù)時間域中時間標(biāo)準(zhǔn)差為σt(s)的高斯窗函數(shù)的頻域標(biāo)準(zhǔn)差(Hz)為:

1.3 仿真分析

構(gòu)造的仿真信號為:

這是一個具有5個階次分量并添加了白噪聲的的信號(信噪比34 dB)，轉(zhuǎn)速:n=60 t;第p階分量的瞬時幅值:Ap(t)=1;而瞬時頻率:fp(t)=p·t。對其采用不同σt的高斯窗進行離散STFT，生成峰值幅值譜以進行對比。圖1是其中的3個結(jié)果。

每個結(jié)果中的圖形區(qū)由上下兩幅圖組成:下方是時變單邊峰值幅值譜，橫坐標(biāo)為時間(s)，縱坐標(biāo)為頻率(Hz)，顏色值為幅值(‘pk’代表峰值之意)。為解釋上的方便，其中有標(biāo)記為“0”與“1”的輔助點及一些輔助線;上方是“階次切片”，是STFT單邊峰值幅值譜上某階次在不同時刻的理論瞬時頻率處的對應(yīng)值，橫坐標(biāo)為時間，縱坐標(biāo)為幅值。

圖1 信號S1在不同窗寬下的STFT單邊峰值幅值譜(L:窗的點數(shù);△M:以點數(shù)表示STFT的時間步長)Fig.1 STFT single-sided peak-amplitude spectra with different window width for signal S1(L:window length in point number，△M:time step in point number)

1.3.1 高斯窗STFT幅值譜的對階次分量的分隔能力

STFT是最早使用的非穩(wěn)態(tài)信號時-頻分析方法，它的思想是直接而有效的。但其固有的局限性是時頻分辨率不能同時提高。此仿真的目的不為了展示時間分辨率與頻率分辨率之間反向變化的關(guān)系，而是想揭示在什么條件下，從STFT幅值譜上可以將各分量分隔開來，在分量間沒有或較少有干擾現(xiàn)象。

信號分量的理論分布在時頻面上越靠近的區(qū)域，在STFT頻譜圖上相應(yīng)的區(qū)間，兩分量間越易產(chǎn)生干擾模糊現(xiàn)象從而使分量難以分辨。變速階段信號的特點是不僅同一時刻存在多個階次分量，同一頻率下也有多個分量。

為便于觀察STFT對階次分量的分隔能力與窗參數(shù)間的關(guān)系，圖1中除了給出了σt、L、ΔM等分析參數(shù)的值外，還給出了6σf、6σt這兩個值。

圖1(a)中，輔助點“0”處的時刻為 6.82 s，在此處，相鄰兩階次分量間的頻率間隔均為6.82 Hz，等于6σf。在大于6.82 s的時段，所有相鄰兩分量之間沒有明顯的干擾現(xiàn)象，而在此時段不存在同一理論瞬時頻率下的兩相鄰階次分量間的時間間隔小于6σt的情況。

當(dāng) σt等于 200 ms，而 6σf等于 2.387 Hz時(圖 1(b)):在大于2.387 s的時段，相鄰階次分量在同一時刻的理論頻率間隔均大于6σf，但與圖1(a)不同的是，在此時段卻有干擾現(xiàn)象存在。4階分量與5階分量間的干擾出現(xiàn)在大約24 Hz以下;在24 Hz處，4階分量對應(yīng)的時刻是6 s，5 階是4.8 s，兩者間的間隔是1.2 s，等于6σt。同樣，3階分量與4階分量間的干擾出現(xiàn)在大約14.4 Hz以下;而在該頻率處3階分量對應(yīng)的時刻是 4.8 s，4 階是3.6 s，兩者間的間隔仍然是1.2 s，等于6σt。對圖1(c)也有相同的觀察，即虛假成分會出現(xiàn)在同一瞬時頻率下時間間隔小于6σt的兩分量之間。

綜上所述，窗寬為σt的高斯窗STFT幅值譜如果可以將階次分量區(qū)分開(各分量在時頻譜上占據(jù)著互不重疊的區(qū)間)，且用fspacing，min(Hz)表示信號中所關(guān)心的時頻范圍內(nèi)同一時刻分量間的最小理論頻率間隔，而用tspacing，min(s)表示同一理論瞬時頻率下的最小時間間隔，則σt應(yīng)滿足式(6):

式(6)等價于式(7):

這一現(xiàn)象可以從理論上進行定性的解釋:高斯窗函數(shù)自身的 FT的主瓣在3σf處的值是其峰值的0.105，衰減 -19.5 dB，因此當(dāng)同一時刻兩分量的理論頻率間隔大于6σf時，信號的理論頻譜與窗頻譜卷積后，兩個分量所分布的主瓣之間交疊部分很少甚至沒有，因而使它們得以區(qū)分。同樣地，在時域中，當(dāng)同一理論瞬時頻率下兩分量的時間間隔大于6σt時，這兩個分量要么不能同時處于STFT過程中的移動窗內(nèi)，即使都在一個窗內(nèi)，也是分布在時域窗的趨于零值的兩端，出現(xiàn)在STFT頻譜上兩者間干擾成分的量值就會很微弱。

在相鄰階次分量的階次差相同的情況下，轉(zhuǎn)速越低的時刻，越易違背式(7)中的最小頻率間隔限值條件;轉(zhuǎn)速越低、階次越高的區(qū)域，相鄰階次分量同一理論瞬時頻率間的時間間隔越小，越易違背式(7)中的最小時間間隔限值條件。

圖2 線性變速段信號中階次分量的理論時頻位置示意圖Fig.2 Schematic diagram for the theoretical time-frequency locations of the order components in a vibroacoustic signal with linearly varying speed

由圖2可確定

1.3.2 高斯窗STFT幅值譜幅值誤差

下面對于STFT幅值譜上能區(qū)分開階次分量的地方，再來考察其階次幅值的精度。從圖1(a)、圖1(b)、圖1(c)中的階次切片來看，這些值不宜做為階次分量瞬時幅值的估計，甚至也不能做為瞬時幅值的相對大小的估計(見圖1(b)與圖1(c)，雖然各階次理論瞬時幅值都為1，但STFT幅值譜上的各階次理論頻率處的值卻不同)。非穩(wěn)態(tài)信號的STFT幅值譜除了與窗有關(guān)，還與信號中分量的自身特征(瞬時幅值變化率、瞬時頻率變化率)密切相關(guān)。從圖1可知，分量在STFT幅值譜中所占據(jù)的頻率寬度受窗寬與分量自身特征的共同決定;同一階次分量的幅值誤差隨寬窗的增加而增加;同一窗寬下瞬時頻率變化率越大的分量，其STFT幅值誤差就越大。

2 基于Gabor階次跟蹤的階次分量瞬時幅值估計

Gabor階次跟蹤是基于離散Gabor變換與展開這一變換對進行階次分量時域重構(gòu)的技術(shù)［6］。式(9)與式(10)分別是離散Gabor變換與展開式。

其中，M是Gabor變換的時間采樣點數(shù);N是頻率分析點數(shù);L是離散窗的點數(shù);ΔM是以點數(shù)表示的時間步長。Gabor變換中的窗γ［i］被稱為分析窗，而Gabor展開中的窗h［i］被稱為合成窗。這一變換對存在的條件是h［i］與 γ［i］互為對偶函數(shù)(它們在式(9)與式(10)中的位置可以互換)，且過采樣率ros需滿足:

2.1 Gabor階次波形重構(gòu)的前提條件

如果對于一個離散分析窗，能夠計算出其對偶窗，則用式(9)與式(10)總能實現(xiàn)采樣信號的完全重構(gòu)，無論這兩個窗的時頻集中性的相似程度如何。但若想實現(xiàn)信號中的特定分量的重構(gòu)，還需要滿足如下的要求:

(2)系數(shù)譜中與想提取出波形的分量有關(guān)的系數(shù)所占據(jù)的區(qū)域與其它分量間無重疊。前面的仿真分析已表明系數(shù)譜對分量的區(qū)分能力由信號中分量特點與分析窗共同決定;文獻［3］中指出了Gabor階次跟蹤方法適用于“信號各階階比成分在時頻面上可以和其它階比成分分開時”，但并沒有指出在什么條件下可以使它們分開，而本文的第一部分回答了這個問題。

(3)合成窗應(yīng)具有與分析窗盡可能相同的時頻集中性。由于Gabor基本函數(shù)通常并不構(gòu)成正交基，所以其對偶函數(shù)不是唯一的。只有尋求到與分析窗的時頻集中性最相近的對偶窗，當(dāng)將m，n中與某分量對應(yīng)的系數(shù)保留而其它位置處的值置為0(這一過程稱為“時頻遮罩”，有關(guān)“時頻遮罩”的具體步驟可參見文獻［3］)得到矩陣，n，并用如下公式:

進行重構(gòu)所得到的波形才是與理論分量波形相近的?；趯ε即昂瘮?shù)之間的最小均方誤差(LMSE)意義上的相似性，已發(fā)展出了似正交的 Gabor展開算法［7］。研究表明，對偶窗之間的相似性與過采樣率有關(guān)，總的來說，過采樣率越大，相似性越高。在給定的離散時間步長ΔM，頻率區(qū)間點N下的高斯窗，其對偶窗與其具有最小均方誤差時，該高斯窗的離散域時間方差為［8］:

由此，可推出對應(yīng)的最佳連續(xù)時間域上的方差為:

2.2 改進的Gabor階次分量重構(gòu)流程

圖3是常規(guī)的Gabor階次波形重構(gòu)流程，它以高斯窗作為合成窗，而不是分析窗。在此過程中，要獲得區(qū)分得開階次分量的Gabor系數(shù)譜，是通過窗參數(shù)的試湊來實現(xiàn)的，沒有結(jié)合信號的特征，試湊的參數(shù)不一定能滿足將相鄰階次分量區(qū)分開的條件。

圖3 常規(guī)的階次波形重構(gòu)流程圖Fig.3 The flowchart for the conventional order waveform reconstruction

基于本文中第1部分的仿真結(jié)論，提出一種改進的階次波形重構(gòu)算法，并對重構(gòu)出的波形用Hilbert變換求包絡(luò)作為階次分量瞬時幅值的估計，全部流程如圖4所示。

圖4 基于改進的階次波形重構(gòu)算法的階次瞬時幅值估計流程圖Fig.4 The flowchart for the order component's instantaneous amplitude estimate based on the improved order waveform reconstruction algorithm

3 仿真驗證與應(yīng)用實例

為驗證此流程對變幅值階次分量的適用性，構(gòu)造了如下的仿真信號:

對此信號，如果想提取800 r/min以上的30階以內(nèi)的階次分量的瞬時幅值，按圖4所示流程，求得fspacing，min是 285.6 ms，fspacing，min是 13.3 Hz，而可用的 σt的范圍是35.8 ms～47.6 ms。

圖5是σt=40 ms時的分析結(jié)果。圖5(a)中的系數(shù)譜中30階以下的階次間沒有出現(xiàn)干擾現(xiàn)象;采用對重構(gòu)的階次波形用Hilbert變換求包絡(luò)做為瞬時幅值估計，其各階誤差在5%以內(nèi)，而用STFT單邊峰值幅值譜做為瞬時幅值估計，其誤差會受階次值影響，階次越高，誤差越大。

圖5 仿真信號S2基于改進的Gabor重構(gòu)流程的階次瞬時幅值估計結(jié)果Fig.5 Result of the instantaneous amplitude estimates of the order components in signal S2 based on the improved Gabor reconstruction

仿真試驗還表明，按圖4所示流程、選用35.8 ms～47.6 ms之間的任何一值做為窗寬得到的階次瞬時幅值估計都有與圖5(e)中一致的誤差水平，而STFT幅值譜做為瞬時幅值估計，窗寬越大，誤差越大。

對將式(15)中的SNR改為5(14 dB)的仿真信號也按圖4所示的流程做了仿真實驗，此時得到的瞬時幅值相比于小噪聲情況下的高頻噪聲要大，對該瞬時幅值進行小波降噪后再做為瞬時幅值估計，得到的各階幅值誤差在10%以內(nèi)。

圖6 某實測信號基于改進的Gabor重構(gòu)流程的階次瞬時幅值估計結(jié)果Fig.6 Result of the instantaneous amplitude estimates of the order components in a real-world signal based on improved Gabor reconstruction flowchart

圖6是一個應(yīng)用實例。試驗發(fā)動機為某直列4缸高壓共軛柴油機，標(biāo)定工況為82 kW/(3 600 r/min)，在一電控電渦流測功機上進行全負(fù)荷加速試驗。發(fā)動機從900 r/min線性加速到3 600 r/min，用時30 s，用B＆K公司數(shù)采設(shè)備同步測量振聲信號。圖6(b)中的信號為距離機體ATS側(cè)中心1 m處的聲壓。對該信號進行階次分析的參數(shù)是:相鄰階次分量的階次差:Δp=0.5;最大分析階次:pmax=10;轉(zhuǎn)速范圍:1 000 r/min～3 400 r/min(圖中只顯示了一部分)。圖6(a)是其Gabor系數(shù)譜，從中可見，主要能量集中在偶數(shù)階階次上，這與直列4缸機的不平衡的2階往復(fù)慣性力及各缸激勵轉(zhuǎn)矩合成后的主要階次為偶數(shù)階的特點是一致的。根據(jù)轉(zhuǎn)速與時間的對應(yīng)關(guān)系，從圖6(d)轉(zhuǎn)換成圖6(e)的圖形對于定位噪聲大的轉(zhuǎn)速、噪聲源識別、聲學(xué)仿真模型的驗證及零部件改進設(shè)計都有指導(dǎo)作用。

4 結(jié)論

(1)提出了工程上可用的高斯窗時間標(biāo)準(zhǔn)差與STFT幅值譜/Gabor系數(shù)譜對階次分量分隔能力之間的定量關(guān)系式。仿真與實際應(yīng)用證明了此關(guān)系式的適用性。

(2)雖然與FT諧波分析一樣采用了幅值校正，但STFT幅值不能做為瞬時幅值估計，其誤差水平會因階次、轉(zhuǎn)速變化率、窗寬而不同。

(3)采用文中提出的改進的階次波形重構(gòu)算法進行瞬時幅值估計的優(yōu)點是:①能判定是否存在對所關(guān)注的階次分量進行理想重構(gòu)的高斯窗時間標(biāo)準(zhǔn)差(見圖4中的判斷分支);②在存在適合窗寬的情況下，在小噪聲情況下，可直接用重構(gòu)的階次波形的包絡(luò)做為瞬時幅值估計;大噪聲情況下，還需對該包絡(luò)進行降噪處理，以保證階次分量瞬時幅值估計的精度。

［1］趙曉平，張令彌，郭勤濤.旋轉(zhuǎn)機械階比跟蹤技術(shù)研究進展綜述［J］.地震工程與工程振動，2008，28(6):213-219.

［2］孔慶鵬，宋開臣，陳 鷹.發(fā)動機變速階段振動信號時頻分析階比跟蹤研究［J］.振動工程學(xué)報，2005，18(4):448-452.

［3］郭 瑜，秦樹人.基于瞬時頻率估計及時頻濾波的階比分量提?。跩］.中國機械工程，2003，14(17):506-1509.

［4］田 昊，唐力偉，田 廣，等.基于核獨立分量分析的齒輪箱故障診斷［J］.振動與沖擊，2009，28(5):163-164.

［5］ Qian S.Gabor expansion for order tracking［J］.Sound and Vibration，2003，37(6):18-22.

［6］Shao H，Jin W，Qian S.Order tracking by discrete Gabor expansion［J］.IEEE Transactions on Instrumentation and Measurement，2003，52(3):754-761.

［7］ Qian S，Chen D.Discrete Gabor transform［J］.IEEE Trans.Signal Processing，1993，41(7):2429-2439.

［8］ Qian S.Introduction to time-frequency and wavelet transforms［M］.Upper Saddle River，NJ:Prentice Hall PTR，2002:41-72.