突變信號的時幅拐點分析法*

戴屹梅 張和生 李東 齊紅梅 方軻

(1.北京交通大學(xué) 電氣工程學(xué)院， 北京 100044； 2.華南理工大學(xué) 軟件學(xué)院, 廣東 廣州 510006)

突變信號的時幅拐點分析法*

戴屹梅1張和生1李東2齊紅梅1方軻1

(1.北京交通大學(xué) 電氣工程學(xué)院， 北京 100044； 2.華南理工大學(xué) 軟件學(xué)院, 廣東 廣州 510006)

為解決發(fā)動機高速軸承振動數(shù)據(jù)時域信號中隱藏的突變信號難以精確時頻定位的問題，提出一種突變信號的時幅拐點分析算法.該算法首先通過離散傅里葉變換對信號進行特征頻率的提取，再針對該特征頻率進行時移傅里葉分析，通過該頻率時移傅里葉分析下的時幅拐點找出突變信號的進入時刻和消失時刻，最后通過對確定時段內(nèi)的信號進行二次傅里葉分析獲得該突發(fā)信號的幅值.理論分析、理論仿真和高速軸承振動數(shù)據(jù)工程仿真驗證均表明，該算法能準確捕捉和提取突變信號的幅值及其在振動信號中的出現(xiàn)時刻和消失時刻.

突變信號；傅里葉變換；時幅拐點；頻率分析

目前，在信號分析中人們普遍認為存在一種測不準關(guān)系，即：信號不能同時在時域和頻域準確定位[1].傳統(tǒng)的博里葉變換是一種基于平穩(wěn)信號的全局性變換，無法反映信號的局部特征，即無法反映非平穩(wěn)信號頻率隨時間變化的特征.

自20世紀80年代以來，隨著現(xiàn)代信號處理理論和方法的逐步發(fā)展，對于非平穩(wěn)信號的分析方法逐漸豐富起來，時頻域分析法成為目前最為有效的非平穩(wěn)信號處理方法.其中短時傅里葉變換(STFT)、小波變換、S變換和廣義S變換等方法目前已廣泛應(yīng)用于非平穩(wěn)信號的分析[2- 7].

短時傅里葉分析是為解決傳統(tǒng)傅里葉分析無法獲取信號的局部特征而提出來的.它通過引入一個時間窗函數(shù)來獲取信號頻率隨時間的變化情況.文獻[8]提出一種基于短時傅里葉變換的能量譜和獨立分量分析的抗干擾滾動軸承包絡(luò)分析方法，該方法首先對獲取的信號進行STFT能量譜分析，然后進行帶通濾波并提取濾波后的包絡(luò)波形，再通過ICA實現(xiàn)包絡(luò)波形分離去噪，最后通過匹配各獨立分量的包絡(luò)頻譜和理論計算的故障特征頻率來實現(xiàn)特征提取.

小波變換是目前使用最為廣泛的時頻分析方法，關(guān)于該方法的研究成果也頗為豐富，文獻[9]中采用小波時頻分析的方法，對航空發(fā)動機實測數(shù)據(jù)進行了分析，發(fā)現(xiàn)其在不同頻率范圍內(nèi)含有不同的沖擊成分；然后采用諧波小波包對不同的頻帶進行分解，成功地得到了其故障特征.文獻[10]利用小波分解重構(gòu)具有故障特征的信號，再通過Hilbert變換進行解調(diào)和細化頻譜分析，能有效地提取特征頻率.文獻[11]中利用小波包分解對軸承的動態(tài)信號進行分析、提取特征，采用RBF神經(jīng)網(wǎng)絡(luò)進行軸承故障診斷.文獻[12]對小波和Hilbert變換、小波包變換、EMD經(jīng)驗?zāi)B(tài)分解方法進行了深入研究，提出小波包和EMD分解相結(jié)合的故障診斷方法.該方法首先利用小波包分解去噪，重組后進行EMD分解，分析前幾個IMF，提取故障脈沖的有用信息.然后找到與原始信息最相關(guān)的IMF，進行小波包分解，最后選擇合適的小波包進行快速傅里葉變換，提取特征頻率.

