摘 要：針對風(fēng)電機組狀態(tài)監(jiān)測中主軸軸承故障特征受到其他傳動部件振動諧波干擾且原始振動信號非平穩(wěn)難以準(zhǔn)確提取故障沖擊特征頻率的問題，提出一種基于諧波干擾消除和譜峭度相結(jié)合的分析方法（HR-Fast-Kurtogram）。利用離線角域重采樣將非平穩(wěn)時域信號轉(zhuǎn)換成平穩(wěn)的角度-階次信號，在階次域通過梳狀陷波器其去除齒箱嚙合成分等干擾階次，應(yīng)用基于短時傅里葉變換（STFT）的快速譜峭度圖提取峭度峰值所在故障子頻帶從而提取軸承故障特征。既可避免非周期性沖擊的干擾，亦可區(qū)分沖擊性較小的非軸承幅值調(diào)制成分和軸承故障沖擊成分。在實際機組軸承故障數(shù)據(jù)應(yīng)用中的對比分析結(jié)果表明：新方法可提高在秒級短時間尺度提取軸承故障特征和識別故障類型的準(zhǔn)確性。

關(guān)鍵詞：風(fēng)電機組；軸承；狀態(tài)監(jiān)測；陷波器；譜峭度

中圖分類號：TH165+.3" " " " " " " " " " 文獻標(biāo)志碼：A

0 引 言

風(fēng)電機組主軸負(fù)責(zé)傳遞葉輪到齒輪箱之間的扭矩，三點布置式雙饋機組主軸承同時承受徑向力和軸向力，是傳動鏈的關(guān)鍵承載部件，其主要失效形式是滾道和滾動體損傷［1］，在振動上表現(xiàn)為規(guī)律性沖擊特征。由于機組運行中轉(zhuǎn)速的實時波動，導(dǎo)致振動的非平穩(wěn)性以及沖擊的非周期性，但所產(chǎn)生的沖擊共振頻帶則相對固定。因此風(fēng)電機組軸承故障診斷的關(guān)鍵在于在非平穩(wěn)信號中尋找對軸承故障敏感的沖擊共振頻帶，并從中提取表征故障信息周期性參數(shù)。

針對這一特點，近年來出現(xiàn)諸多風(fēng)電軸承振動分析方法，固有時間尺度分解（intrinsic time decomposition，ITD）將復(fù)雜的非平穩(wěn)、非線性信號分解成若干個固有旋轉(zhuǎn)分量，可分析調(diào)幅調(diào)頻成分的頻譜特征［2］，但該方法與EMD分解同屬于無監(jiān)督經(jīng)驗信號分解，其不可預(yù)測性對于軸承狀態(tài)監(jiān)測適用性不強。文獻［3］利用盲源分離和流形學(xué)習(xí)分別提取高維和低維特征實現(xiàn)軸承故障的提取和區(qū)分，然而盲源分離對振動信號中獨立源數(shù)目的估計準(zhǔn)確性不高。文獻［4］通過構(gòu)造并行卷積神經(jīng)網(wǎng)絡(luò)實現(xiàn)小波時頻故障特征的增強。之后相繼出現(xiàn)將經(jīng)驗?zāi)B(tài)分解（empirical mode decomposition，EMD）與峭度-相關(guān)系數(shù)準(zhǔn)則［5］、改進的卷積神經(jīng)網(wǎng)絡(luò)［6］、極限學(xué)習(xí)機［7］等訓(xùn)練模型相結(jié)合的方法，獲得了較好的診斷效果，但該類算法對于數(shù)據(jù)量的要求較高，對機組故障診斷的實時性不好。

