張龍,黃婧,吳榮真,宋成洋,王朝兵,2
(1.華東交通大學機電與車輛工程學院,南昌330013;2.中車戚墅堰機車有限公司,江蘇常州213011)
減少停機時間成本和實現(xiàn)接近零的停機時間是預診斷的最終目標[1]。然而,在實際發(fā)生故障之前,如果沒有對剩余使用壽命的準確預測,就不可能實現(xiàn)預診斷的所有優(yōu)點。不準確的預測信息可能會導致不必要的維護,例如部件的早期更換等。性能退化評估(Performance degradation assessment,PDA)是實現(xiàn)預診斷的前提與基礎(chǔ),對充分實現(xiàn)預測維修的潛力起著至關(guān)重要的作用。
齒輪作為旋轉(zhuǎn)機械的關(guān)鍵部件之一,其性能的好壞直接決定著設(shè)備工作性能的優(yōu)劣。其一旦發(fā)生故障,將直接影響到機械設(shè)備的正常安全運行甚至造成重大安全事故[2]。因此如何對齒輪實現(xiàn)在役狀態(tài)監(jiān)測和性能退化評估具有重大意義。
對設(shè)備的信號進行合適的狀態(tài)特征提取是故障診斷與預測的前提[3]。振動信號由于具有信息量大、易采集等優(yōu)點而被廣泛采用。常見的基于振動信號的設(shè)備狀態(tài)特征提取可分為時域統(tǒng)計分析(均方根值、峭度、歪度)、頻域分析(幅值譜分析、倒頻譜分析)、時頻域分析(二次型時頻分布、小波分析)以及時序模型分析(自回歸滑動平均模型、AR 模型)等[4-8]。其中,時序模型分析法,尤其是自回歸時序(Autoregressive,AR)模型分析方法,其模型參數(shù)因具有表征系統(tǒng)狀態(tài)的能力且對系統(tǒng)的狀態(tài)變化敏感程度高,而在故障診斷領(lǐng)域中應用廣泛。如文獻[9]利用經(jīng)驗模態(tài)分解對齒輪信號進行預處理后,再AR 模型表征齒輪振動信號系統(tǒng),將其自回歸系數(shù)的混沌特征參數(shù)導入支持向量機中,完成齒輪故障診斷。文獻[10]先利用小波變換技術(shù)對齒輪原始振動信號進行降噪處理,提取降噪后的振動信號的AR 模型系數(shù)作為特征,結(jié)合主成分分析實現(xiàn)了齒輪故障診斷。
性能退化評估的本質(zhì)為將待測樣本信號與無故障基準模型之間進行相似性度量。近年來,一些概率相似度的性能退化評估模型被相繼提出,如隱馬爾科夫模型(Hidden markov model,HMM)、高斯混合模型(Gaussian mixed models,GMM)等?;诟怕氏嗨贫仍u估模型的核心為建立無故障下的密度模型并以此進行異常檢測。如馮輔周等[11]首先提取正常狀態(tài)下滾動軸承振動信號的小波相關(guān)排列熵,并將其作為特征向量構(gòu)建無故障狀態(tài)下的HMM模型進行性能退化識別。Heyns等[12]將GMM和負對數(shù)似然值相結(jié)合進行齒輪狀態(tài)識別,計算每段信號的負對數(shù)似然值作為度量偏離正常信號密度分布的指標,計算信號的角度同步平均來檢測是否出現(xiàn)故障,完成對齒輪箱故障位置和故障程度的判斷。
基于特征提取的概率相似度評估模型期望通過合適的信號處理方法對相應數(shù)據(jù)進行深層次的信息挖掘,以提高特征對故障程度的敏感性、一致性等。但在實際應用中仍存在一些問題:1)需要足夠大的樣本用以訓練;2)計算復雜,GMM 及HMM 等的訓練和測試過程復雜;3)過早飽和現(xiàn)象,當HMM等概率評估模型方法表明待測樣本與無故障基準模型的相似度為零時,存在設(shè)備并未完全進入真正失效狀態(tài)的狀況,即模型極限值早于真實失效值。
基于重構(gòu)的性能退化評估方法通過建立無故障基準模型對待測信號特征向量進行重構(gòu),將重構(gòu)后的差異性來度量設(shè)備性能劣化程度。如Zhan 等[13-14]以自適應自回歸模型為基礎(chǔ),首先建立一個多工況下健康齒輪振動信號殘差最小的折中模型,最后通過殘差信號的3σ準則或K-S測試得到性能退化程度指標,有效的避免了概率相似度評估模型所存在的過早飽和現(xiàn)象。
稀疏表征最早由Hubel和Wiesel[15]在1951年研究貓的視覺感受細胞中提出。其基本思想是在過完備字典中以盡可能少的原子實現(xiàn)該信號的簡潔表達,近年來大都被運用于機械故障診斷領(lǐng)域當中[16-18],而尚未將其運用于機械設(shè)備PDA當中。字典學習作為稀疏分解的一大關(guān)鍵步驟,其本質(zhì)為一種模式識別方法,將其運用于PDA 的核心思路為通過建立系統(tǒng)或設(shè)備健康狀態(tài)下過完備字典,并以此字典實現(xiàn)對待測特征向量的特征重構(gòu),利用待測特征向量于重構(gòu)特征向量之間的差異性來反映系統(tǒng)或設(shè)備的健康狀態(tài)。
綜上所述,本文首先利用AR 模型提取信號的AR 系數(shù)作為其特征向量,以此來表征信號的狀態(tài)特征。將得到的AR 系數(shù)作為特征向量輸入正常運行狀態(tài)下得到的字典模型中進行最優(yōu)重構(gòu)估計,結(jié)合引用的均方根誤差作為故障程度指標實現(xiàn)齒輪的性能劣化評估。
AR 模型作為一種隨機信號參數(shù)化建模的重要方法,其將隨機噪聲定義為白噪激勵線性系統(tǒng)的響應,通過參數(shù)模型對信號進行描述[19]。AR 模型的參數(shù)對系統(tǒng)狀態(tài)變化敏感,且能夠表征系統(tǒng)狀態(tài)特征,故將其作為多變量狀態(tài)估計的輸入?yún)?shù)。
取時間序列x(t),AR 模型的分析階數(shù)為p,則關(guān)于時間序列的p階AR 模型可以表示為
式中: e(t)為 AR 模型的殘差; φi為第i 項的系數(shù)。
AR 模型參數(shù)估計的本質(zhì)為選取合適的AR 系數(shù)使模型的殘差 e(t)為高斯白噪聲。本文用最小二乘法和BIC準則分別計算模型系數(shù)和選擇模型階數(shù)。具體步驟如下:
1)首先確定一個AR 模型合適的分析階數(shù)p,這里取分析階數(shù)p分別取1,2,···,100;
2)通過最小二乘法分別求的各階次下的自回歸參數(shù) φi(i=1,2,···,100),構(gòu)造式(1)所示的AR 模型,進而得到殘差 e(t);
3)分別計算各階次殘差 e(t)的BIC 值,根據(jù)所得BIC值最小確定各階次所對應的最優(yōu)階數(shù)。
信號稀疏表示是信號處理中一個重要的領(lǐng)域,目的是在已知的冗余字典中找到較少的原子線性組合來表示信號。假定一個信號矩陣Y={y1,y2,···,yn,}, 字典 D={d1,d2,···,dk,},且D 為過完備冗余字典,即k?n,根據(jù)稀疏理論,信號Y 可以分解為
式中:α為Y 在字典D 下的稀疏表示系數(shù),α=[α1,α1,···αk]T。
由于 k ?n,表達式欠定,需要增加約束條件得到最優(yōu)解的系數(shù),在稀疏問題中需要得到最少非零值的系數(shù)[20],其表達為
式中‖·‖0為l0范數(shù),代表一個向量中非零值的個數(shù)[21]。
但是在實際求解中 l0范數(shù)是難以優(yōu)化求解的,會導致NP難問題, l1范 數(shù)是 l0范數(shù)的最優(yōu)凸近似,l1范 數(shù)相較 l0范數(shù)更容易優(yōu)化求解,其優(yōu)化式為
式中 ‖·‖0為 l1范數(shù)。
式(4)可以轉(zhuǎn)化為
式(5)的字典求解部分可以寫成[21]
稀疏表示過程可分為字典學習和信號的稀疏編碼兩部分。最優(yōu)方向法(Method of optimal directions,MOD)[22]和KSVD(K Singular value decomposition)[21]等算法通常用于字典學習,基追蹤(Basic pursuit,BP)、匹配追蹤(Matching pursuit,MP)[23]和正交匹配追蹤(Orthogonal matching pursuit,OMP)[24]等算法通常用于求稀疏編碼。在本文中采用KSVD進行字典學習,OMP進行稀疏編碼。由于字典學習和稀疏編碼屬于同一個優(yōu)化求解過程,所以在本文中兩者統(tǒng)稱為字典學習。
KSVD屬于迭代算法的一種,它通過對字典原子稀疏編碼和更新字典原子兩部分來更好地擬合數(shù)據(jù)。更新字典原子時對誤差矩陣進行SVD分解,根據(jù)誤差最小原則,更新后的字典原子和對應的系數(shù)使用誤差最小的分解項代替。
KSVD的實現(xiàn)步驟如下:
1)字典初始化:初始字典選擇訓練樣本集Y ∈Rm×n中隨機的K 個列向量,D(0)∈Rm×K,K<n。設(shè)置迭代指標J=1。
2)稀疏編碼:利用步驟1得到的字典D(J),對樣本集中每一個樣本 Yi通過OMP算法計算系數(shù)向量 αi。
3)字典更新:每次更新一列,逐列更新字典D={d1,d2,···,dk,}。當更新字典第i 列原子 di時,誤差矩陣 Ei的表達式為
式中 αTj為α 中的第j 行。
取出稀疏矩陣第i 行向量 αTj中不為0的集合xTi和對應的位置索引集 ωi:
從 Ei中取出 ωi對 應的列,組合得到 Ei′,對 Ei′進行奇異值分解為Ei′=UΔVT,取U 的第1列來替換第i 列的字典原子 di,取Δ(1,1)與V 的第1列的乘積x′Ti=Δ(1,1)V(:,1)T來 替換 xTi,完成更新。設(shè)置J=J+1。
OMP的實現(xiàn)步驟如下:
步驟1對于給定的正交字典,計算信號 yi和每個字典原子的內(nèi)積,找到內(nèi)積最大時對應的字典原子 dr1,該原子為 yi最匹配的原子,其表達式為
yi第一次分解可被分解為 dr1的垂直投影分量和殘值兩部分,其表達式為
步驟2對殘差 R1同樣進行步驟1的分解,可得到第二次分解。分解次數(shù)根據(jù)稀疏度或精度要求確定。經(jīng)過K 步分解后,可得殘差 RK的表達式為
步驟3 yi信號分解計算式為
式中R0=yi,即可完成OMP分解。
使用字典學習算法建立性能退化評估模型,訓練階段使用無故障樣本的AR 模型系數(shù)作為字典學習的訓練數(shù)據(jù)進行KSVD計算,可得到對應的基準字典Dnormal,該字典只適用于無故障狀態(tài)信號的稀疏表示,其他故障狀態(tài)的信號在Dnormal上不能很好地進行稀疏表示。當測試信號的AR 系數(shù)與無故障樣本的AR 系數(shù)相似程度較高時,在Dnormal上 進行稀疏表示得到的重構(gòu)AR 系數(shù)和殘差與測試信號的AR 系數(shù)和殘差誤差較小,表明該測試信號故障程度較小。同理當測試信號在基準字典得到的重構(gòu)AR 系數(shù)和殘差的重構(gòu)誤差較大時,表明該測試信號與無故障狀態(tài)信號差別較大,此時故障程度也較大。使用均方誤差(Mean squared error,MSE)來檢測重構(gòu)值和原值之間的偏差,并作為故障程度評估的指標DI。MSE 表達式為
基于自回歸模型和字典學習對齒輪故障程度進行評估的流程如圖1所示。
圖1 齒輪性能退化評估流程圖
齒輪性能退化評估主要過程如下:
1)采集齒輪無故障信號構(gòu)建自回歸模型,得到無故障信號的AR 模型系數(shù)和殘差,將AR 系數(shù)進行KSVD計算得到基準字典 Dnormal。
2)計算待測齒輪信號的AR 模型系數(shù)和殘差es,使用OMP在基準字典Dnormal上計算AR 系數(shù)的稀疏系數(shù)并重構(gòu)AR 系數(shù),帶入待測信號的AR 模型中計算重構(gòu)殘差 er。
3)計算殘差 es和重構(gòu)殘差 er的MSE值,作為退化指標DI來評估齒輪故障程度。
實驗數(shù)據(jù)來源于法國CETIM 的齒輪全壽命實驗數(shù)據(jù)[25]。實驗對象為一級減速齒輪箱,圓柱齒輪齒輪箱詳細參數(shù)見表1。
表1 圓柱齒輪齒輪箱參數(shù)
驅(qū)動馬達轉(zhuǎn)速為1000 r/min,齒輪嚙合頻率為330 Hz,信號采樣頻率為20000 Hz,以主動輪為測試齒輪,每天采樣一次,每次采樣時長3 s,一次采樣60000個數(shù)據(jù)點。每天采樣結(jié)束后停機一次觀測齒輪的狀況,觀測結(jié)果如表2所示。在第6 d 測試齒輪的第二齒出現(xiàn)剝落,但直到失效前剝落范圍并沒有擴大,在第8 d 第十六齒出現(xiàn)早期剝落并在第12 d齒輪全面剝落后停機。圖2a)和圖2b)為主動輪第10 d 和第11 d第二齒發(fā)生剝落故障的齒面照片,可以看出此時故障程度相同(圖中黑色圓圈表示齒面剝落面積,面積的尺寸大小可反映故障的嚴重程度),而從圖2c)看出第11 d第十六齒的剝落面積較大。結(jié)合圖3可知,剝落故障只有在剝落面積較大時才能在時域波形圖上變現(xiàn)出明顯的沖擊特征,所以齒輪的失效主要是由第十六齒的較大面積剝落引起。
表2 齒輪箱疲勞試驗每日停機觀測結(jié)果
圖2 齒輪故障程度
圖3齒輪全壽命時域圖
圖3 是齒輪全壽命的時域波形圖,圖3中第1 d到第10 d 的振動信號幅值變化不大,在第11 d 幅值才增大,有明顯的瞬態(tài)沖擊成分,齒輪進入失效狀態(tài)。圖4是齒輪信號在第1 d,第10 d 和第11 d的包絡(luò)譜,第1天到第10 d 齒輪的嚙合頻率及其諧波分量能量較大(fm為嚙合頻率,2fm為諧波分量),可以反映此時齒輪為正常狀態(tài)。第11 d 的包絡(luò)譜以轉(zhuǎn)頻為主,此時齒輪出現(xiàn)局部異常。
圖4 部分信號包絡(luò)譜
首先計算常用的6個時域指標,將齒輪數(shù)據(jù)中每天的數(shù)據(jù)作為一個樣本,共12個樣本,時域指標圖如圖5所示。在圖5中均方根值、峭度值和波形因子在前10個樣本即前10 d 變化不大,在11個樣本之后發(fā)生較大面積剝落時才發(fā)生顯著變化,而且均方根值在前10 d 還呈現(xiàn)下降趨勢,可能的原因是前幾天齒輪處于磨合狀態(tài),振動能量較大。峰值因子、脈沖因子和裕度因子后期曲線與實際觀測結(jié)果較為一致,但前期的故障趨勢并不明顯。時域指標圖表明常用時域指標不能很好地反映齒輪故障的趨勢,對早期故障不敏感。
圖5 齒輪時域指標圖
使用上述所提模型對齒輪數(shù)據(jù)進行性能退化評估。將原始數(shù)據(jù)按時間順序分割為120個樣本,每個樣本包含6000個數(shù)據(jù)點,其中每10個樣本對應一天的采集數(shù)據(jù)。根據(jù)1.2節(jié)所述,在字典學習中為了得到冗余字典需字典原子個數(shù)大于字典維數(shù),而字典原子維數(shù)與信號段維數(shù)相同,這就需要在實際操作中采集更多的信號樣本。針對本案例數(shù)據(jù)信號樣本較少的情況,可將原始數(shù)據(jù)進行擴充來豐富訓練樣本。
訓練階段使用前35個樣本作為測試樣本,為避免相鄰信號段間的人為邊界干擾并保持因劃分產(chǎn)生的沖擊時移的魯棒性[26],使用循環(huán)移位的方式對樣本數(shù)據(jù)進行擴充,將每個基準樣本段進行10次循環(huán)移位后擴充為10個樣本段,擴充后得到350個測試樣本。圖6為循環(huán)移位過程。其中x1為原始信號,x2為循環(huán)移位后的新信號,n為信號數(shù)據(jù)點數(shù),n1為循環(huán)移位的位數(shù),n2為固定不動的位數(shù)。整個過程為:將樣本信號段x1的前n1個數(shù)據(jù)點順序不變移動到信號段的末尾,生成新的信號段x2,完成一次循環(huán)移位。本文中n1=1。
圖6 循環(huán)位移
AR 模型建立階段,采用BIC準則確定最佳模型階數(shù),BIC值變化曲線如圖7所示,可知最佳模型階數(shù)r=80,使用Burg 算法可計算得到訓練樣本的AR 系數(shù)。
圖7 BIC 值變化曲線
使用所提模型進行評估前需要先對超參數(shù)進行優(yōu)化選擇。模型中涉及到的主要參數(shù)有字典原子維數(shù)n,字典原子個數(shù)K,訓練樣本的樣本數(shù)N,稀疏度L,迭代次數(shù)I。由上述已知字典原子維數(shù)n= 80,訓練樣本原子數(shù)N =350,為對剩余的超參數(shù)進行優(yōu)化,本文采用單因素分析法。即分析某一參數(shù)時,將其他參數(shù)固定,來檢驗此參數(shù)對模型的影響。引入均方根誤差(Root mean square error,RMSE)作為評價標準,RMSE 表達式為
對字典原子個數(shù)K 進行分析時,為保證字典是冗余字典,選擇K 的取值范圍為80~340,間隔為10,其他兩個參數(shù)分別隨機設(shè)置為L=14,I =30,使用KSVD和OMP算法進行字典學習重構(gòu),得到不同字典原子個數(shù)對應的均方根誤差變化如圖8所示。由圖8中可知當K =270時重構(gòu)誤差最小,且滿足冗余字典的要求。而隨著字典原子維數(shù)的增大,K 越趨近于訓練樣本原子數(shù)時重構(gòu)誤差變得越不穩(wěn)定,綜合分析后K 取270。
圖8 字典原子個數(shù)K 值不同時的RMSE
分析稀疏度L對重構(gòu)誤差的影響時,字典原子維數(shù)n,訓練樣本原子數(shù)N,字典原子個數(shù)K 均已確定,選擇稀疏度L的取值范圍為1~40,間隔為1。圖9為不同稀疏度L對應的均方根值誤差變化,圖9中可以看出當L=12時重構(gòu)誤差最小,且隨著稀疏度L的增大重構(gòu)誤差趨于穩(wěn)定。而L的增加,會使計算時間顯著增加,故L選擇12較為合適。
圖9 稀疏度L 不同時的RMSE
最后對迭代次數(shù)I 進行分析,設(shè)置I 的取值范圍為1~ 40,間隔為1。不同迭代次數(shù)對應的均方根值誤差變化如圖10所示。圖10中可以看出,I =32時RMSE 最小,選取最優(yōu)參數(shù)I = 32。
圖10 迭代次數(shù)I 不同時的RMSE
綜上分析,字典學習的主要參數(shù)選擇為:字典原子維數(shù)n=80,字典原子個數(shù)K =270,訓練樣本原子數(shù)N = 350,稀疏度L=12,迭代次數(shù)I =32。
將齒輪全壽命數(shù)據(jù)輸入模型中,得到齒輪性能退化評估結(jié)果如圖11所示。實線為故障程度指標DI 值,點劃線為3σ 自適應報警閾值。從圖11中可以看出,前40個樣本即前4 d 的DI 值較低,齒輪處于正常狀態(tài)。曲線在第41個樣本處超過了一級報警閾值,并在第41個樣本到第80個樣本的DI 值出現(xiàn)反復波動,可認為此時出現(xiàn)早期故障。第81個樣本處出現(xiàn)第二次增幅,超過了二級報警閾值,表現(xiàn)出故障程度的加深。第101個樣本處退化曲線發(fā)生較大增幅,之后一直處于較高的幅值,表明此時齒輪失效。
圖11 前35 個樣本訓練的齒輪性能退化曲線
將退化曲線與表2的停機觀測結(jié)果進行對比,可以看出退化曲線的3個故障階段與觀測結(jié)果大致吻合。觀測結(jié)果顯示第6 d 結(jié)束時觀測到齒輪出現(xiàn)剝落,而退化曲線在第5 d 的樣本處出現(xiàn)波動異常超過報警閾值,評估結(jié)果比肉眼觀測結(jié)果提前一天發(fā)現(xiàn)早期故障,可能的原因有兩種:一是本文的退化評估方法存在問題出現(xiàn)誤報,或者并不適用于當前所用案例;二是在第5 d 時已經(jīng)出現(xiàn)了剝落裂紋,但由于裂紋較小肉眼不容易觀測,故觀測人員認為第5 d 沒有出現(xiàn)故障。退化曲線在第9 d 時有激增,超過第二報警閾值,而觀測結(jié)果為第8 d 出現(xiàn)第十六齒的早期剝落,可能原因為第8天結(jié)束時才產(chǎn)生剝落,隨著剝落程度的加深,第9 d 的剝落增大使退化曲線出現(xiàn)明顯增幅。
針對上述問題進行進一步分析,如果第5 d 時第二齒就出現(xiàn)了剝落,則認為第41個樣本到第80個樣本都處于同一故障水平,如果使用前50個樣本作為訓練樣本訓練字典,則字典中包含了早期故障成分,可能會使前80個樣本都可得到誤差較小的重構(gòu)而被判定為無故障狀態(tài)。分別使用前40個樣本和前50個樣本進行訓練,得到如圖12和圖13所示性能退化曲線。從圖12和圖13可以看出使用前40個樣本訓練得到的結(jié)果和前35個樣本訓練得到的結(jié)果是一致的,而使用前50個樣本為訓練樣本得到的性能退化曲線前80個樣本的DI 值都處于較低水平,到第81個樣本才發(fā)生激增超過報警閾值,因此表明之前推斷正確,即第50個樣本處已經(jīng)產(chǎn)生早期故障。
圖12 前40個樣本訓練的齒輪性能退化曲線
圖13 前50個樣本訓練的齒輪性能退化曲線
文獻[27]也對該齒輪數(shù)據(jù)進行了分析,將數(shù)據(jù)特征進行LLP降維后使用基于概率相似度的耦合隱馬爾可夫模型進行性能退化評估,評估結(jié)果如圖14所示。文獻[27]中將退化過程也分為了4個階段:第1階段第1 d ~5 d,為無故障階段;第2階段為第6 d~8 d,由于第二齒開始剝落指標出現(xiàn)波動,進入輕微故障階段;第3階段為第9 d ~ 10 d,由于第十六齒的剝落故障程度進一步加深;第4階段為第11 d ~12 d,齒輪嚴重剝落進入失效階段。與本文方法進行對比,本文所得結(jié)果在輕微故障階段的指標波動較為明顯,且在第5 d 時監(jiān)測到早期故障的出現(xiàn),表明本文方法對故障更加敏感。在第十六齒剝落的早期階段,兩種方法的評價結(jié)果都顯示在第9 d 出現(xiàn)剝落,晚于觀測結(jié)果,可能的原因是第十六齒的早期剝落響應被第二齒的剝落所淹沒,表明所提方法對故障位置的檢測還有待加強??傮w而言本文所提方法對故障程度的評估更為準確。
圖14 基于LLP和耦合隱馬爾可夫模型的齒輪退化評估結(jié)果
論文針對基于概率估算模型的評估方法存在的過早飽和現(xiàn)象,將目前運用于故障診斷領(lǐng)域的字典學習引入其中,提出一種基于AR 模型和字典學習相結(jié)合的齒輪性能劣化評估方法。通過全壽命疲勞試驗數(shù)據(jù)分析處理,發(fā)現(xiàn)在整個性能退化曲線中,故障程度指標DR 值與故障發(fā)展趨勢的一致性更好;相比與傳統(tǒng)的時域指標,該方法能更為及時地發(fā)現(xiàn)軸承早期故障,為設(shè)備維護提供更為精確的數(shù)據(jù)基礎(chǔ)。