金 陽(yáng)，郝志勇，鄭 旭，景國(guó)璽

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

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

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

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

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

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

1.1 離散STFT幅值譜

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

1.2 高斯窗

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

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

其中:

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

1.3 仿真分析

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

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

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

圖1 信號(hào)S1在不同窗寬下的STFT單邊峰值幅值譜(L:窗的點(diǎn)數(shù);△M:以點(diǎn)數(shù)表示STFT的時(shí)間步長(zhǎng))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幅值譜的對(duì)階次分量的分隔能力

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

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

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

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

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

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

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

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

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

圖2 線性變速段信號(hào)中階次分量的理論時(shí)頻位置示意圖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幅值譜幅值誤差

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

2 基于Gabor階次跟蹤的階次分量瞬時(shí)幅值估計(jì)

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

其中，M是Gabor變換的時(shí)間采樣點(diǎn)數(shù);N是頻率分析點(diǎn)數(shù);L是離散窗的點(diǎn)數(shù);ΔM是以點(diǎn)數(shù)表示的時(shí)間步長(zhǎng)。Gabor變換中的窗γ［i］被稱(chēng)為分析窗，而Gabor展開(kāi)中的窗h［i］被稱(chēng)為合成窗。這一變換對(duì)存在的條件是h［i］與 γ［i］互為對(duì)偶函數(shù)(它們?cè)谑?9)與式(10)中的位置可以互換)，且過(guò)采樣率ros需滿(mǎn)足:

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

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

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

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

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

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

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

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

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

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

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

3 仿真驗(yàn)證與應(yīng)用實(shí)例

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

［2］孔慶鵬，宋開(kāi)臣，陳 鷹.發(fā)動(dòng)機(jī)變速階段振動(dòng)信號(hào)時(shí)頻分析階比跟蹤研究［J］.振動(dòng)工程學(xué)報(bào)，2005，18(4):448-452.

［3］郭 瑜，秦樹(shù)人.基于瞬時(shí)頻率估計(jì)及時(shí)頻濾波的階比分量提取［J］.中國(guó)機(jī)械工程，2003，14(17):506-1509.

［4］田 昊，唐力偉，田 廣，等.基于核獨(dú)立分量分析的齒輪箱故障診斷［J］.振動(dòng)與沖擊，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.