(1.石家莊鐵道大學(xué) a.電氣與電子工程學(xué)院，b.省部共建交通工程結(jié)構(gòu)力學(xué)行為與系統(tǒng)安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，河北 石家莊 050043；2.邯鄲學(xué)院 機(jī)電學(xué)院，河北 邯鄲 056005)

瞬時(shí)頻率作為非平穩(wěn)信號(hào)分析的重要參數(shù)[1]，已經(jīng)廣泛應(yīng)用于橋梁振動(dòng)檢測(cè)、地震勘測(cè)、機(jī)械、電力系統(tǒng)、雷達(dá)、通信、醫(yī)學(xué)等各方面。變轉(zhuǎn)速工況故障軸承振動(dòng)信號(hào)作為典型的非平穩(wěn)信號(hào)[2]，包含機(jī)械設(shè)備的瞬態(tài)狀況以及設(shè)備實(shí)時(shí)運(yùn)行狀態(tài)。在一些大型機(jī)械設(shè)備的故障診斷時(shí)無(wú)法安裝測(cè)速裝置，只能從振動(dòng)信號(hào)中提取瞬時(shí)轉(zhuǎn)頻信號(hào)，而從時(shí)頻分布圖中精確地估計(jì)瞬時(shí)頻率是時(shí)頻分析的重要一步。

近幾年來(lái)，在瞬時(shí)頻率估計(jì)方面提出許多方法，其中最簡(jiǎn)單的時(shí)頻脊提取算法屬于峰值檢測(cè)法[3]，但是該法最容易受噪聲成分干擾。文獻(xiàn)[4-5]中針對(duì)連續(xù)小波變換得到的時(shí)間-尺度平面，采用具有時(shí)頻點(diǎn)幅值和相鄰點(diǎn)頻率波動(dòng)的2項(xiàng)約束建立代價(jià)函數(shù)提取時(shí)頻脊。文獻(xiàn)[6]中針對(duì)時(shí)頻面提出了與其建立代價(jià)函數(shù)相似的脊線搜索法。為了避免對(duì)初始搜索點(diǎn)位置的限制，文獻(xiàn)[3]中也利用相同的代價(jià)函數(shù)建立了雙向時(shí)頻脊搜索算法。代價(jià)函數(shù)通常只利用了相鄰時(shí)頻點(diǎn)信息，一旦時(shí)頻面局部位置出現(xiàn)強(qiáng)噪聲或干擾分量的影響，很容易造成脊線提取失效。小波脊線提取方法[7]也常用于脊線提取，對(duì)單分量或多分量信號(hào)中能量較高的分量提取效果較好，但是對(duì)于多分量信號(hào)中能量較弱的分量難以識(shí)別。

為此，本文中介紹一種新的瞬時(shí)頻率提取方法。該方法利用二維時(shí)頻分布圖中時(shí)頻脊線能量分布特點(diǎn)，引入圖像處理中Steger算子進(jìn)行時(shí)頻譜中瞬時(shí)頻率脊線提?。煌ㄟ^(guò)仿真實(shí)驗(yàn)，驗(yàn)證本方法在多分量信號(hào)瞬時(shí)頻率估計(jì)方面的精度以及適應(yīng)性、抗噪性。

1 原理介紹

1.1 時(shí)頻脊線特點(diǎn)

時(shí)頻分析是研究非平穩(wěn)信號(hào)的主要工具，其中時(shí)頻脊線是瞬時(shí)頻率提取的關(guān)鍵信息。短時(shí)傅里葉變換(STFT)作為時(shí)頻分析最經(jīng)典的方法之一，對(duì)非平穩(wěn)信號(hào)x(t)，在線性空間可測(cè)、平方可積x(t)∈L2(R)，定義式為

(1)

式中:Fx(t,ω)為信號(hào)的X(t)的短時(shí)傅里葉變換,其中t為時(shí)間，ω為角頻率；gt,ω(τ)為窗函數(shù)短時(shí)傅里葉變換，τ為時(shí)間點(diǎn)；g(τ-t)為窗函數(shù)。

