李紀(jì)永，李舜酩，陳曉紅，江星星

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

基于雙樹(shù)復(fù)小波變換的轉(zhuǎn)靜碰摩響應(yīng)分析

李紀(jì)永1，李舜酩1，陳曉紅2，江星星1

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

針對(duì)碰摩故障診斷中轉(zhuǎn)靜碰摩響應(yīng)消噪及特征提取構(gòu)造分段圓弧平滑閾值函數(shù)，利用小波、小波包及雙樹(shù)復(fù)小波（Dual-tree Complex Wavelet）分別對(duì)非平穩(wěn)信號(hào)降噪處理，以信噪比為指標(biāo)判斷降噪效果。結(jié)果表明，構(gòu)造的閾值函數(shù)消噪具有高效性。通過(guò)計(jì)算單轉(zhuǎn)子局部碰摩動(dòng)態(tài)響應(yīng)，用雙樹(shù)復(fù)小波對(duì)位移響應(yīng)進(jìn)行特征提取。與小波及小波包相比，雙樹(shù)實(shí)小波近似解析性及采樣插值具有互補(bǔ)特點(diǎn)，能完整提取動(dòng)態(tài)響應(yīng)特征，減少信息丟失。并用仿真及轉(zhuǎn)子碰摩實(shí)驗(yàn)驗(yàn)證該方法的有效性。

轉(zhuǎn)靜碰摩；降噪；動(dòng)態(tài)響應(yīng)；閾值函數(shù)；雙樹(shù)復(fù)小波

航空發(fā)動(dòng)機(jī)中轉(zhuǎn)靜碰摩有兩類，即轉(zhuǎn)動(dòng)部分碰摩部位不在轉(zhuǎn)軸本身如葉片等與靜止部分發(fā)生的碰摩及轉(zhuǎn)動(dòng)部分發(fā)生在轉(zhuǎn)軸本身的碰摩。碰摩故障通常表現(xiàn)為其它故障的間接結(jié)果［1－2］，如轉(zhuǎn)子質(zhì)量不平衡、熱彎曲、軸系不對(duì)中及轉(zhuǎn)子渦動(dòng)失穩(wěn)等均可能表現(xiàn)為轉(zhuǎn)靜碰摩。碰摩發(fā)生時(shí)一般為先有局部碰摩，即在一個(gè)周期內(nèi)轉(zhuǎn)子與定子發(fā)生一次或數(shù)次碰撞。沖擊產(chǎn)生接觸力及摩擦力，摩擦力將轉(zhuǎn)子轉(zhuǎn)動(dòng)能量轉(zhuǎn)化為橫向振動(dòng)同時(shí)伴有較復(fù)雜的振動(dòng)現(xiàn)象，在某些條件下會(huì)出現(xiàn)混沌運(yùn)動(dòng)。

通過(guò)小波及小波包分析碰摩非平穩(wěn)信號(hào)［3－4］可提取細(xì)節(jié)成分，但仍有3種缺陷限制其應(yīng)用［5］：①平移敏感性，即輸人信號(hào)平移會(huì)在小波分解系數(shù)中產(chǎn)生不可預(yù)計(jì)變化；②方向性差，只能提供3個(gè)空間方位；③相位信息缺乏。小波變換域的能量越集中說(shuō)明小波分析效果越好，而該能量集中程度與小波幅頻、相頻特性有關(guān)。小波濾波器非理想截止特性及分解過(guò)程隔點(diǎn)采樣會(huì)致信號(hào)分解、重構(gòu)時(shí)產(chǎn)生虛假頻率成分。

雙樹(shù)復(fù)小波變換（Dual-tree Complex Wavelet Transform，DTCWT）由Kingsbury［6］提出。Salesnick等［7］對(duì)雙樹(shù)復(fù)小波濾波器構(gòu)造進(jìn)行深入研究。雙樹(shù)復(fù)小波變換具有的良好抗頻帶混疊能力及平移不變特性可保證雙樹(shù)復(fù)小波變換在濾波器特性及周期復(fù)合沖擊特征檢測(cè)方面優(yōu)于經(jīng)典離散小波變換、第二代小波變換及經(jīng)驗(yàn)?zāi)Ｊ椒纸獾葧r(shí)頻域分析方法。

