吳 磊 王家序 張 新 劉治汶

1.西南交通大學(xué)機(jī)械工程學(xué)院，成都，6100312.電子科技大學(xué)自動(dòng)化工程學(xué)院，成都，611731

0 引言

風(fēng)電機(jī)組長(zhǎng)期在陣風(fēng)等復(fù)雜交變載荷作用下工作，惡劣工況以及復(fù)雜機(jī)械結(jié)構(gòu)嚴(yán)重影響機(jī)組運(yùn)行安全與維護(hù)保障[1-2]。風(fēng)電機(jī)組關(guān)鍵部件——齒輪和軸承容易出現(xiàn)各種形式的故障，一旦出現(xiàn)故障可能造成巨大損失，因此研究風(fēng)電機(jī)組的運(yùn)行狀態(tài)監(jiān)測(cè)與故障診斷具有重要意義[3-6]。

得益于振動(dòng)信號(hào)中含有豐富的機(jī)械運(yùn)行狀態(tài)信息，振動(dòng)監(jiān)測(cè)成為旋轉(zhuǎn)機(jī)械狀態(tài)監(jiān)測(cè)與故障診斷的有效技術(shù)[7-9]。然而，測(cè)量信號(hào)中由齒輪或軸承局部缺陷所引起的故障沖擊序列經(jīng)常被環(huán)境噪聲以及其他干擾成分所掩蓋，加上傳遞路徑的傳播效應(yīng)，導(dǎo)致齒輪和軸承的故障診斷存在較大困難[10-12]。盲解卷積方法通過(guò)設(shè)計(jì)一個(gè)有限長(zhǎng)脈沖響應(yīng)(finite impulse response，F(xiàn)IR)濾波器，能實(shí)現(xiàn)在未知振動(dòng)源、傳輸通道特性以及噪聲強(qiáng)度的情況下從測(cè)量信號(hào)中恢復(fù)故障沖擊特征，在旋轉(zhuǎn)機(jī)械故障診斷領(lǐng)域引起了廣泛的關(guān)注[13]。在設(shè)計(jì)FIR濾波器時(shí)，通常是基于信號(hào)的某種數(shù)字特征，且待恢復(fù)的信號(hào)成分與干擾成分在此特征上表現(xiàn)出可區(qū)分性[14-15]。如經(jīng)典的MED(minimum entropy deconvolution)方法利用地震波與環(huán)境噪聲在峭度方面的不同，成功獲取地震相關(guān)信息[16]。MED具有調(diào)參少、收斂快等優(yōu)勢(shì)，但用于旋轉(zhuǎn)機(jī)械振動(dòng)信號(hào)分析時(shí)卻傾向于恢復(fù)單個(gè)主導(dǎo)沖擊而非故障沖擊序列，原因在于相同信號(hào)長(zhǎng)度下單個(gè)沖擊的峭度遠(yuǎn)大于沖擊序列的峭度。在機(jī)械故障診斷領(lǐng)域，MED仍得到了部分應(yīng)用。ENDO等[17]結(jié)合MED與自回歸(auto-regressive，AR)，提出AR-MED濾波方法并用于齒輪剝落和齒角裂紋檢測(cè)；YANG等[18]利用MED對(duì)原始信號(hào)進(jìn)行預(yù)處理，提取被強(qiáng)噪聲干擾的軸承特征信息；LIU等[19]利用MED凸顯故障沖擊特征，使得包絡(luò)分析更準(zhǔn)確有效。然而，值得注意的是，在上述應(yīng)用中，MED傾向于恢復(fù)單個(gè)沖擊的缺陷并未得到有效解決。

