基于HWPT-ZFFT的二維全息譜計(jì)算方法*

李紀(jì)永， 李舜酩， 陳曉紅， 王 勇

(1.南京航空航天大學(xué)能源與動(dòng)力學(xué)院 南京，210016) (2.南京航空航天大學(xué)理學(xué)院 南京，210016)

基于HWPT-ZFFT的二維全息譜計(jì)算方法*

李紀(jì)永1， 李舜酩1， 陳曉紅2， 王 勇1

(1.南京航空航天大學(xué)能源與動(dòng)力學(xué)院 南京，210016) (2.南京航空航天大學(xué)理學(xué)院 南京，210016)

精確的頻率、相位和幅值識(shí)別是進(jìn)行全息譜計(jì)算的必備條件，針對(duì)含3個(gè)及以上的密集頻譜成分，提出一種基于諧波小波包變換的頻譜細(xì)化方法(harmonic wavelet packets transform-zoom fast Fourier transform,簡(jiǎn)稱HWPT-ZFFT)，較傳統(tǒng)的復(fù)調(diào)制細(xì)化傅里葉變換所利用的低通及帶通濾波器相比，其盒型頻譜特性可將感興趣頻段的信號(hào)正交，無(wú)冗余、無(wú)泄漏地提取出，提高了識(shí)別精度。首先，利用諧波小波包對(duì)密集頻譜成分進(jìn)行濾波；然后，頻移進(jìn)而重采樣，進(jìn)行傅里葉變換得到細(xì)化的頻率、幅值及相位；最后，計(jì)算密集頻率下二維全息譜，進(jìn)行雙盤轉(zhuǎn)子全息譜計(jì)算，考慮高次分倍頻，得到更豐富的故障特征。仿真及雙盤轉(zhuǎn)子實(shí)驗(yàn)結(jié)果表明所提出方法的有效性。

諧波小波包； 密集頻譜； 頻譜細(xì)化； 全息譜

引 言

傳統(tǒng)的譜分析應(yīng)用在旋轉(zhuǎn)機(jī)械故障診斷有兩個(gè)缺陷：a.相位與振幅分離，甚至忽略相位信息；b.轉(zhuǎn)子截面兩向振動(dòng)信號(hào)關(guān)系未能體現(xiàn)出來(lái)。從軸心軌跡的形狀、穩(wěn)定性和旋轉(zhuǎn)方向等方面進(jìn)行分析，得到比較全面的機(jī)組運(yùn)行狀態(tài)信息，發(fā)展了基于信息融合的全譜、全息譜和全矢譜技術(shù)[1]。二維全息譜將轉(zhuǎn)子在一個(gè)截面的水平和垂直的振動(dòng)幅值與相位繪制譜圖，反映了轉(zhuǎn)子在一個(gè)支承面的振動(dòng)情況[2-3]，在實(shí)踐中得到廣泛應(yīng)用并積累了許多特征譜圖[4-5]。

由于全息譜技術(shù)是建立在傅里葉變換基礎(chǔ)之上，精確的頻率、幅值與相位是全息譜計(jì)算的關(guān)鍵，經(jīng)典的全息譜方法利用比值內(nèi)插法校正頻率，但難以區(qū)分密集頻譜成分。密集頻譜成分細(xì)化與校正方法可分為兩類：a.包含兩個(gè)密集頻率成分信號(hào)的頻譜細(xì)化方法;b.包含3個(gè)及以上密集頻率成分信號(hào)的頻譜細(xì)化方法。含3個(gè)及以上的密集頻率成分細(xì)化思路有：進(jìn)行頻移、濾波和重采樣提高頻率分辨率；在窗長(zhǎng)度不變的情況下進(jìn)行局部細(xì)化；拓展號(hào)序列空間，增加數(shù)據(jù)長(zhǎng)度等[6]。文獻(xiàn)[7]提出復(fù)解析帶通濾波器的復(fù)調(diào)制細(xì)化選帶方法，采用復(fù)解析帶通濾波器濾波，具有速度快精度高特點(diǎn)。文獻(xiàn)[8]利用基于小波變換的頻譜細(xì)化方法細(xì)化密集頻譜成分，即在進(jìn)行細(xì)化處理之前利用組合的高斯小波進(jìn)行濾波，在組合高斯小波邊緣依然不光滑，并指出盒型頻譜其濾波性能好[9]。