降噪可視為經(jīng)典的非參數(shù)回歸，現(xiàn)有函數(shù)估計(jì)方法包括樣條曲線擬合估計(jì)、核函數(shù)估計(jì)及傅里葉展開(kāi)估計(jì)等［8］。Donolo等［9－10］在小波變換基礎(chǔ)上提出的閾值降噪方法優(yōu)勢(shì)在于小波基函數(shù)空間多樣性。小波變換稀疏性及去相關(guān)性保證真實(shí)信號(hào)小波變換能量集中于有限的小波系數(shù)，而如何選取閾值及閾值函數(shù)為實(shí)現(xiàn)小波去噪算法中關(guān)鍵問(wèn)題［11］。閾值計(jì)算方法［12］有：stein無(wú)偏估計(jì)、固定閾值及啟發(fā)式閾值；閾值函數(shù)有硬閾值、軟閾值。硬閾值含間斷點(diǎn)造成重構(gòu)信號(hào)振蕩，重構(gòu)信號(hào)不光滑；軟閾值所得小波系數(shù)與原小波系數(shù)存在固定偏差，影響重構(gòu)信號(hào)逼近程度。為此，研究者構(gòu)建不同閾值函數(shù)，包括折衷閾值函數(shù)、指數(shù)閾值函數(shù)等，但大多與真實(shí)小波系數(shù)有一定偏差。

本文采用stein無(wú)偏估計(jì)計(jì)算閾值，利用構(gòu)造的閾值函數(shù)對(duì)雙樹(shù)復(fù)小波變換后非平穩(wěn)信號(hào)進(jìn)行降噪處理，計(jì)算單轉(zhuǎn)子碰摩響應(yīng)，并對(duì)碰摩位移響應(yīng)用雙樹(shù)復(fù)小波進(jìn)行特征提?。焕脝无D(zhuǎn)子及多轉(zhuǎn)子碰摩實(shí)驗(yàn)驗(yàn)證本文所提方法。

1 雙樹(shù)復(fù)小波分解

DTCWT由兩實(shí)值小波φh（t）與φg（t）構(gòu)成復(fù)小波，即

式中：φg（t）＝H［φh（t）］，H為Hilbert算子。

兩正交小波函數(shù)組成Hilbert變換對(duì)的充要條件［5］為兩低通濾波器滿足半采樣延遲條件。其幅值與相位存在關(guān)系為

構(gòu)造關(guān)鍵環(huán)節(jié)為使兩小波函數(shù)濾波器滿足式（2），構(gòu)造方法有線性正交法、共同因子法及Q平移方法（quarter shift）。Q平移法設(shè)計(jì)的濾波器為

式（3）精確滿足幅值條件，近似滿足相位條件，其相位表達(dá)式為

對(duì)稱中心在點(diǎn)n＝0．5（N－1）－0．25附近。利用Q平移方法設(shè)計(jì)的雙樹(shù)復(fù)小波變換整體為精確線性相位，由于φc（t）的解析性，DTCWT具有近似平移不變與抗頻帶混疊等優(yōu)良特性。前者有益于多重諧波信號(hào)檢測(cè)；后者對(duì)提取周期性沖擊特征非常有利。

2 雙樹(shù)復(fù)小波降噪閾值及閾值函數(shù)

利用小波降噪核心為閾值及閾值函數(shù)選擇。含噪信號(hào)表示為：s＝f＋se，e為服從標(biāo)準(zhǔn)正態(tài)分布的高斯白噪聲，其強(qiáng)度可用穩(wěn)健估計(jì)方法計(jì)算，即σ＝中值運(yùn)算利用stein無(wú)偏似然估計(jì)，將某次層系數(shù)平方按由小到大排列，獲得某向量W＝［w1，w2，…，wn］，其中w1≤…wn，n為小波系數(shù)個(gè)數(shù)。由此估算向量R＝［r1，r2，…，rn］，其中：

式中：0≤α≤1，不同a值對(duì)應(yīng)不同閾值函數(shù)，實(shí)用性強(qiáng)，缺點(diǎn)亦明顯，即具有間斷點(diǎn)且與真實(shí)系數(shù)有一定偏差。