為克服MED這一缺陷，多種盲解卷積方法被提出并應(yīng)用到機(jī)械故障診斷中，其中最常見的方法包括MCKD(maximum correlated kurtosis deconvolution)[20]、MOMEDA(multipoint optimal minimum entropy deconvolution adjusted)[21]以及CYCBD(maximum second-order cyclostationarity blind deconvolution)[22]。在恢復(fù)故障沖擊序列方面，這三種方法表現(xiàn)出優(yōu)異的性能，但它們卻依賴于待恢復(fù)沖擊序列頻率這一先驗(yàn)知識(shí)，即方法并非全“盲”。實(shí)際應(yīng)用中，受速度波動(dòng)等因素的影響，難以預(yù)先獲知準(zhǔn)確的故障特征頻率，嚴(yán)重影響了三種方法的有效性，從而也限制了方法的適用性[23-25]。近年來(lái)，涌現(xiàn)出不少解卷積方法并被用于機(jī)械故障診斷中。DUAN等[26]提出利用形態(tài)濾波來(lái)提高M(jìn)ED的濾波效果；MA等[27]基于稀疏子空間重編碼提出了一種有效的齒輪箱故障診斷方法；ZHANG等[28]利用連續(xù)振動(dòng)分離這一方法克服調(diào)制效應(yīng)并結(jié)合MED提出了一種齒輪故障診斷方法；MIAO等[29]采用基尼系數(shù)作為盲解卷積指標(biāo)，在旋轉(zhuǎn)機(jī)械故障診斷中表現(xiàn)出較好性能；LIANG等[30]提出的最大平均峭度解卷積方法對(duì)典型干擾具有較好的魯棒性，但不足之處在于所提方法的有效性仍依賴先驗(yàn)故障特征頻率。

為有效恢復(fù)故障沖擊序列、準(zhǔn)確診斷齒輪故障，本文提出了一種新的盲解卷積方法——最大重加權(quán)峭度盲解卷積方法。該方法利用重加權(quán)峭度對(duì)故障信號(hào)中單個(gè)或少量強(qiáng)沖擊干擾具有較好魯棒性，以及不依賴待恢復(fù)故障沖擊序列先驗(yàn)知識(shí)的特性，能有效地解決基于峭度最大化方法傾向于恢復(fù)單個(gè)主導(dǎo)沖擊的問(wèn)題。同時(shí)，相較于常見依賴故障特征頻率先驗(yàn)知識(shí)的非全“盲”方法，所提方法在齒輪故障診斷方面具有更強(qiáng)的適用性。通過(guò)仿真信號(hào)分析以及風(fēng)電機(jī)組故障診斷應(yīng)用研究來(lái)驗(yàn)證所提方法的有效性。

1 常見解卷積指標(biāo)及其主要缺陷

1.1 問(wèn)題描述

在旋轉(zhuǎn)機(jī)械故障診斷中，受多振動(dòng)源以及復(fù)雜傳遞路徑和噪聲干擾等多因素的影響，振動(dòng)信號(hào)中由齒輪局部缺陷所引起的故障沖擊成分往往被其他干擾所掩蓋。由傳感器測(cè)量得到的振動(dòng)信號(hào)x一般可表示為

x=s0*gs+n*gn

(1)

式中，s0為局部缺陷產(chǎn)生的故障沖擊成分；n為干擾成分(如高斯噪聲、電流產(chǎn)生的強(qiáng)沖擊干擾等)；*表示卷積運(yùn)算；gs、gn分別為故障沖擊成分與干擾成分的脈沖響應(yīng)函數(shù)。

盲解卷積旨在設(shè)計(jì)一FIR濾波器h，從振動(dòng)信號(hào)中恢復(fù)得到故障沖擊成分s0，其過(guò)程如圖1所示，可表述為

s=x*h=(s0*gs+n*gn)*h≈s0

(2)

其中，濾波信號(hào)s為s0的估計(jì)值。在此過(guò)程中，利用衡量指標(biāo)(如峭度、相關(guān)峭度、D-范數(shù)等)來(lái)評(píng)估濾波信號(hào)以及作為目標(biāo)函數(shù)求解濾波器是必要的，而不同的指標(biāo)表現(xiàn)出不同的性能，甚至直接決定方法的有效性。

