丁嘉凱, 王 義, 穆志國, 張光耀, 李 城

(1.重慶大學 機械與運載工程學院,重慶 400030; 2.重慶大學機械傳動國家重點實驗室,重慶 400030;3.中航西安飛機工業(yè)集團股份有限公司,西安 710089)

旋轉機械被廣泛應用在航空發(fā)動機、風力發(fā)電機和直升機等傳動系統(tǒng)。由于在實際運行過程中,其工況經(jīng)常發(fā)生變化,常常處于變載荷、變轉速的工況下,工作環(huán)境及其惡劣,極容易發(fā)生故障[1-5]。目前,越來越多的旋轉機械處在變轉速工況下進行工作,如汽車的升降速、飛機的起落過程和航空發(fā)動機運行過程等等。當旋轉機械處在變轉速工況下,所測得的振動信號將失去周期性變化的規(guī)律,并呈現(xiàn)出非平穩(wěn)性出現(xiàn)調頻、調幅與調相等現(xiàn)象[6]。傳統(tǒng)的頻譜分析方法僅能適用于平穩(wěn)信號的分析,所以針對變轉速下的旋轉機械故障診斷基于恒定轉速的分析方法無法達到目的。因此,如何從非平穩(wěn)振動信號中提取旋轉機械是進行故障診斷的一個重要問題[7-9]。

近年來,國內外許多學者對非平穩(wěn)振動信號特征提取方法進行了廣泛的研究,其主要方法分別為階次跟蹤方法與時頻分析方法[10-13]。其中階次跟蹤方法主要是對原始信號進行等角度重采樣,轉化為平穩(wěn)信號后再進行頻譜分析得到信號階次譜。但是傳統(tǒng)階次分析中轉速需要通過轉速計測取,由于測取轉速受到安裝成本以及安裝空間限制,故許多學者在原始信號中進行轉速提取,構成無鍵相階次分析方法[14-16]。雖然無鍵相階次分析可以將時域非平穩(wěn)信號轉化為角域平穩(wěn)信號,恢復原始信號的周期性,但是等角度重采樣過程較為復雜,計算效率較低[17],并且從階次分析結果可得其存在包絡畸變現(xiàn)象,故階次分析存在一定的誤差[18]。

時頻分析方法對于提取非平穩(wěn)信號的時變特性有很大的優(yōu)勢,它是一種能在時頻域中定量又能較為直觀地描述信號故障特征的分析方法,所以針對變工況下非平穩(wěn)振動信號時頻分析是一種有效的分析手段[19]。但是傳統(tǒng)時頻分析方法各有優(yōu)缺點,如短時傅里葉變換(short-time Fourier transform, STFT)和小波變換(wavelet transform, WT)本質上是一種以時頻平面上的靜態(tài)分辨率為特征的線性變換,它被細分為一個恒定面積的基本單元,雖然不存在交叉項干擾的問題,但是其時頻分辨率受到Heisenberg不確定性原理的限制,不能在時頻分辨率上同時達到最佳[20]。Wigner-Ville分布(Wigner-Ville distribution, WVD)雖然在時頻分辨率上較STFT和WT高,但是WVD存在固有的交叉項干擾問題,不適用于分析復雜多分量非平穩(wěn)信號[21]。以WVD為基礎擴展的雙線性時頻分布如Cohen類分布和仿射類分布,它通過各類核函數(shù)進行平滑處理,雖然抑制了交叉項干擾,但是其會造成原始信號發(fā)生畸變,從而造成時頻分辨率降低[22]。Hilbert-Huang變換對信號的時變具有自適應性,并且具有良好的局部時頻聚集能力,但是針對信號中的奇異點較為敏感,容易產生模態(tài)混淆等影響,造成產生虛假的本征模態(tài)函數(shù)(intrinsic mode function, IMF),最終影響分析結果[23]。為了處理以上問題,調頻變換(chirplet transform, CT)應運而生,它是針對用于分析線性瞬時頻率(instantaneous frequency, IF)恒定的調頻信號[24]。但是,當待分析信號具有非線性IF軌跡時,CT也不能保證時頻表征的精度。為了解決這個問題,Peng等[25]提出了一種多項式chirplet變換(polynomial chirplet transform, PCT)時頻分析方法,它通過使用多項式函數(shù)代替線性調頻核,可以逼近非線性IF軌跡,最終得到時頻能量集中較好的時頻表征。雖然PCT在分析強非線性IF信號時表現(xiàn)出了較強的時頻分辨率,但是在變工況下,實際振動信號的復雜性和多樣性導致我們不能使用一個單獨的多項式模型來描述所有的信號分量。此外,由于PCT需要在傳統(tǒng)的時頻表征方法上進行IF估計,因此需要初始化的時頻表征具有一定的能量集中度,否則無法達到精細化時頻表征[26]。國內外針對IF脊線精細化表征做了許多研究[27],但是在強噪聲和多分量的情況下,由于噪聲和無關分量的影響對于弱分量提取較為困難。