在時(shí)間方向，取窗函數(shù)g(τ-t)與原始信號(hào)x(t)相乘，然后進(jìn)行傅里葉變換；窗函數(shù)在時(shí)間軸上平移，完成對(duì)整個(gè)信號(hào)的短時(shí)傅里葉變換。傅里葉變換將時(shí)域信號(hào)轉(zhuǎn)化為頻域信號(hào)，通過(guò)幅值譜、相位譜和功率譜來(lái)描述信號(hào)的頻率特征。

短時(shí)傅里葉變換將一維時(shí)域信號(hào)轉(zhuǎn)換成二維時(shí)頻空間，表征時(shí)間與頻率之間的關(guān)系。若在時(shí)頻中加入能量譜密度，即可將二維空間變?yōu)槿S空間，這樣可以直觀地描述能量譜密度隨時(shí)間的變化趨勢(shì)，以及頻率與能量譜密度之間關(guān)系。圖1為一維時(shí)域信號(hào)短時(shí)傅里葉變換頻譜圖。由圖可以看出，瞬時(shí)頻率的能量呈山脊?fàn)顟B(tài)，山脊的最頂端為能量最大點(diǎn)，即峰值點(diǎn)。能量譜密度是一種概率統(tǒng)計(jì)方法，因此時(shí)頻脊能量服從正態(tài)分布。圖2為信號(hào)時(shí)頻分布圖，它是短時(shí)傅里葉變換譜圖在時(shí)頻面上的投影。由圖可以看出，在投影過(guò)程中，山脊峰值和山脊在時(shí)頻域中變換為時(shí)頻脊線，并且時(shí)頻分布也符合能量譜密度分布。

圖1 一維時(shí)域信號(hào)短時(shí)傅里葉變換頻譜圖

圖2 短時(shí)傅里葉變換時(shí)頻分布

1.2 Steger算法原理

在時(shí)頻聚集性和分辨率較好的時(shí)頻分布面上，信號(hào)的頻率切片能量在某一時(shí)刻呈現(xiàn)山脈狀態(tài)，并且服從正態(tài)分布。時(shí)頻脊線邊緣能量分布局部最小，脊線峰頂上能量分布局部最大，并且能量最大處為瞬時(shí)頻率真實(shí)值。對(duì)時(shí)頻譜瞬時(shí)頻率提取是對(duì)時(shí)頻圖像中瞬時(shí)頻率脊線分布局部導(dǎo)數(shù)求解，并沿時(shí)間方向進(jìn)行時(shí)頻脊線的提取，脊線所在位置即為瞬時(shí)頻率。Steger算法提取時(shí)頻脊線的具體過(guò)程如下。

1.2.1 Gaussian算子平滑濾波

Gaussian算子具有低通濾波特性，是對(duì)時(shí)頻分布圖像瞬時(shí)頻率脊線邊緣進(jìn)行平滑濾波，因此可以通過(guò)與Gaussian函數(shù)卷積進(jìn)行平滑濾波。Gaussian算子為

(2)

式中σ為控制濾波平滑性的參數(shù)。

將時(shí)頻圖像進(jìn)行Gaussian算子卷積計(jì)算，可使時(shí)頻分布圖像平滑濾波，并且調(diào)整參數(shù)σ可以控制濾波的平滑性，計(jì)算公式為

F(x,y)=G(x,y)*f(x,y)，

(3)

式中：f(x,y)為時(shí)頻分布圖；F(x,y)為卷積后的時(shí)頻分析圖；*為卷積運(yùn)算符。

1.2.2 構(gòu)造Hessian矩陣求取時(shí)頻脊線方向法向量

二維時(shí)頻分析圖的Hessian矩陣H(x,y)可以表示為

(4)

式中rxx、rxy、ryy為時(shí)頻圖像f(x,y)的任意一點(diǎn)的二階偏導(dǎo)數(shù)，可以用Gaussian核函數(shù)與時(shí)頻分析圖卷積運(yùn)算得到。

(5)

式中g(shù)xx(x,y)、gxy(x,y)、gyy(x,y)為Gaussian核函數(shù)的二階偏導(dǎo)數(shù)，計(jì)算公式為

(6)