筆者提出HWPT-ZFFT方法，首先，利用諧波小波包進(jìn)行濾波[10-11]，其盒型頻譜特性可將信號(hào)正交、無(wú)冗余和無(wú)遺漏的分解在相應(yīng)的頻帶上[12-13]；然后，細(xì)化密集頻譜成分，得到修正的頻率、幅值和相位；最后，進(jìn)行二維全息譜分析。由于考慮了密集頻譜成分，所以體現(xiàn)出豐富的故障特征。

1 改進(jìn)的二維全息譜

二維全息譜具體實(shí)現(xiàn)過(guò)程如下:a.將一個(gè)支承截面上的x，y兩方向的振動(dòng)信號(hào)進(jìn)行快速傅里葉變換得到頻率幅值相位；b.將變換之后的幅值與相位以階次為單位進(jìn)行融合，根據(jù)幅值和相位信息計(jì)算出每個(gè)階次處的離心率、長(zhǎng)軸和短軸等橢圓信息。傳統(tǒng)的全息譜高頻只包含整數(shù)倍頻成分。在實(shí)際故障信號(hào)中不可忽略高次諧頻，筆者將幅值明顯的諧倍頻率考慮進(jìn)去，進(jìn)行全息譜計(jì)算。

含n個(gè)階次頻率的兩路互相垂直信號(hào)為

(1a)

(1b)

其中：ω為頻率；X與Y為信號(hào)幅值；φx與φy為相位。

不失一般性，令n=1，展開(kāi)可得

(2)

將式(1)展開(kāi)，消去t，得到全息譜方程為

(3)

式(3)為階次頻率對(duì)應(yīng)的橢圓軌跡方程，示意圖如圖1所示。

圖1 二維全息譜示意圖Fig.1 The sketch of 2-D holospectrum

根據(jù)幾何關(guān)系，其橢圓長(zhǎng)半軸Ra、短半軸Rb、長(zhǎng)半軸與x軸傾角α及轉(zhuǎn)子中心沿橢圓軌跡運(yùn)動(dòng)的相位角為φα。

(4)

由以上分析可知，構(gòu)成全息譜的3個(gè)參量為階次頻率、相位及其幅值。精確的階次頻率，相位及幅值是構(gòu)成準(zhǔn)確全息譜的必要條件，筆者利用諧波小波包分解信號(hào)，然后利用密集頻譜法對(duì)信號(hào)頻譜、相位及幅度進(jìn)行校正。

2 基于諧波小波包頻譜細(xì)化方法

諧波小波具有零相移濾波器，可以鎖定相位且諧波小波變換在不同尺度下數(shù)據(jù)點(diǎn)數(shù)不變，可精確進(jìn)行頻域定位。諧波小波與小波類似，其頻譜帶寬以二進(jìn)的方式擴(kuò)大，不能夠任意選擇頻段，而諧波小波包則可將信號(hào)分解到相互獨(dú)立的任意頻段上且各頻段分辨率一致，如圖2所示。

圖2 諧波小波變換后頻率分布示意圖Fig.2 Signal HWPT frequency distribution

分析頻帶帶寬為(n-m)2π，分析時(shí)間中心在t=k/(n-m)處的諧波小波的一般形式為

(5)

頻域形式為

(6)

信號(hào)f(t)的諧波小波變換為

(7)

(8)

式(7)與式(8)為信號(hào)f(t)諧波小波變換(harmonic wavelet packets transform,簡(jiǎn)稱HWPT)后時(shí)域及頻域表達(dá)式。

設(shè)s(t)為諧波小波包分解后某一頻段信號(hào)，進(jìn)行平移、重采樣后再做傅里葉變換，即可得到細(xì)化頻率，其基本流程如下。