譜峭度法（Spectral Kurtosis）結(jié)合包絡(luò)譜分析是一種利用故障沖擊在分布上的本質(zhì)特征確定信號中共振頻帶并進行包絡(luò)解調(diào)的方法，最早由Dwyer［8］提出。Antoni［9］在此基礎(chǔ)上改進算法的實時性，提出可用于工業(yè)診斷的快速譜峭度法。通過濾波器組對信號分段層級分解搜尋譜峭度峰值形成譜峭度圖（Kurtogram）直接尋找譜峭度圖中的峰值所代表的帶寬（Bw）和中心頻率。實際上，機組傳動鏈時常會受到瞬時風(fēng)載和風(fēng)切變的影響產(chǎn)生非周期性沖擊，由此衍生出改進由于非故障沖擊造成峭度異常偏大的方法［10-11］。但狀態(tài)監(jiān)測系統(tǒng)中信號長度的限制使譜峭度平均法［11］區(qū)分低頻非周期沖擊特征的效果不佳，并且傳動鏈中的旋轉(zhuǎn)部件除了主軸承還有齒輪箱，齒輪箱的嚙合頻率成分受到行星輪通過效應(yīng)和載荷變化的影響會產(chǎn)生幅值調(diào)制現(xiàn)象，干擾傳統(tǒng)譜峭度法的分析。

針對以上不足，本文提出通過離線角域重采樣和梳狀陷波器去除高次諧波，結(jié)合基于短時傅里葉變換（short-time Fourier transform，STFT）的快速譜峭度分析軸承故障沖擊的方法，實現(xiàn)主軸軸承的故障特征提取和診斷。

1 譜峭度分析

1.1 譜峭度指標(biāo)

峭度是隨機變量的歸一化四階中心矩，衡量信號偏離正態(tài)分布的程度，也是在機械振動信號時域處理中用于檢測信號中沖擊特征的無量綱指標(biāo)，其定義為：

[K=EX-μ4EX-μ22-3] （1）

式中：X——隨機變量；μ——變量均值。然而時域峭度指標(biāo)作為全局性指標(biāo)，不能準(zhǔn)確反映噪聲環(huán)境下的微弱故障沖擊［11］。

1.2 基于STFT的譜峭度

譜峭度法分為對特定頻帶時域復(fù)包絡(luò)的峭度和基于STFT的譜峭度，前者是在對信號子頻帶分解后直接求其復(fù)包絡(luò)的峭度值，后者表示子頻帶每根時頻譜線的峭度值。若非平穩(wěn)信號[xt]可表示為：

[xt=-∞+∞ej2πftGt， fdxf] （2）

式中：t——時間；f——頻率；[Gt， f]——時變沖擊響應(yīng)[gt，s]的傅里葉變換，s為時間延遲。

則離散形式的STFT變換為：

[XwmP，f=n=-∞∞Xnwn-kPe-j2πnf] （3）

式中：w——窗函數(shù)；m——自然數(shù)，m=1，2，3，…；P——時間步長；[Xn]——[xt]的離散采樣。

定義STFT變換后信號的二階經(jīng)驗譜矩為：

[S2nXf=XwkP， f2nk] （4）

式中：[·]——時間平均算子；[·]——取模算子。

基于STFT的譜峭度定義［12］為：

[KYf=S4YfS22Yf-2] （5）

式中：Y——隨機過程。

由于復(fù)數(shù)效應(yīng)，式（1）中的常數(shù)3改為2［9］。

1.3 快速譜峭度算法

快速譜峭度［9］提供了一種對信號快速分層的算法，利用1/3-二叉樹濾波器組實現(xiàn)信號的層級分解。設(shè)[hn]為歸一化截止頻率為[fc=18+ε]的低通濾波器，其中偏移量[ε≥0]。每一組濾波器分別由如式（6）和式（7）所示的低通和高通基礎(chǔ)濾波器構(gòu)成：

[h0n=hnejπn/4] （6）

[h1n=hnej3πn/4] （7）

式中：[h0n]——由[hn]頻移1/8的準(zhǔn)解析濾波器；[h1n]——由[hn]頻移3/8的準(zhǔn)解析濾波器。

若令分解層級為k，以此類推則每一層有2k個分解子頻帶。分別把原始信號分解為式（8）和式（9）所示的子信號：

[c2i（k+1）n=h0j*cik2n] （8）

[c（2i+1）（k+1）n=h1j*cik2n] （9）

式中：[cikn]——原始信號經(jīng)第k層濾波器分解后的信號；[i=0，1，2，…，2k-1]；*——卷積符號；當(dāng)k=0時，[c0n]即為原始信號[xt]。最終譜峭度圖上每一個子頻帶的峭度為：