為克服軟、硬閾值缺點(diǎn)，構(gòu)造分段圓弧平滑閾值函數(shù)，設(shè)圓1半徑r1，圓2半徑r2，幾何關(guān)系見(jiàn)圖1。由圖1看出，圓弧1半徑為r1，圓弧2半徑為r2，圓1圓心坐標(biāo)為O1（1－r1，r1），圓2圓心坐標(biāo)為O2（1＋r2，r1），半徑滿足幾何關(guān)系r1＋r2＝βt。其中0≤β≤1。

圖1 幾何關(guān)系示意圖Fig．1 Geometric relation schematic

分段圓弧平滑閾值函數(shù)可表示為

改進(jìn)的閾值函數(shù)在各點(diǎn)處一階可導(dǎo)，利于數(shù)學(xué)處理，避免信號(hào)重構(gòu)時(shí)出現(xiàn)振蕩。調(diào)節(jié)參數(shù)b及半徑r1即可得不同折衷閾值函數(shù)。b＝1時(shí)即可平滑過(guò)渡至真實(shí)小波系數(shù)。a＝0．5，b＝1時(shí)示意見(jiàn)圖2。

圖2 閾值函數(shù)Fig．2 Threshold function

利用典型非平穩(wěn)信號(hào)Heavy sine（初始信噪比6．4 dB）、Doppler（初始信噪比6．4 dB）、Bump（初始信噪比7．3 dB）、Blocks（初始信噪比9．6 dB）驗(yàn)證分段圓弧平滑閾值函數(shù)降噪的有效性。信噪比依次（每組SNR順序從左至右，先第一行后第二行）見(jiàn)表1。綜合比較可知，雙樹(shù)復(fù)小波去噪性能較好，從而驗(yàn)證分段光滑圓弧平滑閾值函數(shù)降噪的有效性。

表1 Heavy sine，Doppler，Bump，Blocks降噪后SNR值Tab.1 SNR of denoised Heavy sine Dopp ler，Bum p，Blocks

3 單轉(zhuǎn)子仿真碰摩響應(yīng)雙樹(shù)復(fù)小波分析

Jeffcott單轉(zhuǎn)子碰摩模型見(jiàn)文獻(xiàn)［2］，本文重設(shè)參數(shù)為

式中：M為轉(zhuǎn)子質(zhì)量；c為軸阻尼；k為軸剛度；ks為非線性剛度系數(shù)；f為摩擦系數(shù)；kc為定子徑向剛度。

計(jì)算轉(zhuǎn)子動(dòng)態(tài)位移響應(yīng)其分叉圖見(jiàn)圖3。令s＝ω／ωc，ωc為臨界轉(zhuǎn)速。分析s＝5．65時(shí)動(dòng)態(tài)響應(yīng)，利用Wolf方法求得lyapunov指數(shù)為0．022，時(shí)間歷程及龐加萊映射見(jiàn)圖4（a），響應(yīng)呈非線性；其頻譜見(jiàn)圖4（b），存在的高頻率成分未完全顯示，如s＝84、95．2、112，所取信號(hào)頻率成份較豐富，能較好反映轉(zhuǎn)子碰摩故障現(xiàn)象。

圖3 轉(zhuǎn)子系統(tǒng)分叉特性Fig．3 Rotor system bifurcation characteristic

圖4 混沌信號(hào)時(shí)域及頻域特征Fig．4 Chaos response time and frequency characteristic

圖5 不同分解方法所得頻譜Fig．5 Spectrum via different extraction methods

利用小波、小波包及雙樹(shù)復(fù)小波對(duì)混沌信號(hào)進(jìn)行分析，所得頻譜見(jiàn)圖5。由圖5（a）看出，小波分解各層存在嚴(yán)重頻譜泄漏；而圖5（b）中除頻譜泄漏外亦存在較嚴(yán)重頻率混疊，如第2層存在5個(gè)虛假頻率，第3層節(jié)點(diǎn)在高頻處存在虛假頻率，第4層在中高頻處存在虛假頻率，小波及小波包分解的頻率泄漏由濾波器非理想截止特性造成，小波包分解產(chǎn)生的虛假頻率由分解過(guò)程中下抽樣造成采樣率降低，從而導(dǎo)致頻率折疊及在重構(gòu)過(guò)程中隔點(diǎn)插零操作形成。由圖5（c）看出，第2層提取的響應(yīng)高頻成份包括s＝84，95．2，112，得益于波重構(gòu)過(guò)程中雙樹(shù)分解系數(shù)互補(bǔ)，減小信息丟失。各層幾乎不存在虛假頻率。此因雙樹(shù)復(fù)小波帶通頻率響應(yīng)只在正半頻率軸上（除第一層外），能減小頻率混疊。而由于濾波器非理想截止特性，第4層存在頻譜泄漏。