圖1 盲解卷積技術(shù)基本原理Fig.1 The basic principle of blind deconvolution techniques

1.2 峭度

在MED方法中，峭度被用作盲解卷積指標(biāo)。盡管此方法名為“最小熵解卷積”，但它并非使濾波信號(hào)的熵最小，而是通過(guò)最大化濾波信號(hào)的峭度來(lái)求解濾波器h，即

(3)

(4)

式中，Kurt(s)為零均值濾波信號(hào)s的峭度；N為信號(hào)長(zhǎng)度；s[i]為s的第i個(gè)分量。

BUZZONI等[22]證明了在一定條件下，峭度與微分熵成相反關(guān)系。齒輪、軸承在正常運(yùn)行狀態(tài)下，其振動(dòng)信號(hào)近似于高斯白噪聲，信號(hào)峭度接近3，而當(dāng)振動(dòng)信號(hào)中包含由局部缺陷產(chǎn)生的非高斯成分時(shí)，信號(hào)的峭度會(huì)超過(guò)3，可見峭度能在一定程度上反映設(shè)備的健康狀態(tài)。但正如引言中所述，使用峭度作為解卷積指標(biāo)的缺陷在于，方法傾向于恢復(fù)得到單個(gè)主導(dǎo)沖擊而非由局部故障引起的沖擊序列，尤其是當(dāng)測(cè)量信號(hào)中含有因電流等導(dǎo)致的強(qiáng)沖擊干擾時(shí)，這一缺陷會(huì)更加明顯。

1.3 其他常見解卷積指標(biāo)

為解決上述問(wèn)題，多種解卷積指標(biāo)被提出并形成新的方法，最常見的包括相關(guān)峭度[20]、多點(diǎn)D-范數(shù)[21]以及二階循環(huán)平穩(wěn)指標(biāo)[22]等。

基于故障沖擊的周期性，相關(guān)峭度定義為

(5)

其中，Ms為位移數(shù)，Ts為故障信號(hào)周期，當(dāng)Ms=1且Ts=0時(shí)，相關(guān)峭度退化為峭度。MCKD通過(guò)最大化濾波信號(hào)的CK來(lái)求解濾波器，即

(6)

相比峭度，CK能夠更好地表征沖擊序列，因此MCKD在故障診斷中被廣泛提及。

通過(guò)引入目標(biāo)向量t，設(shè)定目標(biāo)沖擊的位置以及權(quán)重，MOMEDA定義了多點(diǎn)D-范數(shù)：

(7)

該指標(biāo)充分考慮了故障沖擊特征的周期性，此外，作為非迭代的求解方法，MOMEDA效率更高。

考慮到旋轉(zhuǎn)機(jī)械局部故障產(chǎn)生的沖擊序列會(huì)表現(xiàn)出循環(huán)平穩(wěn)特性，CYCBD采用二階循環(huán)平穩(wěn)指標(biāo)，定義如下：

(8)

(9)

(10)

其中，k為樣本指數(shù)，L為濾波器長(zhǎng)度。CYCBD在恢復(fù)故障沖擊序列方面表現(xiàn)出非常優(yōu)異的性能。

由上述指標(biāo)計(jì)算公式可見，這三種常見方法都極大程度上依賴于先驗(yàn)故障特征頻率，而實(shí)際應(yīng)用中難以預(yù)先獲知準(zhǔn)確的故障特征頻率，致使方法的適用性受到極大限制。因此，探究不依賴先驗(yàn)知識(shí)的全“盲”方法意義重大。

2 最大重加權(quán)峭度盲解卷積方法

為有效恢復(fù)故障沖擊序列、準(zhǔn)確診斷齒輪故障，本節(jié)介紹一種全“盲”方法——最大重加權(quán)峭度盲解卷積。信號(hào)的重加權(quán)峭度計(jì)算過(guò)程如下：