本文針對CT在分析非平穩(wěn)振動信號時其核函數(shù)固定不變導致IF軌跡無法達到精細化表征的問題,本文提出了一種MSBCT精細化時頻表征方法,理論分析結果表明該方法可利用高階相位算子將核函數(shù)精確逼近IF軌跡,使得該方法在變轉速工況下針對非線性IF軌跡可得到精細化時頻表征。并通過仿真和試驗驗證該方法可適用于處理強噪聲、多分量的信號,并且可應用于變轉速旋轉機械故障診斷,能夠直觀地提取故障特征頻率。

1 多縮放基調頻變換

1.1 CT理論

信號x(t)∈L2(R)的chirplet變換基礎理論公式可表示為

(1)

式中:s(t)為信號x(t)經(jīng)Hilbert變換得到的解析信號,s(t)=x(t)+jH(x(t));φ(t)為相位函數(shù);h(t)為一個歸一化的時不變高斯窗函數(shù)并且h(t)∈L2(R)。采用的高斯窗函數(shù)可表示為

(2)

式中,σ代表高斯窗函數(shù)的標準差。

根據(jù)Ville[28]提出的理論可知,信號的瞬時頻率可以由解析信號的相位導數(shù)來計算得到,則解析信號s(t)可以表示為

(3)

(4)

式中:t0∈R為選擇的時間中心;φ(t0)∈R為t0時刻的信號頻率大小,令C=φ′(t0),在CT中C∈R為調頻率。CT核心為通過旋轉頻率來匹配信號中的瞬時頻率(IF)軌跡,從而獲得具有較高能量濃度的時頻分布(time-frequency representation, TFR)。將式(4)進行二次求導可得調頻率C的具體表達式如下所示

φ(t0)+C(t-t0)

(5)

(6)

式中:α∈[-π/2,π/2]為旋轉角度;C為一個常數(shù)。

當旋轉角度的正切值tanα等于信號的IF軌跡的斜率dφ(t)/dt時,對應的TFR的能量濃度最高。CT調頻率與信號的IF軌跡斜率匹配原理如圖1所示。圖1時頻圖中包含信號分量φ1(t)。φ1(t)IF軌跡斜率為-tanα,CT采用的是一個時不變高斯窗進行時頻變換,其中窗函數(shù)的窗長固定不變?yōu)閃L。當C為0時,此時的CT即為STFT,如圖1中虛線所示,此時的CT窗函數(shù)與信號分量存在一個夾角α,得到的TFR會存在一定的能量泄露造成頻譜模糊效應,得不到具有較高能量濃度的TFR。當CT中的C為-tanα時,此時C等于信號IF軌跡的斜率,造成的能量泄露較少可得到信號較高能量濃度的TFR,如圖1中實線部分。

圖1 傳統(tǒng)CT調頻率與信號的IF軌跡斜率匹配原理Fig.1 Matching principle of traditional CT kernel function and signal IF trajectory

1.2 CT方法的局限性

由于CT變換的調頻率C是固定值,所以它在處理具有相同IF軌跡斜率的多分量信號具有優(yōu)勢,它能得到具有較高的能量濃度的TFR。但是實際信號中通常不會存在一致的IF軌跡斜率多分量信號,比如機械裝備中軸承在變轉速的條件下運行時,出現(xiàn)故障后會產生不同的故障頻率,這些故障IF軌跡斜率通常不一致。所以,CT在處理不同IF軌跡斜率多分量信號也會產生能量泄露導致頻譜模糊效應。

因此,在具有不同IF軌跡斜率的多分量信號中進行CT變換時,要在同一時刻不同頻率分量上達到最佳的能量濃度是不可能實現(xiàn)的,這是由于CT的理論基礎決定的。故需要對傳統(tǒng)CT變換進行改進,需要針對具有不同IF軌跡斜率的多分量信號進行C的自適應選取,以至于得到最佳能量濃度的TFR,為后續(xù)的變工況下機械裝備故障診斷與非平穩(wěn)信號時頻脊線精確提取奠定基礎。

