李紀永， 李舜酩， 陳曉紅， 王 勇

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

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

李紀永1， 李舜酩1， 陳曉紅2， 王 勇1

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

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

諧波小波包； 密集頻譜； 頻譜細化； 全息譜

引 言

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

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

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

1 改進的二維全息譜

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

含n個階次頻率的兩路互相垂直信號為

(1a)

(1b)

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

不失一般性，令n=1，展開可得

(2)

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

(3)

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

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

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

(4)

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

2 基于諧波小波包頻譜細化方法

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

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

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

(5)

頻域形式為

(6)

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

(7)

(8)

式(7)與式(8)為信號f(t)諧波小波變換(harmonic wavelet packets transform,簡稱HWPT)后時域及頻域表達式。

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

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

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

3 仿真信號頻譜細化

構(gòu)造仿真信號中信號采樣頻率為fs=1 024;傅里葉變換點數(shù)為1 024，信號長度點數(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。

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

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

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

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

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

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

項目f/Hz幅值相位/(°)原信號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 實驗分析

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

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

圖5 實驗裝置示意圖Fig.5 Experimental device schematic diagram

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

圖7 故障信號頻譜Fig.7 Fault signal spectrum

圖8 細化頻率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

階次長軸短軸無量綱比值離心率長軸傾角旋轉(zhuǎn)方向1-19.07033.950-68.4700.85786.400負向2-12.93014.87025.6200.927124.800正向31.180-1.6500.8100.888121.400負向40.8601.9500.6300.99749.600負向0.05-0.530-0.2903.1200.99096.400負向0.58-1.1501.2300.1900.954114.500負向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ā)生碰摩故障。高次分倍頻全息譜的計算豐富了故障信息，綜合倍頻、低次分倍頻及高次分倍頻的結(jié)果可知轉(zhuǎn)子發(fā)生了碰摩及不對中故障。

5 結(jié) 論

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

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

Tab.5 Holospectrum information of high order fractional frequency

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

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

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

[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)機械動態(tài)信號全息譜分析[J].振動、測試與診斷,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] 屈梁生.機械故障的全息診斷原理[M].北京：科學出版社，2007:156-174.

[6] 李宏坤, 周帥, 朱泓, 等. 基于經(jīng)驗模式分解的全息譜故障識別方法[J]. 振動、測試與診斷, 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] 毛育文, 涂亞慶，肖瑋，等. 離散密集頻譜細化分析與校正方法研究進展[J]. 振動與沖擊, 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] 馬建倉, 吳啟彬. 基于小波變換的頻譜細化分析方法[J]. 信號處理, 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] 曹豫寧, 李永麗, 梅云, 等. 基于小波變換的頻譜細化方法在電動機故障檢測中的應用[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]李舜酩. 諧波小波包方法及其對轉(zhuǎn)子亞頻軸心軌跡的提取[J]. 機械工程學報, 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

*機械結(jié)構(gòu)強度與振動國家重點實驗室開放課題資助項目(SV2015-KF-01)；中央高?；究蒲袠I(yè)務費專項基金資助項目(NZ2015103)

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

V233.1； TH132

李紀永，男，1985年1月生，博士研究生。主要研究方向為旋轉(zhuǎn)機械故障信號處理。 E-mail：ljynav@163.com