(1)對(duì)濾波信號(hào)s進(jìn)行M等分，得到各子段信號(hào)sm(m=1，2，…，M)。

(2)計(jì)算各子段信號(hào)的峭度Kurtm。

(3)對(duì)Kurtm進(jìn)行升序排列，并將其表示為一行向量的形式，即Kurtasc。

(4)計(jì)算Kurtm對(duì)其和的權(quán)重，即

(11)

(5)對(duì)Wm進(jìn)行降序排列，同樣將其表示為一行向量的形式，即Wdesc。

(6)利用重排后的權(quán)重向量Wdesc對(duì)重排后的峭度向量Kurtasc進(jìn)行加權(quán)，得到信號(hào)的重加權(quán)峭度RK：

(12)

根據(jù)上述計(jì)算過(guò)程可見，當(dāng)M=1時(shí)，RK退化為峭度。為直觀呈現(xiàn)RK在表征故障沖擊序列方面的優(yōu)勢(shì)，計(jì)算了圖2所示的單個(gè)沖擊與周期性沖擊序列(信號(hào)長(zhǎng)度相同)的RK值。如表1所示，單個(gè)沖擊的峭度遠(yuǎn)大于沖擊序列的峭度，而單個(gè)沖擊的RK遠(yuǎn)小于沖擊序列的RK(M>1)，同時(shí)沖擊序列的峭度幾乎等于其RK。因此，以重加權(quán)峭度作為盲解卷積指標(biāo)，可有效地解決基于峭度最大化方法傾向于恢復(fù)單個(gè)主導(dǎo)沖擊的問(wèn)題。

(a)單個(gè)沖擊

表1 不同M值下單個(gè)沖擊與周期性沖擊序列的RK值

M的取值不依賴具體的公式，理論上不超過(guò)信號(hào)長(zhǎng)度的任意正整數(shù)都能使算法正常運(yùn)行，但每一個(gè)子段信號(hào)長(zhǎng)度不能太小，否則子段信號(hào)的峭度Kurtm將失去物理意義。對(duì)多組仿真信號(hào)和實(shí)際齒輪振動(dòng)信號(hào)分析后發(fā)現(xiàn)，當(dāng)M=4時(shí)，大部分?jǐn)?shù)據(jù)分析均能取得理想結(jié)果，因此在使用所提方法時(shí)推薦M取值為4。

基于上述分析，最大重加權(quán)峭度盲解卷積求解FIR濾波器可表示為如下所示的最大化問(wèn)題：

(13)

為求得RK最大值，令RK對(duì)濾波器系數(shù)的偏導(dǎo)數(shù)等于0，即

(14)

首先，

(15)

其中，Kurtθ為第θ段濾波信號(hào)的峭度，可見，式(14)是高度非線性的，難以直接求解。因此，采用迭代求解法，通過(guò)不斷更新濾波器系數(shù)的方式求得方程的有效近似解。為簡(jiǎn)化運(yùn)算，令?Kurtθ/?h[k]=0得到M個(gè)濾波器，然后對(duì)這M個(gè)濾波器加權(quán)求和即可得到更新后的濾波器，過(guò)程如下：

(16)

根據(jù)式(2)，有

(17)

(18)

則可得到下式：

(19)

進(jìn)而，得到其中一個(gè)濾波器，則

(20)

(21)

其中，xθ、sθ分別對(duì)應(yīng)第θ段子段信號(hào)及其濾波信號(hào)，Xθ為xθ對(duì)應(yīng)的Toeplitz矩陣。從而，得到濾波器更新公式：

(22)

綜上，所提方法整體流程如下：

(1)輸入測(cè)量信號(hào)x，指定參數(shù)M=4，隨機(jī)初始化濾波器h系數(shù)，如h=[0 0 … 1 … 0 0]T。

(2)根據(jù)式(2)計(jì)算濾波信號(hào)s，根據(jù)式(12)計(jì)算RK(s)。