S變換和廣義S變換是近年來新興起的信號特征提取方法，該方法是小波變換和短時傅里葉變換的繼承和發(fā)展.文獻[13]基于STFT和連續(xù)小波變換的理論基礎(chǔ)推導(dǎo)出連續(xù)S變換，并利用快速傅里葉變換(FFT)實現(xiàn)S變換離散化計算，仿真結(jié)果表明該方法能夠提取低信噪比信號中的沖擊特征.文獻[14]對采集到的軸承信號進行S變換，并對變換結(jié)果s進行奇異值分解，成功提取了軸承信號特征.

為了彌補傳統(tǒng)頻譜分析在局部時域上的局限性，本研究基于快速傅里葉變換，提出了一種突變信號的時幅拐點分析方法.該方法先通過降低閾值來獲取全時域下傅里葉變換的特征頻率，然后針對某一個頻率進行時移傅里葉分析，即在全時域(T1，T2)內(nèi)，從T1時刻開始，以某一時間步長為單位逐步后移逼近T2，獲取每一個時段內(nèi)分析得到的頻率幅值，最后將所有獲得的頻率幅值與分析時段形成一條時幅曲線.若該頻率信號為突發(fā)信號，則可根據(jù)其時幅曲線的變化特征及其拐點來判斷該信號的存在時間，最后對該頻率在出現(xiàn)時段內(nèi)進行FFT分析，獲得該頻率對應(yīng)的信號幅值，從而較為準確地獲得突發(fā)信號的特征參數(shù).

1 時幅拐點分析算法

1.1 信號疊加的基本原理

在討論突變信號的頻率幅值在時移傅里葉分析中形成時幅拐點之前，首先分析在一信號序列上線性疊加另一個頻率為k的信號序列后，該合成序列中頻率k所對應(yīng)的幅值變化情況.

現(xiàn)假設(shè)在(T1，T2)時間內(nèi)，有一采樣頻率為Fs的信號序列q[x]={q[0]，q[1]，q[2]，…，q[N-1]}，在z時刻，q[x]上線性疊加另一序列p[x]={p[0]，p[1]，p[2]，…，p[M-1]}，其中，0≤Z≤N-M+1，M≤N.則合成后的信號序列f[x]如下式：

f[x]={q[0],q[1],…,q[Z]+p[0],q[Z+1]+

p[1],…,q[Z+M-1]+p[M-1],…,

q[N-1]}

(1)

由離散傅里葉變換可知頻率k在原序列q[x]中的變換表示為式(2)，其中N為原序列q[x]的總采樣點數(shù)，記其變換系數(shù)為E，則有

(2)

同理，頻率k在序列p[x]中的離散傅里葉變換系數(shù)表示為式(3)，其中M為突變信號q[x]的總采樣點數(shù)，記其變換系數(shù)為D，則有

(3)

那么，由式(1)可得頻率k在合成序列f[x]中的離散傅里葉變換系數(shù)F[k]表示為式(4)，記其變換系數(shù)為Y，則有

(4)

1.2 時幅拐點分析算法介紹

(5)

由式(5)可知E和D均為常數(shù)，其中E是頻率k在原始序列中的變換系數(shù).若原始序列不包含頻率k，則E=0，此時為頻率突變；若包含頻率k，則E為一個正數(shù)，該值只能上下平移時幅曲線，不影響曲線的變化趨勢及時幅拐點的對應(yīng)時刻，此時為幅值突變.下面針對突變信號出現(xiàn)的不同情況進行相應(yīng)分析，并找出其時幅拐點與突變信號發(fā)生時間的對應(yīng)關(guān)系.

設(shè)突變信號進入到原始信號的時刻為TZ，消失時刻為Tm.則突變信號出現(xiàn)的時間可能有以下3種情況.

(1)TZ=T1且T1

此時由式(5)可得合成序列下頻率k的幅值與遞進時間之間的關(guān)系式：