Hessian矩陣的特征值應(yīng)與時(shí)頻圖像灰度函數(shù)的二階方向?qū)?shù)的極值對(duì)應(yīng)，并且極值與Hessian矩陣的特征向量方向相同，因此時(shí)頻圖像的法向量(nx,ny)與Hessian矩陣的絕對(duì)值最大的特征值的特征向量相等，同時(shí),時(shí)頻圖像函數(shù)二階方向?qū)?shù)也與Hessian矩陣的絕對(duì)值最大的特征值相等。

1.2.3 基于亞像素提取時(shí)頻圖脊線

時(shí)頻二維圖像中的某像素(x0,y0)可以用二階泰勒展開(kāi)多項(xiàng)式來(lái)表示，即

f(x,y)=r(x0,y0)+((x-x0)(y-y0))

(7)

時(shí)頻二維圖像f(x,y)在脊線邊緣方向上一階方向?qū)?shù)值為0，二階方向?qū)?shù)取極小值點(diǎn)為脊線峰值。通常用(nx,ny)來(lái)表示n(x,y)時(shí)頻邊緣方向，因此沿時(shí)頻脊線邊緣方向也可以用(nx,ny)表示，即

f[(tnx+x0),(tny+y0)]=r(x0,y0)+

(8)

針對(duì)時(shí)頻脊線邊緣，令

可得

(9)

利用條件(tnx,tny)∈[-0.5,0.5]×[-0.5,0.5]判斷一階導(dǎo)數(shù)零點(diǎn)是否位于當(dāng)前時(shí)頻像素點(diǎn)，再加上該方向上二階導(dǎo)數(shù)的極小值是否小于閾值，通過(guò)這2個(gè)判斷，就能比較準(zhǔn)確地提取時(shí)頻瞬時(shí)脊線的峰值所在時(shí)頻分布圖的位置。最后通過(guò)瞬時(shí)頻率脊線位置與時(shí)頻分布圖位置關(guān)系，實(shí)現(xiàn)瞬時(shí)頻率精確估計(jì)。

2 算法流程

根據(jù)上述分析，本文中提出一種基于Steger算子的非平穩(wěn)信號(hào)瞬時(shí)頻率估計(jì)方法，算法流程如圖3所示。具體計(jì)算過(guò)程如下：

1)將時(shí)域非平穩(wěn)信號(hào)進(jìn)行時(shí)頻變換得到二維時(shí)頻分布圖；2)對(duì)二維時(shí)頻分布圖進(jìn)行Gaussian算子平滑濾波，增強(qiáng)時(shí)頻分布圖；3)構(gòu)造Hessian矩陣并求取時(shí)頻脊線方向法向量；4)對(duì)時(shí)頻脊線進(jìn)行亞像素提取，得到瞬時(shí)頻率脊線；5)根據(jù)時(shí)頻分布圖像矩陣位置與瞬時(shí)頻率脊線峰值位置的關(guān)系估計(jì)瞬時(shí)頻率。

圖3 基于Steger算子的非平穩(wěn)信號(hào)瞬時(shí)頻率估計(jì)算法流程圖

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

3.1 多分量仿真信號(hào)設(shè)計(jì)

為了驗(yàn)證Steger算子可以提取多分量信號(hào)瞬時(shí)頻率的適應(yīng)性，設(shè)計(jì)復(fù)雜的多分量信號(hào)s，其中包括2個(gè)相位不同的正弦調(diào)頻信號(hào)s1、s2和一個(gè)非線性時(shí)變信號(hào)s3，3個(gè)信號(hào)疊加組成多分量信號(hào)，表達(dá)為s(t)=s1(t)+s2(t)+s3(t)，然后再添加Gaussian白噪聲，時(shí)域波形如圖4所示。

(10)

圖4 多分量仿真信號(hào)時(shí)域波形

3.2 適應(yīng)性

為了驗(yàn)證Steger算子在不同時(shí)頻譜提取瞬時(shí)頻率的適應(yīng)性，信號(hào)模型選用多分量信號(hào)，并添加Gaussian白噪聲，然后對(duì)多分量信號(hào)分別進(jìn)行短時(shí)傅里葉變換、同步壓縮(SST)[8]、時(shí)頻重排(RS)變換[9]、小波變換得到二維時(shí)頻分布，如圖5所示。