(3)根據(jù)式(22)更新濾波器h系數(shù)。

(4)重復(fù)步驟(2)和步驟(3)使RK(s)達(dá)到最大。

(5)步驟(4)中最大RK(s)對(duì)應(yīng)的濾波信號(hào)即為目標(biāo)濾波信號(hào)。

3 仿真信號(hào)分析

本節(jié)進(jìn)行仿真信號(hào)分析以驗(yàn)證所提方法對(duì)恢復(fù)故障沖擊序列的有效性，同時(shí)將本領(lǐng)域常見方法(MED、MCKD、MOMEDA以及CYCBD)作為對(duì)比研究來(lái)探討所提方法的優(yōu)勢(shì)。

考慮到實(shí)際測(cè)量信號(hào)中包含的環(huán)境噪聲，以及可能含有強(qiáng)沖擊干擾(如電流干擾)和諧波分量成分，仿真信號(hào)x(t)的組成如下：

(23)

式中，s0(t)為故障沖擊序列(頻率為取100 Hz)；b(t)為強(qiáng)沖擊干擾；c(t)為諧波分量；w(t)為信噪比-10 dB的高斯白噪聲；ρ為阻尼系數(shù)。

(a)故障沖擊s0(t)

各方法的濾波信號(hào)如圖4所示，可見，本文方法與MCKD、MOMEDA以及CYCBD均準(zhǔn)確恢復(fù)了故障沖擊序列，而MED只得到單個(gè)主導(dǎo)沖擊。然而，正如1.2節(jié)所述，MCKD、MOMEDA和CYCBD作為非全“盲”方法，需要提前指定待恢復(fù)沖擊序列的頻率，而圖4中的分析結(jié)果正是建立在指定頻率完全等于待恢復(fù)沖擊序列頻率基礎(chǔ)之上的。為進(jìn)一步探究這三種方法對(duì)故障頻率這一先驗(yàn)的依賴特性，圖5給出了它們?cè)谥付l率為103 Hz時(shí)的濾波信號(hào)?？梢?，即便是指定頻率與故障沖擊序列實(shí)際頻率存在較小偏差時(shí)，這三種方法均無(wú)法有效恢復(fù)故障沖擊序列。

(a)本文方法

由上述分析可見，與同為全“盲”的經(jīng)典MED相比，所提方法解決了MED傾向于恢復(fù)單個(gè)主導(dǎo)沖擊的問(wèn)題，在恢復(fù)故障沖擊序列中表現(xiàn)出更好的潛力。而相較于常見非全“盲”的MCKD、MOMEDA以及CYCBD，本文方法在保證分析結(jié)果準(zhǔn)確的基礎(chǔ)上，不需要有關(guān)故障沖擊的先驗(yàn)條件，在實(shí)際應(yīng)用中具有更強(qiáng)的適用性。

(a)MCKD

由于實(shí)際場(chǎng)景中測(cè)量信號(hào)還可能同時(shí)包含多個(gè)沖擊干擾，因此在上述仿真信號(hào)x(t)的基礎(chǔ)上再添加多個(gè)沖擊干擾b′(t)(具體參數(shù)如表2所示)，得到新的仿真信號(hào)x′(t)。b′(t)的表達(dá)式為

(24)

表2 b′(t)的具體參數(shù)

b′(t)和x′(t)波形分別如圖6a、圖6b所示，可見，仿真信號(hào)中的沖擊干擾較為明顯。如圖6c所示，所提方法仍準(zhǔn)確恢復(fù)了故障沖擊序列，驗(yàn)證了本文方法對(duì)分析存在多個(gè)沖擊干擾的信號(hào)的有效性。

(a)多個(gè)沖擊干擾b′(t)

4 風(fēng)電機(jī)組故障診斷應(yīng)用研究