1) 濾波：按實(shí)際選取頻段頻率區(qū)間為[flfh]，則其中心頻率可表示為f0=(fl+fh)/2。

其思路與復(fù)調(diào)制細(xì)化傅里葉變換及基于傳統(tǒng)小波變換頻譜細(xì)化方法一致，不同之處在于筆者采用諧波小波包進(jìn)行濾波，根據(jù)需要設(shè)置分解層數(shù)。鑒于信號(hào)噪聲往往存在于高頻段，選取低頻段則降低了噪聲對(duì)細(xì)化的影響。

3 仿真信號(hào)頻譜細(xì)化

構(gòu)造仿真信號(hào)中信號(hào)采樣頻率為fs=1 024;傅里葉變換點(diǎn)數(shù)為1 024，信號(hào)長(zhǎng)度點(diǎn)數(shù)為20 480。

f(t)=cos(63.774 3t+0.523 6)+cos(64.406t+

0.698 1)+cos(66.915 9t+1.047 2)+

cos(125.663 7t+0.733)+cos(251.327 4t)+

(9)

其中：n(t)為噪聲；信噪比為3.9。

仿真時(shí)間歷程與傅里葉變換如圖3所示。由于傅里葉變換時(shí)分辨率為1 Hz，密集頻譜成分并未顯示出。

圖3 仿真時(shí)間歷程與頻域幅值Fig.3 Simulation signal time history and spectrum graph

利用線性調(diào)頻Z變換(chirp Z transform，簡(jiǎn)稱CZT)細(xì)化，結(jié)果如圖4(a)所示。復(fù)調(diào)制ZFFT細(xì)化，分別利用低通及帶通濾波，結(jié)果如圖4(b)及圖4(c)所示。諧波小波包濾波變換后(分解層數(shù)為4，頻段間隔為32 Hz)細(xì)化結(jié)果如圖4(d)所示。可以看出，帶通濾波時(shí)出現(xiàn)虛假頻率成分，這是由于濾波器設(shè)計(jì)不當(dāng)造成頻譜泄漏，從而造成細(xì)化頻譜失真。比較去噪效果可知，HWPT去噪效果優(yōu)于帶通及低通濾波，而CZT方法由于未涉及濾波，故無(wú)去噪效果。表1為細(xì)化后的頻率、幅值及相位。結(jié)合圖4可知，諧波小波包校準(zhǔn)精度最高，這說(shuō)明了諧波小波包濾波及細(xì)化的有效性及優(yōu)越性。

圖4 細(xì)化結(jié)果Fig.4 Zoom results

表1 頻率、幅值及相位細(xì)化結(jié)果

Tab.1 Detail zoom results of frequency, amplitude and phase

項(xiàng)目f/Hz幅值相位/(°)原信號(hào)10.1501.00030.00010.2501.00040.00010.6501.00060.000低通濾波10.1500.97729.72810.2500.99639.14110.6501.10160.244CZT10.1401.05056.70910.2501.04025.43110.6500.99548.797HWPT10.1500.98029.72810.2500.99839.21410.6501.10260.062

4 實(shí)驗(yàn)分析

圖5為雙盤轉(zhuǎn)子實(shí)驗(yàn)裝置示意圖。采用柔性聯(lián)軸器連接，利用渦流傳感器采集位移信號(hào)。其渦流傳感器水平方向與x軸夾角為0，如圖6所示，并按觀察傳感器方向觀察軸心軌跡。

雙跨轉(zhuǎn)子含碰摩及不對(duì)中故障，信號(hào)采樣頻率為2 048 Hz，點(diǎn)數(shù)為32 768。對(duì)其進(jìn)行傅里葉變換，結(jié)果如圖7所示，傅里葉變換點(diǎn)數(shù)為2 048，密集頻率成分沒(méi)有計(jì)算出來(lái)，采用HWPT-ZFFT算法，分解層數(shù)設(shè)為4，細(xì)化倍數(shù)設(shè)為30。對(duì)x方向基頻及二倍頻細(xì)化結(jié)果如圖8所示。細(xì)化的具體頻率、幅值與相位如表2所示。y方向細(xì)化后的頻率、幅值與相位如表3所示。利用傳統(tǒng)的二維全息譜方法計(jì)算，得到倍頻及低次分倍頻成分下的二維全息譜如圖9所示。橢圓具體參數(shù)如表4所示。橢圓軌跡主要在基頻及二倍頻成分上呈現(xiàn)出復(fù)合故障，旋向不一致表明了不對(duì)稱故障。