(a)短時(shí)傅里葉變換(b)同步壓縮(c)時(shí)頻重排(d)小波變換圖5 經(jīng)過(guò)不同方式處理得到的多分量仿真信號(hào)時(shí)頻分布

對(duì)上述多分量的二維時(shí)頻分布圖進(jìn)行Steger算子瞬時(shí)頻率估計(jì)，結(jié)果如圖6所示。由圖可以看出：從短時(shí)傅里葉變換時(shí)頻分布圖圖6(a)提取的瞬時(shí)頻率結(jié)果來(lái)看，可以清晰地看到時(shí)頻脊線受到噪聲的影響，但是在瞬時(shí)頻率估計(jì)時(shí)對(duì)二維時(shí)頻分布圖進(jìn)行Gaussian平滑濾波，使得時(shí)頻脊線銳化，減小了噪聲影響，提高瞬時(shí)頻率估計(jì)精度；在同步壓縮時(shí)頻分布圖圖6(b)提取的瞬時(shí)頻率結(jié)果中，同步壓縮是對(duì)時(shí)頻脊線進(jìn)行擠壓，擠壓過(guò)程受到噪聲的影響，使得時(shí)頻脊線擠壓程度發(fā)生了變化，瞬時(shí)頻率曲線出現(xiàn)了變形；在時(shí)頻重排時(shí)頻分布圖圖6(c)提取的瞬時(shí)頻率結(jié)果中，受噪聲的影響時(shí)頻分布圖出現(xiàn)毛刺導(dǎo)致瞬時(shí)頻率曲線出現(xiàn)小的波浪起伏；在小波變換時(shí)頻分布圖圖6(d)提取的瞬時(shí)頻率結(jié)果中，時(shí)頻分布圖受小波核函數(shù)尺度的影響，在信號(hào)瞬時(shí)頻率變化時(shí)，時(shí)頻分辨率發(fā)生明顯變化，能量出現(xiàn)發(fā)散，使得瞬時(shí)頻率始端和末端無(wú)法估計(jì)，但總體上可以精確提取。通過(guò)上述分析可知，Steger算子在提取瞬時(shí)頻率的適應(yīng)范圍比較廣泛。

(a)短時(shí)傅里葉變換(b)同步壓縮(c)時(shí)頻重排(d)小波變換圖6 對(duì)不同方法處理得到的時(shí)頻分布圖進(jìn)行Steger算子瞬時(shí)頻率估計(jì)結(jié)果

3.3 抗噪性與估計(jì)精度

為了進(jìn)一步驗(yàn)證基于Steger算子的非平穩(wěn)信號(hào)瞬時(shí)頻率估計(jì)方法在強(qiáng)噪聲下估計(jì)瞬時(shí)頻率的抗噪性與精度，并與峰值搜索法進(jìn)行對(duì)比，選用線性調(diào)頻單分量信號(hào)作為仿真信號(hào)，仿真信號(hào)角頻率和模型表達(dá)式分別見(jiàn)式(11)、(12)，添加Gaussian白噪聲，得到的時(shí)域波形如圖7所示。

ω(t)=2.5πt+30π，

(11)

(12)

圖7 單分量仿真信號(hào)時(shí)域波形

式中：η(t)為Gaussian白噪聲；t為時(shí)間，0≤t≤20 s。 將仿真信號(hào)進(jìn)行短時(shí)傅里葉變換得到二維時(shí)頻分布，如圖8所示，分別使用本文中提出的方法與峰值搜索法[10-12]進(jìn)行瞬時(shí)頻率估計(jì)，結(jié)果如圖9所示。

圖8 單分量仿真信號(hào)時(shí)頻分布

圖9 不同算法的瞬時(shí)頻率估計(jì)結(jié)果