[Kik=cikn4cikn22-2] （10）

進一步地，為提高頻率分辨率，在每一層之間再插入3×2k個濾波器，形成1/3-二叉樹濾波器組，經(jīng)其濾波后的子信號組成結(jié)構(gòu)前3層如圖1所示。

2 諧波消除

2.1 角域重采樣

風(fēng)電機組運行過程中受到風(fēng)速變化的影響，轉(zhuǎn)速往往處于波動狀態(tài)，使原本周期性的故障沖擊變?yōu)榉侵芷谛詻_擊，并且振動信號中由齒輪箱傳遞到軸承上的變轉(zhuǎn)速諧波信號在頻譜上會產(chǎn)生嚴(yán)重的“頻率模糊”，影響譜峭度分析的故障子頻帶提取結(jié)果。為避免變轉(zhuǎn)速諧波對軸承故障診斷的干擾，利用主控的發(fā)電機轉(zhuǎn)速信號將時域振動信號重采樣到角域，從而將諧波成分平穩(wěn)化。令主控每分鐘轉(zhuǎn)速信號為N，則秒級角速度信號[ω]為：

[ω=N60×2π] （11）

時間-角度關(guān)系可通過角速度[ω]積分得到：

[θt=0tωdt] （12）

根據(jù)計算精度設(shè)定每圈的等角度重采樣間隔q，則離散角度序列[θi]如式（14）所示：

[θi=iq， i=0，1，2，…，L] （13）

其中L的取值為：

[L=floorθtmax/q] （14）

式中：[floor]——向下取整函數(shù)；[tmax]——最大采樣時間。

利用角度序列[θi]對原角度-時間信號進行插值得到新時間信號[tθi]：

[tθi=tiq] （15）

最后利用新時間信號對原始時域振動信號[xt]進行插值獲得角度域信號[xθi]。

2.2 梳狀陷波

梳狀濾波器是一種由等頻率間隔排列的通帶或阻帶組成，只通過或阻止特定中心頻率及其倍頻信號的數(shù)字濾波器。由于機組振動信號中的齒輪箱嚙合頻率和旋轉(zhuǎn)軸轉(zhuǎn)頻總是以基頻及其倍頻關(guān)系出現(xiàn)，因此梳狀陷波器可較好去除倍頻諧波成分。設(shè)計以基頻階次為中心頻率，-3 dB帶寬為0.05的16階IIR梳狀陷波器Hc，若陷波中心頻率分別為0.20和0.25，其幅頻響應(yīng)特性如圖2所示。

從幅頻響應(yīng)特性可看出，梳狀陷波器可濾除任意特定中心頻率和帶寬的基頻以及高次諧波，且阻帶非常狹窄，符合對于消除振動信號中干擾諧波的需要。

3 主軸承故障診斷

3.1 基本步驟

實際應(yīng)用中針對機組主軸承的狀態(tài)監(jiān)測和故障診斷，本方法的基本步驟如下：

1）根據(jù)主控同步轉(zhuǎn)速數(shù)據(jù)對采集的振動數(shù)據(jù)段進行離線角域重采樣；

2）設(shè)計梳狀陷波器在階次域去除重采樣信號中齒輪箱各級嚙合頻率和高速軸基頻及其高次諧波；

3）對干擾諧波消除后的信號應(yīng)用基于STFT的譜峭度分析，提取最大譜峭度所在子頻帶；

4）對提取的子頻帶進行包絡(luò)譜分析，診斷軸承故障類型。

3.2 軸承故障數(shù)據(jù)

以某型號雙饋機組狀態(tài)監(jiān)測系統(tǒng)所采集的額定轉(zhuǎn)速1750 r/min下分別含有主軸承不同部位損傷的兩段軸向振動數(shù)據(jù)為例，原始信號的采樣頻率2560 Hz，采樣時間50 s，信號1不僅含有周期性故障沖擊，還含有偶發(fā)性非故障沖擊，信號2中的故障沖擊更為微弱，淹沒在噪聲和諧波中，案例較為典型。其時域曲線和STFT時頻譜如圖3和圖4所示。