1.3 MSBCT方法

針對不同IF軌跡斜率的多分量信號,其IF函數(shù)利用泰勒公式進行展開,文獻[26]與文獻[29]均將其展開至低階形式,本文利用泰勒公式將IF函數(shù)展開至高階形式,可針對多個信號分量IF軌跡斜率進行匹配,并根據(jù)實際信號在分析頻率范圍內存在多個分量信號的特性,得到MSBCT方法中高階相位算子具體形式如下所示

φ(t0)+φ′(t0)(t-t0)+

(7)

聯(lián)合式(3)和式(7)可得MSBCT方法的具體表達式為

exp(-j2πφ(t))dt

(8)

由于信號的IF不能提前知道,但是IF軌跡斜率的角度可以知道其范圍,其角度范圍從-π/2變化到π/2。則多信號IF軌跡斜率角度可設置為

(9)

(10)

(11)

C1=Atanα

(12)

C2=Btanαtanβ

(13)

C3=Ctanαtanβtanθ

(14)

式中,A,B,C分別為等式系數(shù),即存在系數(shù)A,B,C使得式(12)～(14)成立。

將式(12)～(14)代入式(8)中,可得:

Atanα(t-t0)2+…+Btanαtanβ(t-t0)3+…+

Ctanαtanβtanθ(t-t0)4))·…·h(t-t0)·

exp(-j2πφ(t))dt

(15)

由于解析信號擴展為多項式形式,則MSBCT方法中的核函數(shù)需進行泰勒展開并且與解析信號中多項式抵消項,則可得到TFR的最佳峰度,其最佳峰度表達式為

(16)

式中,V為TFR的最大范圍。

此時得到的即為最佳濃度的TFR,所以最終MSBCT表達式為

exp(-j2π(t+Atanα(t-t0)2+…+

Btanαtanβ(t-t0)3+…+

Ctanαtanβtanθ(t-t0)4))dt

(17)

本文確定MSBCT參數(shù)的有效方法為通過生成一系列IF軌跡斜率夾角α,β和θ的組合來生成TFR。對于每個時間點找出TFR能量濃度最高的參數(shù)組合作為最佳參數(shù)組合。其中通過式(16)來尋找TFR最優(yōu)能量濃度,然后依次迭代確定每個時間點的最優(yōu)參數(shù)組合,從而得到整段信號的最優(yōu)MSBCT參數(shù)組合。其中,MSBCT最優(yōu)參數(shù)組合表達式為

(18)

2 性能評估

該部分本文利用仿真信號對MSBCT算法進行性能評估,驗證其是否具有在時變工況下對時頻脊線精細化表征的能力,并且利用最新的時頻變換方法與MSBCT算法進行性能比較,如VSLCT,SBCT,CWT和STFT。

本文利用一個IF軌跡斜率時變的多分量信號進行MSBCT時頻表征,從而驗證MSBCT算法可針對時變工況下IF軌跡進行精確化時頻表征。本章采用的仿真信號s1(t)如下所示

(19)

式中,φ1(u)為基頻。其表達式如下

(20)

在此案例中仿真信號s1(t)的采樣頻率為100 Hz,采樣時間為50 s。并且其他倍頻與基頻之間的關系如下

(21)

仿真信號s1(t)的時域圖和頻域圖如圖2所示。s1(t)為調頻調幅信號,所以利用MSBCT算法對仿真信號s1(t)進行精確時頻表征,其TFR如圖3所示。MSBCT算法可將仿真信號s1(t)進行精確時頻表征,在A部分IF軌跡緊鄰分量,MSBCT將4條緊鄰分量進行精確的時頻表征,沒有產生任何時頻模糊現(xiàn)象。在B部分IF軌跡斜率突變處也得到了清晰的時頻表征圖。所以圖3的結果可以看出,MSBCT可在時變工況下將時頻脊線進行精細化表征。

(a)

(b)圖2 仿真信號s1(t)的時域圖和頻域圖Fig.2 Time-domain and frequency-domain diagrams of the simulation signal s1(t)

圖3 利用MSBCT算法對仿真信號s1(t)的TFRFig.3 The TFR of s1(t) by using MSBCT algorithm