從圖8中可以看出，由于受到噪聲的影響，時(shí)頻脊線的能量分布發(fā)生明顯的變化。當(dāng)使用峰值搜索法時(shí)，如果同一時(shí)刻受到噪聲能量的影響，使時(shí)頻脊線的峰值能量偏離真實(shí)值，將會(huì)把時(shí)頻脊線偏離值作為瞬時(shí)頻率。采用本文中提出的瞬時(shí)頻率估計(jì)方法，首先進(jìn)行Gaussian平滑濾波，可以濾除噪聲的影響，從而減小誤差。從圖9中可以明顯看出，在受到噪聲干擾時(shí)，本文中提出的方法的瞬時(shí)頻率估計(jì)結(jié)果更接近真實(shí)值，表明該方法的抗噪性能明顯優(yōu)于峰值搜索法的。

以下定量分析2種方法在不同信噪比噪聲下的瞬時(shí)頻率估計(jì)精度。計(jì)算瞬時(shí)頻率誤差率η[8]，

(13)

通過(guò)大量實(shí)驗(yàn)的測(cè)試，將上述單分量線性調(diào)頻信號(hào)中加入不同信噪比的Gaussian白噪聲，使信噪比為-20～10 dB，并將信號(hào)進(jìn)行短時(shí)傅里葉變換得到二維時(shí)頻分布圖。然后分別采用2種方法進(jìn)行瞬時(shí)頻率估計(jì)，計(jì)算2種方法在不同信噪比下瞬時(shí)頻率估計(jì)誤差，結(jié)果如圖10所示。從圖中可以看出，本文中提出的方法的瞬時(shí)頻率誤差率小于峰值搜索方法的，且誤差率均小于5%，說(shuō)明該方法具有較高的精度，同時(shí)還具有較好的抗噪性。

圖10 不同信噪比下不同算法的瞬時(shí)頻率誤差

Rényi熵可以定量分析時(shí)頻聚集性，是評(píng)價(jià)時(shí)頻分辨率的重要指標(biāo)[13]，Rényi熵越小，時(shí)頻能量聚集性越好，同時(shí)時(shí)頻分辨率越高。對(duì)于維數(shù)為α的Rényi熵Hα(X)，當(dāng)α≥0和α≠1時(shí)，定義式為

(14)

式中X為隨機(jī)變量，對(duì)應(yīng)的概率密度為pi(X=i)，i=1,2,…,n。

對(duì)不同Rényi熵的時(shí)頻分布圖，分別采用本文中提出的算法和峰值搜索法進(jìn)行瞬時(shí)頻率估計(jì)，然后計(jì)算2種方法在不同Rényi熵時(shí)瞬時(shí)頻率估計(jì)誤差[9]，結(jié)果如圖11所示。從圖中可以看出，本文中提出的方法的瞬時(shí)頻率誤差率均明顯小于5%，遠(yuǎn)遠(yuǎn)小于峰值搜索方法的。

圖11 采用不同算法時(shí)不同Rényi熵的時(shí)頻分布的瞬時(shí)頻率誤差

4 實(shí)測(cè)信號(hào)分析

為了進(jìn)一步驗(yàn)證本文中提出的基于Steger算子的非平穩(wěn)信號(hào)瞬時(shí)頻率估計(jì)方法的實(shí)用性和有效性，將該方法應(yīng)用于實(shí)驗(yàn)采集的滾動(dòng)軸承振動(dòng)信號(hào)來(lái)進(jìn)行瞬時(shí)轉(zhuǎn)頻估計(jì)。采用QPZZ-Ⅱ型旋轉(zhuǎn)機(jī)械振動(dòng)及故障模擬實(shí)驗(yàn)平臺(tái)進(jìn)行實(shí)驗(yàn)，采用NI PXIe 4496型數(shù)據(jù)采集設(shè)備進(jìn)行數(shù)據(jù)采集，加速度振動(dòng)傳感器采樣頻率為25.6 kHz，參考軸轉(zhuǎn)速傳感器為激光轉(zhuǎn)速計(jì)，采樣頻率為1 kHz。根據(jù)激光轉(zhuǎn)速計(jì)采集的信號(hào)，采用五點(diǎn)式公式求解真實(shí)轉(zhuǎn)頻，并與估計(jì)的瞬時(shí)轉(zhuǎn)頻進(jìn)行比較。為了充分驗(yàn)證方法的有效性，對(duì)升速振動(dòng)信號(hào)進(jìn)行分析。實(shí)驗(yàn)滾動(dòng)軸承的型號(hào)為N205EM，有外圈故障，該滾動(dòng)軸承有關(guān)參數(shù)見(jiàn)表1。