E，Tm

(6)

由上式得出頻率k的時幅曲線如圖1所示，該曲線可明顯看出一個時幅拐點出現(xiàn)在Tm時刻.在該時刻，頻率幅值由遞減狀態(tài)進入常值狀態(tài)，可判斷出該頻率信號的消失時刻.

圖1 頻率k的時幅曲線(TZ=T1且T1

(2)T1

此時由式(5)可得

(7)

由上式得出頻率k的時幅曲線如圖2所示，該曲線可明顯看出兩個時幅拐點：一個為曲線峰值，出現(xiàn)在TZ時刻；一個為曲線遞減進入常值時刻Tm.由此可判斷出該頻率信號的進入時刻和消失時刻.

圖2 頻率k的時幅曲線(T1

(3)T1

此時由式(5)可得

(8)

由上式得出頻率k的時幅曲線如圖3所示，該曲線可明顯看出一個時幅拐點出現(xiàn)在TZ時刻.在該時刻，頻率幅值由遞增狀態(tài)進入常值狀態(tài)，可判斷出該頻率信號的出現(xiàn)時刻.

圖3 頻率k的時幅曲線(T1

2 數(shù)學(xué)模型

以Δt為時間步長，在(T1+nΔt，T2)時間段內(nèi)逐次取信號序列fj的傅里葉變換中頻率k的幅值A(chǔ)kn，則有

(9)

其中L為分析點數(shù)，L=N-nΔtFs；N為原序列的總采樣點數(shù)，N=(T2-T1)Fs，n≤N.

實現(xiàn)該算法的流程圖如圖4所示.

圖4 算法流程圖Fig.4 Flow chart of the algorithm

3 仿真校驗

3.1 理論仿真

3.1.1 仿真函數(shù)f1

該函數(shù)采樣頻率為12 800 Hz，采樣時間50 s.f1為以下4個正弦函數(shù)的疊加：

(10)

對f1的采樣信號進行時移傅里葉分析，結(jié)果如圖5所示.

圖5 f1信號的時幅曲線Fig.5 Amplitude-time curve of signal f1

圖5結(jié)果表明：①664 Hz信號的時幅曲線表現(xiàn)為幅值突變，符合上述圖2情況，該頻率的時幅曲線存在一個峰值時刻t=10 s，一個曲線遞減進入常值時刻t=15 s，說明該信號的進入時刻為10 s，消失時刻為15 s；曲線最終進入一個常數(shù)狀態(tài)，幅值為2，說明在該信號中始終存在一個幅值為2的664 Hz的正弦信號.②1 195 Hz信號的時幅曲線最終進入幅值為0狀態(tài)，表現(xiàn)為頻率突變，同樣符合上述圖2情況，存在一個峰值時刻t=10 s，一個曲線遞減進入常值時刻t=12 s，說明該頻率信號的進入時刻為10 s，消失時刻為12 s.③1 328 Hz信號的時幅曲線仍符合圖2情況，可判斷出該信號的出現(xiàn)時刻為20 s，消失時刻為22 s.該結(jié)果準確地捕捉到了上述第2種情況下突變信號的進入時刻和消失時刻.

3.1.2 仿真函數(shù)f2

該函數(shù)采樣頻率為12 800 Hz，采樣時間為50 s.f2為以下3個正弦函數(shù)的疊加：

(11)

對f2的采樣信號進行時移傅里葉分析，結(jié)果如圖6所示.

圖6 f2信號的時幅曲線Fig.6 Amplitude-time curve of signal f2

圖6結(jié)果表明：①1 195 Hz信號的時幅曲線符合圖3情況，存在一個幅值遞增進入常值時刻t=10 s，說明該信號的出現(xiàn)時刻為10 s；②1 328 Hz信號的時幅曲線同樣符合圖3情況，存在一個幅值遞增進入常值時刻t=20 s，說明該信號的出現(xiàn)時刻為20 s.該結(jié)果準確地捕捉到了上述第3種情況下突變信號的進入時刻.

3.1.3 仿真函數(shù)f3

該函數(shù)采樣頻率為12 800 Hz，采樣時間為50 s.f3為以下兩個正弦函數(shù)的疊加：

(12)

對f3的采樣信號進行時移傅里葉分析，結(jié)果如圖7所示.

圖7結(jié)果表明：1 195 Hz信號的時幅曲線符合圖1情況，存在一個幅值遞減進入常值時刻t=20 s，說明該信號的消失時刻為20 s.該結(jié)果準確地捕捉到了上述第1種情況下正弦突變信號的消失時刻.

圖7 f3信號的時幅曲線Fig.7 Amplitude-time curve of signal f3

3.2 工程仿真

根據(jù)參考資料[9],可以構(gòu)造以下函數(shù)來仿真軸承外圈單點故障特征頻率時域波形:

(13)

圖8 外圈單點故障特征波形Fig.8 Characteristic waveform of single point fault of outer ring

對該函數(shù)進行128 000 Hz采樣，采樣時間為1 s，其頻譜圖如圖9所示.從圖上可看出，基頻3 585 Hz的幅值為0.1.

圖9 外圈單點故障頻譜圖Fig.9 Single point fault spectrum of outer ring

采用某火箭發(fā)動機高速滾動軸承正常試驗振動數(shù)據(jù)，將仿真波形采樣數(shù)據(jù)從第378 952個數(shù)據(jù)點疊加到真實的試驗數(shù)據(jù)上，構(gòu)成一個新的數(shù)據(jù)序列.真實數(shù)據(jù)與合成數(shù)據(jù)的波形如圖10和圖11所示，從圖中可以看出，在時域上很難發(fā)現(xiàn)該信號，而直接進行頻譜分析又會大大削弱該頻率信號的幅值.

針對外圈單點故障特征波形的基頻3 585 Hz,采用文中提出的分析方法對合成數(shù)據(jù)進行定頻分析，分析結(jié)果如圖12所示，顯然符合上述圖2情況，可以找出該頻率出現(xiàn)和消失的時點分別為2.983 61 s和3.983 61 s.

圖10 真實試驗數(shù)據(jù)的時域圖Fig.10 Time domain figure of real test data

圖11 合成數(shù)據(jù)的時域圖Fig.11 Time domain figure of synthetic data

圖12 合成數(shù)據(jù)的時幅曲線Fig.12 Amplitude-time curve of synthetic data

在找出頻率的存在時間后，在該時間區(qū)域內(nèi)進行二次頻譜分析便可得到該信號的實際幅值.根據(jù)圖12的結(jié)果可知，3 585 Hz信號出現(xiàn)在2.983 61 s到3.983 61 s之間.對2.983 61 s到3.983 61 s之間的數(shù)據(jù)進行頻譜分析則可得到信號的實際幅值，結(jié)果如圖13所示，從圖中可以看出該信號的幅值為0.1，與圖9添加前的信號幅值相等.

圖13 合成數(shù)據(jù)的區(qū)間分析Fig.13 Interval analysis of synthetic data

4 與小波分析的比較

用Matlab中“cmor5- 5”對3.2節(jié)的工程數(shù)據(jù)采用‘db4’三層分解，取第三層低頻段數(shù)據(jù)進行小波變換，得出的時頻圖如圖14所示.

由圖14可見，在2.961～3.957 s時間段內(nèi)出現(xiàn)了一個在3 250～3 750 Hz頻帶范圍內(nèi)的頻率.將該時間段進行局部放大，可以看出該條紋呈齒輪狀，說明信號幅值是周期變化的.

與時幅拐點分析方法相比，兩者都準確地捕獲到突發(fā)信號的出現(xiàn)和消失時刻，區(qū)別在于時幅拐點分析的是一個頻率，小波分析的是一個頻帶.在突發(fā)信號的幅值確定方面，小波分析方法并不能直接獲取該突發(fā)信號的幅值.圖中通過顏色來區(qū)分的僅代表該突發(fā)信號與小波基函數(shù)的相關(guān)程度，不代表頻率幅值情況.

圖14 小波分析頻譜圖Fig.14 Wavelet transform spectrum

5 結(jié)論

針對時域信號中突變信號難以準確定位的問題，提出一種突變信號的時幅拐點分析算法.該算法彌補了傳統(tǒng)頻譜分析在局部時域上的局限性，與短時傅里葉變換和小波變換相比，它以充足的信息量作為傅里葉變換的基礎(chǔ)，對單頻率進行時移傅里葉分析，獲得其時幅曲線，通過時幅曲線拐點能較為準確地判別特征頻率信號的出現(xiàn)與消失時刻，對分析時段內(nèi)的信號進行二次傅里葉分析可獲得較為準確的信號幅值.最后通過理論分析、理論仿真和工程仿真驗證了文中模型的精確性以及時幅拐點算法在捕捉突變信號中的有效性.

[1] Ingrid Daubechies.小波十講 [M].李健平,譯.北京：國防工業(yè)出版社，2011：1- 5.

[2] 隆軍,吳金強.STFT和HHT在風(fēng)力機軸承故障診斷中的應(yīng)用 [J].噪聲與振動控制,2013(4):219- 222.

LONG Jun,WU Jin-qiang.Application of short time fourier transform and hilbert-huang transform in fault diagnosis of rolling bearing of windmill [J].Noise and Vibration Control,2013(4):219- 222.

[3] 蘇文勝,王奉濤,朱泓,等.基于小波包樣本熵的滾動軸承故障特征提取 [J].振動、測試與診斷,2011,31(2):162- 166.

SU Wen-sheng,WANG Feng-tao,ZHU Hong,et al.Feature extraction rolling element bearing fault using wavelet packet sample entropy [J].Journal of Vibration Measurement & Diagnosis,2011,31(2):162- 166.

[4] 蘇濤,夏均忠,李樹珉,等.基于小波降噪和希爾伯特黃變換的滾動軸承故障特征提取 [J].軍事交通學(xué)院學(xué)報,2014,16(3):52- 56.

SU Tao,XIA Jun-zhong,LI Shu-min,et al.Fault features extraction of rolling bearing based on wavelet De-noising and HHT [J].Journal of Military Transportation University,2014,16 (3):52- 56.

[5] 張云鵬,蓋強.S變換在滾動軸承故障診斷上的應(yīng)用 [J].應(yīng)用科技,2011,38(7):26- 28.

ZHANG Yun-peng,GAI Qiang.Application of S-Transform in rolling bearing fault diagnosis [J].Applied Science and Technology,2011,38(7):26- 28.

[6] 陳換過,易永余,陳文華,等.基于廣義S變換的齒輪箱軸承故障診斷方法 [J].中國機械工程,2017,28(1):51- 56.

CHEN Huan-guo,YI Yong-yu,CHEN Wen-hua,et al.Fault diagnosis method of gearbox bearings based on gene-ralized S-transform [J].China Mechanical Engineering,2017,28(1): 51- 56.

[7] 張云強,張培林,吳定海,等.基于最優(yōu)廣義S變換和脈沖耦合神經(jīng)網(wǎng)絡(luò)的軸承故障診斷 [J].振動與沖擊,2015,34(9):26- 31.

ZAHNG Yun-qiang,ZHANG Pei-lin,WU Ding-hai,et al.Bearing fault diagnosis based on optimal generalized S transform and pulse couped neural network [J].Journal of Vibration and Shock,2015,34(9):26- 31.

[8] 唐先廣,郭瑜,丁彥春,等.基于短時傅里葉變換和獨立分量分析的滾動軸承包絡(luò)分析 [J].機械強度,2012,34(1):5- 9.

TANG Xian-guang,GUO Yu,DING Yan-chun,et al.App-lication of rolling element bearing envelope analysis based on short time fourier transition and Independent components analysis [J].Journal of Mechanical Strength,2012,34(1):5- 9.

[9] 何衛(wèi)鋒,程禮,李應(yīng)紅.某雙轉(zhuǎn)子發(fā)動機沖擊振動的分頻小波提取與故障分析 [J].航空動力學(xué)報,2006,21(2):361- 365.

HE Wei-feng,CHENG Li,LI Ying-hong.Wavelet extraction of Impact vibration frequency component and fault analysis for a twine spool aero-engine [J].Journal of Aerospace Power,2006,21(2):361- 365.

[10] 陳東雷.高速牽引電機軸承試驗與測試研究 [D].北京：北京交通大學(xué),2009.

[11] 吳德華,魯五一,熊紅云,等.基于RBF神經(jīng)網(wǎng)絡(luò)的航空發(fā)動機軸承故障診斷 [J].長沙航空職業(yè)技術(shù)學(xué)院學(xué)報,2005,5(2):23- 26.

WU De-hua,LU Wu-yi,XIONG Hong-yun,et al.Aeroplane engine bearing fault testing by RBF neural network [J].Journal of Changsha Aeronautical Vocational & Technical College,2005,5(2):23- 26.

[12] 吳勇.基于小波分析的高速牽引電機軸承故障診斷研究 [D].北京：北京交通大學(xué),2011.

[13] 郭遠晶,魏燕定,周曉軍,等.S變換用于滾動軸承故障信號沖擊特征提取 [J].振動、測試與診斷,2014,34(5):818- 822.

GUO Yuan-jing,WEI Yan-ding,ZHOU Xiao-jun,et al.Impact feature extraction from rolling bearing fault signal by S transform [J].Journal of Vibration Measurement & Diagnosis,2014,34(5):818- 822.

[14] 徐宏武,范玉剛,吳建德.基于S變換特征提取和隱馬爾科夫模型的故障診斷方法研究 [J].計算機與應(yīng)用化學(xué),2016,33(2):152- 156.

XU Hong-wu,FAN Yu-gang,WU Jian-de.The research on fault diagnosis method based on feature extraction by S-transform and hidden Markov Model [J].Computers and Applied Chemistry,2016，33(2):152- 156.

[15] 楊國安.滾動軸承故障診斷實用技術(shù) [M].北京： 中國石化出版社，2012：32- 54.

AMethodtoAnalyzeAmplitude-TimeInflectionPointofMutationSignals

DAIYi-mei1ZHANGHe-sheng1LIDong2QIHong-mei1FANGKe1

(1. School of Electrical Engineering, Beijing Jiaotong University, Beijing 100044, China; 2. School of Software Engineering, South China University of Technology, Guangzhou 510640, Guangdong, China)

As the mutation signals hidden in the vibration data of engine’s high-speed bearing are difficult to pinpoint in time domain, a new algorithm to obtain the amplitude-time inflection point of mutation signals is proposed. Firstly, the characteristic frequency of the signal is extracted via discrete Fourier transform. Secondly, a time-lapse Fourier analysis of the frequency is executed, which helps obtain amplitude-time inflection points. Then, the appearing and disappearing moments of the mutation signal are discovered according to the inflection points. Finally, the amplitude of the mutation signal is obtained via the quadratic Fourier analysis of the signal in a certain period. The results of theoretical analysis, theoretical simulation and engineering simulation all indicate that the proposed algorithm can capture the amplitude of mutation signals accurately, and obtain the appearing and disappearing moments of mutation signals in vibration data.

mutation signal; Fourier transform; amplitude-time inflection point; frequency analysis

2016- 10- 28

廣東省科技計劃項目(2015B010103002，2016B09098062);北京航天動力研究所委托項目(E14GY500010)

*Foundationitems： Supported by the Science and Technology Plan of Guangdong Province(2015B010103002，2016B09098062)

戴屹梅(1970-)，女，博士生，高級工程師，主要從事信號采集、數(shù)值分析、高速軸承故障診斷等研究.E-mail:daiym@spacechina.com

1000- 565X(2017)07- 0077- 07

TH 133.3

10.3969/j.issn.1000-565X.2017.07.011