利用不同的時頻變換方法對仿真信號s1(t)進行時頻表征,其TFR如圖4所示。其中,SBCT與STFT的窗口長度設置為128。如圖4(a)所示,從整體看來,VSLCT算法得到的TFR均存在時頻模糊問題,尤其是在A和B兩個部分,其中在A處發(fā)生了頻譜模糊效應。在B處發(fā)生了頻散現(xiàn)象,使得此處的能量發(fā)生了泄露。如圖4(b)所示,SBCT在A部分IF軌跡相隔較近處發(fā)生了時頻模糊現(xiàn)象,在其余部分獲得了較好的TFR,但是其時頻能量集中度較差,其整體效果不如MSBCT算法。如圖4(c)所示,從整體可得CWT算法得到的時頻能量集中較VSLCT和MSBCT效果較差,并且在A部分發(fā)生了頻譜模糊現(xiàn)象,這導致了CWT無法對IF軌跡緊鄰分量信號進行準確的時頻表征。如圖4(d)所示,STFT所得到的時頻圖和CWT的類似,并且STFT所得到的時頻能量集中較CWT更差,其中第一條和第二條IF軌跡在整個時間段均發(fā)生了頻譜模糊現(xiàn)象,故STFT也無法對IF軌跡緊鄰分量信號進行精確的時頻表征。表1為利用MSBCT和不同時頻變換方法對仿真信號s1(t)分析所得到的運行時間。由表1可知,本文所提MSBCT算法相比于CWT和STFT傳統(tǒng)時頻變換方法所用時間較長,但是CWT和STFT所得到的時頻表征效果不佳,無法得到精細化時頻表征,并且時頻能力不集中。但是,MSBCT所運行時間相比于SBCT和VSLCT等精細化時頻表征方法較短,可突出MSBCT算法在精細化時頻表征方法中的優(yōu)越性。

(a) VSLCT

(b) SBCT

(c) CWT

(d) STFT圖4 仿真信號s1(t)的TFRFig.4 TFR of simulation signal s1(t)

表1 不同時頻變換方法分析仿真信號s1(t)所用時間Tab.1 The running time of simulation signal s1(t) obtained by different time-frequency transform method

為了進一步得到MSBCT算法對噪聲抑制的效果,本文在式(19)的基礎上加入SNR=-5 dB的高斯白噪聲進行時頻表征,其公式如下所示

(22)

式中,n(t)為信噪比為-5 dB的高斯白噪聲。

對于加入SNR=-5 dB的高斯白噪聲仿真信號s1(t),本文利用MSBCT和VSLCT對其進行時頻表征,所得到的時頻圖如圖5所示。如圖5(a)所示,對于加入-5 dB噪聲再利用MSBCT對其進行時頻表征后,從整體所看四條IF軌跡得到了精確表征,并且時頻能量集中較好,并且在A和B兩個部分的IF軌跡時頻表征均較為準確,充分說明MSBCT對噪聲抑制有著較強的優(yōu)勢。如圖5(b)所示,VSLCT算法在A和B兩個部分均發(fā)生了時頻模糊現(xiàn)象,并在整個時頻面上對于噪聲抑制沒有MSBCT算法效果好。

(a) MSBCT

3 試驗驗證

在本章內容中MSBCT算法被應用至變轉速軸承內圈故障信號精細化時頻表征,本文采用渥太華大學SpectraQuest機械故障模擬器去生成變轉速軸承內圈故障信號[30],其試驗臺如圖6所示。該試驗臺由電機、編碼器、測試軸承、健康軸承、AC和加速度傳感器組成。其中測試軸承型號為ER16K,軸承節(jié)圓直徑為38.52 mm,滾子直徑為7.94 mm。加速度傳感器安裝在測試軸承上方,便于采集變轉速軸承故障信號。其中,本試驗采樣頻率為200 kHz,采樣時間為10 s。在此案例中,本文選取軸承轉速先升后降的工況,具體為軸承轉速從14.8 Hz升高至21.7 Hz,然后降至13.6 Hz。由文獻[31]提供該軸承內圈故障階次為5.43,故該軸承內圈一階故障頻率fI具體變化為從80.36 Hz升高至117.83 Hz,然后降至73.85 Hz。本文選取的變轉速滾動軸承內圈故障信號時域圖和頻域圖如圖7所示。由于滾動軸承的轉速是隨著時間變化的,導致其發(fā)生頻譜模糊效應,所以利用MSBCT算法對其進行時頻表征,在時頻域上對其進行滾動軸承內圈故障IF軌跡進行分析,從而進行故障診斷。

圖6 變轉速滾動軸承內圈故障試驗臺Fig.6 Rolling bearing inner ring fault test rig under variable speed

(a)

(b)圖7 變轉速滾動軸承內圈故障信號的時域圖和頻域圖Fig.7 Time-domain and frequency-domain diagrams of variable speed rolling bearing inner ring fault signal