3.3 重采樣結(jié)果

取重采樣頻率[fres]=1/q為80次/轉(zhuǎn)，并使用樣條插值方法進行插值重采樣。重采樣后的STFT時頻譜結(jié)果如圖5所示。

由圖5得知，重采樣后諧波階次譜線不再波動，周期性沖擊變?yōu)榈乳g隔，消除了頻譜模糊，有利于后續(xù)階次消除和譜峭度處理。

3.4 梳狀陷波結(jié)果

根據(jù)齒輪箱各級齒輪的傳動方式和齒數(shù)得出，傳動鏈中含有齒輪箱Ⅰ、Ⅱ級行星齒輪傳動嚙合頻率、高速平行級齒輪傳動嚙合頻率及其高次諧波成分。為濾除主軸承之外的干擾成分，對每個干擾成分都分別應(yīng)用梳狀陷波器。在機組額定轉(zhuǎn)速下的干擾諧波基頻、階次及其相應(yīng)的陷波器參數(shù)如表1所示。

表1中陷波器中心頻率和帶寬參數(shù)均為歸一化頻率。經(jīng)梳狀陷波器處理后的振動信號STFT時頻譜如圖6所示。

消除齒輪箱等傳動鏈其他旋轉(zhuǎn)部件的諧波激勵成分后，軸承故障沖擊更為突出，且陷波帶寬足夠窄，幾乎不會影響故障沖擊共振頻帶本身的質(zhì)量。

3.5 譜峭度分析方法結(jié)果對比

設(shè)定譜峭度圖的最大分解層數(shù)為10層，對原始振動數(shù)據(jù)分別應(yīng)用時域快速譜峭度、基于STFT的快速譜峭度兩種經(jīng)典方法以及本文所提出的基于諧波消除的STFT快速譜峭度（HR-Fast-Kurtogram）方法，以比較各方法的故障提取能力。

3.5.1 時域快速譜峭度診斷結(jié)果

原始信號1時域快速譜峭度圖計算出的最大譜峭度Kmax為29.2，所在的中心頻率[fc]為1120.00 Hz，分解層數(shù)k為3.6，原始信號2時域快速譜峭度圖計算出的最大譜峭度為29.2，所在的中心頻率[fc]為1209.38 Hz，分解層數(shù)k為10，如圖7所示。根據(jù)最大譜峭度子頻帶參數(shù)提取的濾波信號1并非規(guī)律性故障沖擊，如圖8a所示。子頻帶信號的平方包絡(luò)譜顯示，無明顯突出譜線，提取的成分是高頻偶發(fā)性非故障沖擊，如圖8b所示。濾波信號2同樣不存在主軸承故障沖擊和包絡(luò)譜特征譜線，根據(jù)表1可知該信號為齒輪箱高速平行級嚙合2倍頻諧波，諧波的幅值調(diào)制導(dǎo)致計算子頻帶的譜峭度偏大，干擾譜峭度分析結(jié)果，時域快速譜峭度未能提取出故障沖擊，如圖8c和圖8d所示。

3.5.2 STFT快速譜峭度診斷結(jié)果

圖9顯示，原始信號1和2基于STFT的快速譜峭度圖計算出的最大譜峭度值Kmax為5.3和1.5，所在的中心頻率[fc]為613.13和615.00 Hz，分解層數(shù)k均為10。

根據(jù)表1和圖10可知，原始信號1和2所提取的子頻帶信號均是齒輪箱高速平行級齒輪嚙合頻率基頻成分，不存在規(guī)律性故障沖擊。子頻帶信號的平方包絡(luò)譜顯示，不存在軸承故障特征譜線，經(jīng)典的STFT快速譜峭度亦未能提取出軸承故障沖擊成分。

3.5.3 HR-Fast-Kurtogram診斷結(jié)果

原始信號1和2中的干擾諧波被消除后應(yīng)用基于STFT的快速譜峭度圖算法所計算的故障頻帶最大譜峭度值Kmax分別為1.7和0.4，中心階次[oc]均為1.875，分解層數(shù)k均為7，如圖11所示。