4 轉(zhuǎn)子實(shí)驗(yàn)故障信號(hào)雙樹(shù)復(fù)小波分析

選單轉(zhuǎn)子及三轉(zhuǎn)子碰摩位移響應(yīng)信號(hào)進(jìn)行分析，其中三轉(zhuǎn)子實(shí)驗(yàn)裝置示意圖見(jiàn)圖6。單轉(zhuǎn)子臨界轉(zhuǎn)速4 216 r／min、2 754 r／min時(shí)的碰摩位移及頻譜見(jiàn)圖7。由圖7看出，基頻45．9 Hz及2倍頻91．8 Hz清晰可辨，而若干高次倍頻成份淹沒(méi)在噪聲中。利用雙樹(shù)復(fù)小波進(jìn)行降噪處理，選構(gòu)造分段光滑閾值函數(shù)對(duì)信號(hào)重構(gòu)，頻譜見(jiàn)圖8。由圖8看出，已提取信號(hào)的高頻細(xì)節(jié)信息，其中第二層含505．1 Hz，593．2 Hz，643 Hz等；第三層含328．3 Hz（7倍頻），367．4 Hz等；第四層含183．8 Hz（4倍頻），204．1 Hz；第5層為2倍頻，頻帶間有輕微混疊，與相同分解層下小波包分解混疊頻譜成份相比（圖9，線框內(nèi)），頻率混疊成份大幅降低。

圖6 轉(zhuǎn)子裝置示意圖Fig．6 Rotor unit schematic

圖7 轉(zhuǎn)子位移響應(yīng)時(shí)間歷程及頻譜Fig．7 Rotor displacement time history and its spectrum

圖8 位移頻譜Fig．8 Displacement spectrum

圖9 位移頻譜Fig．9 Displacement spectrum

三轉(zhuǎn)子臨界轉(zhuǎn)速為1 867 r／min、轉(zhuǎn)速為5 700 r／min時(shí)轉(zhuǎn)子振動(dòng)較大，軸與軸套發(fā)生碰摩。位移傳感器測(cè)得碰摩信號(hào)及頻譜見(jiàn)圖10，除基頻95．8 Hz外，亦產(chǎn)生2倍頻及若干高次倍頻。

圖10 轉(zhuǎn)子位移響應(yīng)時(shí)間歷程及其頻譜Fig．10 Rotor displacement time history and its spectrum

利用雙樹(shù)復(fù)小波進(jìn)行降噪處理，選構(gòu)造分段光滑閾值函數(shù)對(duì)信號(hào)重構(gòu)所得頻譜見(jiàn)圖11，已提取信號(hào)的高頻細(xì)節(jié)信息，其中第二層含669．5 Hz，765．3 Hz，860 Hz，分別為7、8、9倍頻；第三層含382 Hz，為4倍頻；第四層含191．2 Hz，為2倍頻。各頻帶間無(wú)混疊，有用信息獲得充分顯示。

圖11 位移頻譜Fig．11 Disp lacement spectrum

5 結(jié) 論

（1）構(gòu)造的分段圓弧平滑閾值函數(shù)在整個(gè)區(qū)間一階導(dǎo)數(shù)連續(xù)，能避免重構(gòu)信號(hào)振蕩，可有效消除非平穩(wěn)信號(hào)中噪聲。

（2）利用雙樹(shù)復(fù)小波變換分析碰摩非線性動(dòng)態(tài)響應(yīng)信號(hào)提取的轉(zhuǎn)子動(dòng)態(tài)響應(yīng)高頻信號(hào)與小波及小波包變換相比，雙樹(shù)復(fù)小波變換能抑制模式混疊，并能提取更多細(xì)節(jié)信息。

（3）由兩平行但不同的離散小波組成的雙樹(shù)復(fù)小波變換其雙樹(shù)分解系數(shù)信息互補(bǔ)，可減少由下采樣造成的信息丟失，能完整提取碰摩故障信號(hào)特征；其雙樹(shù)實(shí)小波解析特性可抗頻帶混疊并具有平移不變特性，減少信號(hào)信息特征丟失，在故碰摩故障信號(hào)特征提取中效果較好。