在進行時頻表征之前,本試驗進行以下操作,由于采樣頻率過高,故對原始信號進行帶通濾波和降采樣操作。首先本文利用‘db10’小波基的5層小波包算法對原始振動信號進行帶通濾波操作,選取第一層濾波信號進行降6倍采樣操作得到最終分析信號,再利用MSBCT對信號進行時頻表征如圖8所示。滾動軸承內圈故障fI,2fI,3fI和4fIIF軌跡均被MSBCT算法精確表征,并且在整個時頻面上MSBCT算法對噪聲有一定程度上的抑制作用,對于該方法得到的2fI來說,MSBCT有著較好的時頻能量集中。對于圖8中的A和B兩部分,該處的IF軌跡斜率發(fā)生了突變,MSBCT算法也可將其精確表征。再利用不同的時頻變換算法與MSBCT進行效果對比,其時頻圖如圖9所示。

圖8 利用MSBCT算法對變轉速滾動軸承內圈故障信號的TFRFig.8 TFR of rolling bearing inner ring fault signal under variable speed by using MSBCT algorithm

(a) VSLCT

(b) SBCT

(d) STFT圖9 變轉速滾動軸承內圈故障信號的TFRFig.9 TFR of rolling bearing inner ring fault signal under variable speed

在本案例中SBCT和STFT的窗長均設置為128。如圖9(a)所示,滾動軸承內圈故障fI,2fI和3fI可以從時頻圖中得到,但是4fI較為模糊無法準確得到。從整體看來,VSLCT算法得到的TFR均存在時頻模糊現(xiàn)象,尤其是在初始時間端和結尾時間段均發(fā)生了較為嚴重的時頻模糊現(xiàn)象。其中在A和B處時頻表征發(fā)生了時頻模糊問題,造成了頻譜模糊效應。如圖9(b)所示,SBCT算法可將滾動軸承內圈故障fI,2fI和3fI從時頻圖中大致得到,但是4fI較為模糊無法準確得到。并且SBCT得到的TFR中IF軌跡時頻能量濃度較差,時頻面模糊問題較為嚴重,噪聲抑制效果較差。如圖9(c)所示,滾動軸承內圈故障fI和2fI可以從時頻圖中大致得到,但是其IF軌跡也發(fā)生了模糊,3fI和4fI較為模糊無法準確得到。從整體可得CWT算法得到的時頻能量集中較VSLCT和MSBCT效果較差,整體受噪聲影響較為嚴重。并且在A和B處CWT無法對其進行準確的時頻表征。如圖9(d)所示,滾動軸承內圈故障fI,2fI和3fI可以從時頻圖中大致得到,但是4fI較為模糊無法準確得到。并且STFT所得到的時頻圖較CWT的時頻能量較為集中,但是STFT所得到的時頻圖效果受到噪聲影響較為嚴重,在整個時頻面上均發(fā)生了時頻模糊現(xiàn)象,時頻能量集中較MSBCT更差。表2為利用MSBCT和不同時頻變換方法對變轉速滾動軸承內圈故障信號分析所得到的運行時間。由表2可知,MSBCT所運行時間較VSLCT和SBCT等精細化時頻表征方法較短。

表2 不同時頻變換方法分析變轉速滾動軸承內圈故障信號所用時間Tab.2 The running time of variable speed rolling bearing inner ring fault signal obtained by different time-frequency transform method

綜上所述,針對變轉速軸承內圈故障信號利用MSBCT算法可對其進行精細化時頻表征,并且在IF軌跡斜率突變處的時頻能量較為集中,時頻表征整體沒有發(fā)生模糊現(xiàn)象,與其他算法進行對比突出了它在變工況下對時頻脊線精細化表征的優(yōu)勢。

4 結 論

本文提出了一種新的時頻分析方法,稱為MSBCT。它克服了傳統(tǒng)CT變換的缺點,主要針對在變工況下對非線性IF軌跡分量有著較高的時頻能量集中度,并且在時頻域內可一定程度上抑制噪聲信號的影響,克制時頻模糊現(xiàn)象。在MSBCT算法中,理論分析結果表明該方法可利用高階相位算子將MSBCT方法核函數(shù)精確逼近IF軌跡,使得MSBCT算法可針對變轉速工況下的IF軌跡得到精確的時頻表征。仿真結果和試驗結果均得到非線性IF軌跡的精細化時頻表征,并且與最新的時頻分析方法進行比較驗證了MSBCT算法的優(yōu)勢性。