其提取出的子頻帶信號為周期性沖擊成分，如圖12a和圖12c所示。根據(jù)圖12b和圖12d中子頻帶的平方包絡(luò)譜，該故障成分的特征階次為0.075和0.087，轉(zhuǎn)換成機組額定轉(zhuǎn)速下的特征頻率為2.19和2.53 Hz，對照表2所示額定轉(zhuǎn)速時該型號軸承的故障特征頻率得知，兩個被測主軸承分別存在外圈和內(nèi)圈滾道損傷。從沖擊峰值和譜峭度值上看，內(nèi)圈故障沖擊相比外圈更加微弱，但沖擊共振頻帶均相同。

現(xiàn)場打開兩個主軸承端蓋排查后發(fā)現(xiàn)信號1對應(yīng)的主軸承外圈滾道表面存在一處材料剝落損傷，信號2對應(yīng)的主軸承內(nèi)滾道表面存在2～3 mm寬的輕微壓痕，如圖13a和圖13b所示。排查結(jié)果驗證了本方法對主軸承不同部位、不同失效形式故障特征提取和診斷的有效性。

4 結(jié) 論

針對風(fēng)電機組主軸承振動信號易受偶發(fā)性非故障沖擊和齒輪箱諧波激勵干擾的問題，提出基于諧波消除和STFT快速譜峭度的故障診斷方法（HR-Fast-Kurtogram）。

在含有偶發(fā)性沖擊和轉(zhuǎn)速波動的故障振動數(shù)據(jù)的應(yīng)用中表明：

1）利用角域重采樣和梳狀陷波剔除齒輪箱的諧波干擾階次可有效避免尋找譜峭度峰值時故障子頻帶的丟失。

2）能提取風(fēng)電機組主軸承內(nèi)圈和外圈失效產(chǎn)生的沖擊時域信號及其包絡(luò)譜特征頻率2.53和2.19 Hz。

3）在敏感性上，對于早期滾道表面寬度2～3 mm的輕微壓痕產(chǎn)生的沖擊依然有效。經(jīng)過與經(jīng)典方法診斷結(jié)果的對比，其故障提取能力和準(zhǔn)確性優(yōu)于時域譜峭度法和單獨的STFT譜峭度法。

[參考文獻]

［1］ 蘇連成， 李興林， 李小俚， 等. 風(fēng)電機組軸承的狀態(tài)監(jiān)測和故障診斷與運行維護［J］. 軸承， 2012（1）： 47-53.

SU L C， LI X L， LI X L， et al. Condition monitoring and fault diagnosis and operating maintenance systems for wind turbine bearings［J］. Bearing， 2012（1）： 47-53.

［2］ 安學(xué)利， 蔣東翔， 劉超， 等. 基于固有時間尺度分解的風(fēng)電機組軸承故障特征提取［J］. 電力系統(tǒng)自動化， 2012， 36（5）： 41-44， 102.

AN X L， JIANG D X， LIU C， et al. Bearing fault feature extraction of wind turbine based on intrinsic time-scale decomposition［J］. Automation of electric power systems， 2012， 36（5）： 41-44， 102.

［3］ 趙洪山， 李浪， 王穎. 一種基于盲源分離和流形學(xué)習(xí)的風(fēng)電機組軸承故障特征提取方法［J］. 太陽能學(xué)報， 2016， 37（2）： 269-275.

ZHAO H S， LI L， WANG Y. Fault feature extraction method of wind turbine bearing based on blind source separation and manifold learning［J］. Acta energiae solaris sinica， 2016， 37（2）： 269-275.

［4］ 孟良， 蘇元浩， 許同樂， 等. 并行卷積神經(jīng)網(wǎng)絡(luò)的風(fēng)電機組故障診斷方法［J］. 太陽能學(xué)報， 2023， 44（5）： 449-456.

MENG L， SU Y H， XU T L， et al. Wind turbine fault diagnosis method based on parallel convolutional neural network［J］. Acta energiae solaris sinica， 2023， 44（5）： 449-456.

［5］ 彭進， 王維慶， 王海云， 等. 基于EEMD峭度-相關(guān)系數(shù)準(zhǔn)則的多特征量風(fēng)電機組軸承故障診斷［J］. 可再生能源， 2016， 34（10）： 1481-1490.