表1 N205EM型滾動(dòng)軸承的基本參數(shù)

模擬啟動(dòng)工況，轉(zhuǎn)速由升速過(guò)程過(guò)渡到勻速過(guò)程，采樣時(shí)間為10 s，升速過(guò)程轉(zhuǎn)頻從11.4 Hz增大到24.6 Hz。原始信號(hào)時(shí)域波形如圖12所示。

圖12 實(shí)測(cè)振動(dòng)信號(hào)的時(shí)域波形

實(shí)驗(yàn)采集到的振動(dòng)信號(hào)包含軸承故障沖擊分量，且轉(zhuǎn)速信號(hào)為低頻信號(hào)。為了提高運(yùn)算速率并有效提取低頻轉(zhuǎn)速信息，首先對(duì)采集的振動(dòng)信號(hào)進(jìn)行低通濾波處理，低通濾波截止頻率為1 300 Hz，然后進(jìn)行降采樣處理，其倍數(shù)為5。對(duì)預(yù)處理的振動(dòng)信號(hào)進(jìn)行小波變換得到二維時(shí)頻分布，如圖13所示。

將得到小波時(shí)頻分布圖進(jìn)行Steger算子提取瞬時(shí)轉(zhuǎn)頻，由于受到倍頻的影響，因此從圖13中可以明顯看出倍頻分量；但是，由于倍頻的能量低于實(shí)際轉(zhuǎn)頻能量，因此在使用Steger算子時(shí)，倍頻的虛假分量可以被Gaussian平滑濾波濾除，從而減少倍頻分量的影響。最終得到瞬時(shí)轉(zhuǎn)頻估計(jì)結(jié)果如圖14所示。從圖中可以直觀地看出，采用本文中提出的方法得到的轉(zhuǎn)頻曲線完全與實(shí)際轉(zhuǎn)頻曲線重合，計(jì)算得到的估計(jì)誤差率為1.32%，表明該方法在瞬時(shí)頻率估計(jì)方面具有較高的精度，并且對(duì)倍頻的抗噪性效果非常好。

圖13 實(shí)測(cè)振動(dòng)信號(hào)的小波時(shí)頻分布

圖14 實(shí)測(cè)振動(dòng)信號(hào)的瞬時(shí)轉(zhuǎn)頻估計(jì)結(jié)果

5 結(jié)論

1)本文中理論和實(shí)驗(yàn)驗(yàn)證了Steger算子能夠準(zhǔn)確進(jìn)行瞬時(shí)頻率估計(jì)，提供了一種新穎的瞬時(shí)頻率估計(jì)方法，并且驗(yàn)證了該方法在瞬時(shí)頻率估計(jì)時(shí)對(duì)時(shí)頻分布的適應(yīng)性、抗噪性、精確性良好。

2)本文中提出的方法對(duì)時(shí)頻分布圖適應(yīng)性強(qiáng)，抗噪性優(yōu)，具有獨(dú)特的優(yōu)勢(shì)。該方法可用于無(wú)法安裝轉(zhuǎn)速傳感器的機(jī)械傳動(dòng)裝置，通過(guò)采集振動(dòng)信息進(jìn)行瞬時(shí)轉(zhuǎn)頻估計(jì)，將瞬時(shí)轉(zhuǎn)頻作為后期進(jìn)行多階比分析重要環(huán)節(jié)，對(duì)滾動(dòng)軸承故障診斷具有重要意義。

3)在惡劣工況下，機(jī)械振動(dòng)信息中包含強(qiáng)背景噪聲和很多的干擾成分，在利用本文中提出的方法進(jìn)行瞬時(shí)頻率估計(jì)時(shí)，需要先得到高分辨率和清晰的時(shí)頻分布圖，這樣才能獲得高精度的瞬時(shí)頻率。