工程應(yīng)用中，由于環(huán)境噪聲以及其他干擾的影響，測(cè)量信號(hào)組成十分復(fù)雜，使得故障診斷具有一定難度。本節(jié)采用風(fēng)電機(jī)組故障診斷實(shí)際案例分析，來(lái)驗(yàn)證本文方法在實(shí)際復(fù)雜工程應(yīng)用中的有效性，同時(shí)與MED進(jìn)行對(duì)比研究。因上述其他方法均依賴故障頻率這一先驗(yàn)條件，而實(shí)際信號(hào)的故障頻率難以準(zhǔn)確獲知，因此本案例中不再將它們作為對(duì)比方法。

本工程數(shù)據(jù)來(lái)自于某風(fēng)電機(jī)組振動(dòng)監(jiān)測(cè)，振動(dòng)監(jiān)測(cè)中，安裝在齒輪箱輸出小齒輪軸上方殼體表面的加速度傳感器(圖7)監(jiān)測(cè)到了異常振動(dòng)。信號(hào)采樣頻率為97 656 Hz，轉(zhuǎn)速計(jì)測(cè)得齒輪箱輸出小齒輪軸轉(zhuǎn)速為1800 r/min。截取信號(hào)51 200點(diǎn)以分析異常振動(dòng)出現(xiàn)的原因，如圖8a所示，由于測(cè)量信號(hào)波形中觀察不到任何故障特征信息，故無(wú)法對(duì)異常振動(dòng)做出有效診斷。

圖7 風(fēng)電機(jī)組傳動(dòng)系統(tǒng)及齒輪箱輸出軸小齒輪的加速度計(jì)測(cè)點(diǎn)位置示意圖Fig.7 Schematic diagram of wind turbine transmission system and accelerator location of the output pinion gear

本文方法與MED的濾波信號(hào)分別如圖8b、圖8c所示，同時(shí)濾波信號(hào)對(duì)應(yīng)的包絡(luò)譜分別如圖9a、圖9b所示。根據(jù)濾波信號(hào)，本文方法恢復(fù)得到的周期性沖擊序列的頻率為29.6 Hz(1/0.0338 Hz)，十分接近輸出小齒輪理論故障特征頻率(30 Hz)。此外，在其濾波信號(hào)的包絡(luò)譜中，可以清楚地觀察到小齒輪故障特征頻率及其多組倍頻。因此，可初步認(rèn)為小齒輪存在故障。在之后的停機(jī)開箱檢測(cè)中發(fā)現(xiàn)輸出小齒輪某一輪齒確已出現(xiàn)損傷(圖10)，驗(yàn)證了本文方法的有效性。相比之下，MED的濾波信號(hào)卻得到的是單個(gè)主導(dǎo)沖擊，如圖8c所示。并且，在其濾波信號(hào)的包絡(luò)譜中只有故障頻率處的譜線比較明顯，無(wú)法觀察到其他倍頻成分，如圖9b所示，因此仍難以對(duì)異常振動(dòng)做出準(zhǔn)確可靠的診斷。

(a)測(cè)量信號(hào)

(a)本文方法濾波信號(hào)的包絡(luò)譜

圖10 齒輪箱輸出小齒輪輪齒損傷Fig.10 The localized damage on the output pinion gear

5 結(jié)論

(1)提出了最大重加權(quán)峭度盲解卷積方法，解決了基于峭度最大化解卷積方法傾向于恢復(fù)單個(gè)主導(dǎo)沖擊的問(wèn)題，并且作為一種全“盲”方法，不需要有關(guān)故障沖擊序列的先驗(yàn)知識(shí)，較非全“盲”方法具有更強(qiáng)的適用性。

(2)仿真信號(hào)分析以及風(fēng)電機(jī)組故障診斷應(yīng)用研究驗(yàn)證了所提方法的有效性，同時(shí)與傳統(tǒng)盲解卷積方法(MED、MCKD、MOMEDA和CYCBD)的對(duì)比結(jié)果進(jìn)一步凸顯了所提方法的優(yōu)勢(shì)。