圖5 實(shí)驗(yàn)裝置示意圖Fig.5 Experimental device schematic diagram

圖6 渦流傳感器安裝示意圖Fig.6 Eddy current sensor installition schematic diagram

圖7 故障信號(hào)頻譜Fig.7 Fault signal spectrum

圖8 細(xì)化頻率Fig.8 Zoom frequency

f/Hz幅值/(mV·Hz-1)相位/(°)50.65024.55010.53051.41015.220152.900100.20022.540-110.854101.30011.920-149.597102.9008.65113.798

表3 y方向的頻率、幅值與相位

Tab.3 y direction signal frequency, amplitude and phase

f/Hz幅值/(mV·Hz-1)相位/(°)50.6541.840170.13151.4130.10048.772100.23.99864.818101.314.51-16.0445102.98.943-125.3941

圖9 倍頻及低次分倍頻二維全系譜圖Fig.9 Traditional holospectrum graph

階次長(zhǎng)軸短軸無(wú)量綱比值離心率長(zhǎng)軸傾角旋轉(zhuǎn)方向1-19.07033.950-68.4700.85786.400負(fù)向2-12.93014.87025.6200.927124.800正向31.180-1.6500.8100.888121.400負(fù)向40.8601.9500.6300.99749.600負(fù)向0.05-0.530-0.2903.1200.99096.400負(fù)向0.58-1.1501.2300.1900.954114.500負(fù)向0.951.1102.74013.2500.98392.100正向0.981.020-0.770-6.1300.98797.100正向

高次分倍頻二維全息譜如圖10所示，包括1.015次倍頻，1.98次倍頻，2.03次倍頻。橢圓軌跡參數(shù)如表5所示。1.98次倍頻的橢圓軌跡扁平，表示發(fā)生碰摩故障。高次分倍頻全息譜的計(jì)算豐富了故障信息，綜合倍頻、低次分倍頻及高次分倍頻的結(jié)果可知轉(zhuǎn)子發(fā)生了碰摩及不對(duì)中故障。

5 結(jié) 論

1) 利用諧波小波包將信號(hào)正交，無(wú)冗余地分解在所需頻段上，其對(duì)噪聲抑制性能優(yōu)于帶通及低通濾波。

圖10 高次分倍頻二維全息譜Fig.10 High order fractional frequency holospectrum graph

Tab.5 Holospectrum information of high order fractional frequency

階次長(zhǎng)軸短軸無(wú)量綱比值離心率長(zhǎng)軸傾角旋轉(zhuǎn)方向1.0153.472-6.77111.3200.87799.200負(fù)向1.98-9.471-3.7121.8100.987169.200正向2.030.9703.961-3.6500.929131.500正向

2) 對(duì)HWPT分解后的信號(hào)利用ZFFT方法細(xì)化頻譜，得到校正的密集頻譜成分，與傳統(tǒng)的復(fù)調(diào)制ZFFT及線性調(diào)頻CZT變換相比，其頻率、幅值及相位校正精度得到提高，表明了所提出HWPT-ZFFT方法的有效性。

3) 考慮細(xì)化后密集的頻譜進(jìn)行二維全息譜計(jì)算，以及雙跨轉(zhuǎn)子實(shí)驗(yàn)計(jì)算中考慮高次諧倍頻成分，能夠得到更豐富的故障特征，表明考慮高次分頻進(jìn)行二維全息譜計(jì)算的必要性。

[1] Goldman P, Muszynska A. Application of full spectrum to rotating machinery diagnostics[J]. Orbit,1999,20(1): 17-21.

[2] Qu L, Liu X, Peyronne G, et al. The holospectrum: a new method for rotor surveillance and diagnosis[J]. Mechanical Systems and Signal Processing, 1989, 3(3): 255-267.