［1］Han Q，Zhang Z，Wen B．Periodic motions of a dual-disc rotor system with rub-impact at fixed limiter［J］．Proceedings of the Institution of Mechanical Engineers，Part C：Journal of Mechanical Engineering Science，2008，222（10）：1935－1946．

［2］Chu F，Zhang Z．Bifurcation and chaos in a rub-impact Jeffcott rotor system［J］．Journal of Sound and Vibration，1998，210（1）：1－18．

［3］Lin J，Qu L．Feature extraction based on Morletwavelet and its application for mechanical fault diagnosis［J］．Journal of Sound and Vibration，2000，234（1）：135－148．

［4］Peng Z，He Y，Lu Q，et al．Feature extraction of the rubimpact rotor system by means of wavelet analysis［J］．Journal of Sound and Vibration，2003，259（4）：1000－1010．

［5］Kingsbury N．Complex wavelets for shift invariant analysis and filtering of signals［J］．Applied and Computational Harmonic Analysis，2001，10（3）：234－253．

［6］Selesnick IW，Baraniuk R G，Kingsbury N C．The dualtree complex wavelet transform［J］．Signal Processing Magazine，2005，22（6）：123－151．

［7］Selesnick IW，Baraniuk R G，Kingsbury N C．The dual-tree complex wavelet transform［J］．Signal Processing Magazine，IEEE，2005，22（6）：123－151．

［8］Chang SG，Yu B，VetterliM．Adaptivewavelet thresholding for image denoising and compression［J］．Image Processing，IEEE Transactions on，2000，9（9）：1532－1546．

［9］Donoho D L．De-noising by soft-thresholding［J］．Information Theory，IEEE Transactions on，1995，41（3）：613－627．

［10］Donoho D L，Johnstone I M．Adapting to unknown smoothness via wavelet shrinkage［J］．Journal of the American Statistical Association，1995，432（90）：1200－1224．

［11］吳定海，張培林，任國(guó)全，等．基于雙樹(shù)復(fù)小波包的發(fā)動(dòng)機(jī)振動(dòng)信號(hào)特征提取研究［J］．振動(dòng)與沖擊，2010，29（4）：160－163．

WU Ding-hai，ZHANG Pei-lin，REN Guo-quan，et al．Feature extraction of an engine vibration signal based on dualtree wavelet package transformation［J］．Journal of Vibration and Shock，2010，29（4）：160－163．

［12］Stein C M．Estimation of themean of a multivariate normal distribution［J］．The Annals of Statistics，1981，9（6）：1135－1151．

Rotor-stator rub-impact response analysis based on DTCWT

LI Ji-yong1，LIShun-ming1，CHEN Xiao-hong2，WANG Xing-xing1
（1．College of Energy and Power Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China；2．College of Science，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China））

The de-noising of rub-impact response signals of a rotor-stator system and its characteristics extraction are key problems of rub fault diagnosis．A piecewise smooth threshold function was constructed forwavelet de-noising，and non-stationary signals were de-noised separately by using Wavelet，Wavelet packets and dual-tree complex wavelet transform（DTCWT）．SNR（signal to noise ratio）was taken as an index to evaluate the de-noising effect．The results reveal that the constructed threshold function is highly active in de-noising．As an example，the dynamic response of a single rotor under rub-impactwas calculated and decomposed by DTCWT．The simulation and experiment results show that the dynamic response can be exactly derived due to the complementary property between approximate analyticity and sampling interpolation of DTCWT．It is also proved that the proposedmethod is effective．

rotor-stator rub；de-noising；dynamic response；threshold function；DTCWT

V233．1

：A

10．13465／j．cnki．jvs．2014．22．019

航空自然科學(xué)基金（2012ZD52054）；國(guó)家自然科學(xué)基金項(xiàng)目（61403193）；南航基本科研業(yè)務(wù)費(fèi)科研項(xiàng)目（NS2014081）

2013－07－24 修改稿收到日期：2013－11－21

李紀(jì)永男，博士生，1984年生

李舜酩男，教授，1962年生郵箱：smli＠nuaa．com