PENG J， WANG W Q， WANG H Y， et al. Fault diagnosis method of wind turbine’s bearing based on EEMD kurtosis-correlation coefficients criterion and multiple features［J］. Renewable energy resources， 2016， 34（10）： 1481-1490.

［6］ 劉洋， 程強， 史曜煒， 等. 基于注意力模塊及1D-CNN的滾動軸承故障診斷［J］. 太陽能學(xué)報， 2022， 43（3）： 462-468.

LIU Y， CHENG Q， SHI Y W， et al. Fault diagnosis of rolling bearings based on attention module and 1D-CNN［J］. Acta energiae solaris sinica， 2022， 43（3）： 462-468.

［7］ 胡超， 沈?qū)殗?謝中敏. 基于EMD-FastICA與DGA-ELM網(wǎng)絡(luò)的軸承故障診斷方法［J］. 太陽能學(xué)報， 2021， 42（10）： 208-219.

HU C， SHEN B G， XIE Z M. Fault diagnosis method of bearing based on EMD-FastICA and DGA-ELM network［J］. Acta energiae solaris sinica， 2021， 42（10）： 208-219.

［8］ DWYER R. Detection of non-Gaussian signals by frequency domain Kurtosis estimation［C］//ICASSP’83. IEEE International Conference on Acoustics， Speech， and Signal Processing. Boston， MASS， USA， 1983.

［9］ ANTONI J. Fast computation of the kurtogram for the detection of transient faults［J］. Mechanical systems and signal processing， 2007， 21（1）： 108-124.

［10］ 張龍， 毛志德， 楊世錫， 等. 基于包絡(luò)譜帶通峭度的改進譜峭度方法及在軸承診斷中的應(yīng)用［J］. 振動與沖擊， 2018， 37（23）： 171-179.

ZHANG L， MAO Z D， YANG S X， et al. An improved Kurtogram based on band-pass envelope spectral Kurtosis with its application in bearing fault diagnosis［J］. Journal of vibration and shock， 2018， 37（23）： 171-179.

［11］ 代士超， 郭瑜， 伍星， 等. 基于子頻帶譜峭度平均的快速譜峭度圖算法改進［J］. 振動與沖擊， 2015， 34（7）： 98-102， 108.

DAI S C， GUO Y， WU X， et al. Improvement on fast kurtogram algorithm based on sub-frequency-band spectral kurtosis average［J］. Journal of vibration and shock， 2015， 34（7）： 98-102， 108.

［12］ ANTONI J. The spectral kurtosis： a useful tool for characterising" " non-stationary" " signals［J］." " Mechanical systems and signal processing， 2006， 20（2）： 282-307.

FAULT DIAGNOSIS OF WIND TURBINE MAIN BEARINGS BASED ON HARMONIC REMOVAL AND SPECTRAL KURTOSIS

Chen Qiang，Qian Xun

（Guodian United Power Technology Co.， Ltd.， Beijing 100039， China）

Abstract：Aiming at the problems that the fault features of the main bearings are interfered by the harmonics of other transmission components in wind turbines’ condition monitoring and the raw vibration signal is non-stationary， making it difficult to accurately extract the characteristic frequency of fault impacts， a combined method of harmonic removal and spectral kurtosis analysis （HR-Fast-Kurtogram） is proposed. Offline angle domain resampling is used to convert non-stationary time domain signals into stationary angle-order signals. Comb notch filters are introduced in the order domain to remove gear meshing frequencies， and Fast-Kurtogram based on short-time Fourier transfer （STFT） is applied to extract the fault sub-bands where the kurtosis peak is located in order to extract bearing fault features. The combined method is able to avoid the interference of non-periodic impacts and distinguish between non-bearing amplitude modulation components with lower impact and bearing fault impact components. The comparative analysis results in the application of real main bearing fault data shows that the proposed method can improve the accuracy of fault features extraction and fault type identification in a short time scale of seconds.

Keywords：wind turbines； bearings； condition monitoring； notch filters； spectral kurtosis