[4] 杜永祚,秦志英.旋轉(zhuǎn)機(jī)械動(dòng)態(tài)信號(hào)全息譜分析[J].振動(dòng)、測(cè)試與診斷,2002，22(2):81-88.

Du Yongzuo, Qin Zhiying. Holo-spectrum analysis of rotating machinery dynamic signals[J]. Journal of Vibration, Measurement & Diagnosis, 2002,22(2): 81-88. (in Chinese)

[5] 屈梁生.機(jī)械故障的全息診斷原理[M].北京：科學(xué)出版社，2007:156-174.

[6] 李宏坤, 周帥, 朱泓, 等. 基于經(jīng)驗(yàn)?zāi)Ｊ椒纸獾娜⒆V故障識(shí)別方法[J]. 振動(dòng)、測(cè)試與診斷, 2011, 31(1): 20-22.

Li Hongkun, Zhou Shuai, Zhu Hong, et al. Rotor fault classification using empirical mode decomposition based holospectrum[J]. Journal of Vibration, Measurement & Diagnosis, 2011, 31(1):20-22. (in Chinese)

[7] 毛育文, 涂亞慶，肖瑋，等. 離散密集頻譜細(xì)化分析與校正方法研究進(jìn)展[J]. 振動(dòng)與沖擊, 2012,31(21): 112-119.

Mao Yuwen, Tu Yaqing, Xiao Wei, et al. Advances and trends of study on discrete intensive frequency spectrum zooming analysis and correction methodology[J]. Journal of Vibration and Shock, 2012, 31(21):112-119. (in Chinese)

[8] 馬建倉(cāng), 吳啟彬. 基于小波變換的頻譜細(xì)化分析方法[J]. 信號(hào)處理, 1997, 13(3): 276-279.

Ma Jiancang, Wu Qibin. Spectrum zoom analysis method based on wavelet transform [J]. Signal Processing,1997, 13(3): 276-279. (in Chinese)

[9] 曹豫寧, 李永麗, 梅云, 等. 基于小波變換的頻譜細(xì)化方法在電動(dòng)機(jī)故障檢測(cè)中的應(yīng)用[J].繼電器, 2002, 30(6): 1-3.

Cao Yuning, Li Yongli, Mei Yun, et al. Diagnosis of rotor faults of inductive motors using spectrum zoom analysis method based on wavelet transform[J].Relay,2002,30(6):1-3. (in Chinese)

[10]李舜酩. 諧波小波包方法及其對(duì)轉(zhuǎn)子亞頻軸心軌跡的提取[J]. 機(jī)械工程學(xué)報(bào), 2004, 40(9): 133-137.

Li Shunming. Harmonic wavelet packetes method and used on accurate obtaining the orbit of rotor sub-frequency signal[J]. Chinese Journal of Mechanical Engineering, 2004,40(9): 133-137. (in Chinese)

[11]Yan R, Gao R X. An efficient approach to machine health diagnosis based on harmonic wavelet packet transform[J]. Robotics and Computer-Integrated Manufacturing, 2005, 21(4): 291-301.

[12]Jiang C, Weng X T, Lou J J. Gear fault diagnosis based on harmonic wavelet packet and BP neural network[J]. Applied Mechanics and Materials, 2012, 21(7): 2683-2687.

[13]Bouzida A. Fault diagnosis in industrial induction machines through discrete wavelet transform[J]. Industrial Electronics, IEEE Transactions on, 2011,58(9): 4385-4395.

10.16450/j.cnki.issn.1004-6801.2016.01.011

*機(jī)械結(jié)構(gòu)強(qiáng)度與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室開(kāi)放課題資助項(xiàng)目(SV2015-KF-01)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)基金資助項(xiàng)目(NZ2015103)

2014-01-15；修回日期：2014-03-03

V233.1； TH132

李紀(jì)永，男，1985年1月生，博士研究生。主要研究方向?yàn)樾D(zhuǎn)機(jī)械故障信號(hào)處理。 E-mail：ljynav@163.com