基于威布爾分布及最小二乘支持向量機的滾動軸承退化趨勢預(yù)測

陳 昌，湯寶平，呂中亮

（重慶大學(xué) 機械傳動國家重點實驗室，重慶 400030）

準(zhǔn)確預(yù)測滾動軸承退化趨勢對預(yù)防設(shè)備性能退化及失效意義重大。滾動軸承趨勢預(yù)測關(guān)鍵即建立合適的性能退化指標(biāo)及預(yù)測模型［1］。單個時頻域特征指標(biāo)［2－4］對早期故障敏感程度低，不能較好確定初始損傷時間。融合時頻域指標(biāo)［5－7］雖能全面反應(yīng)軸承退化信息，但會增加運算的復(fù)雜程度。因不同類型故障的發(fā)生時間具有間歇性，對振動數(shù)據(jù)進(jìn)行統(tǒng)計分析可避免故障類型誤診斷［8］。文獻(xiàn)［8－9］提取的威布爾分布均值、方差及極大似然函數(shù)負(fù)值作為狀態(tài)特征信息識別軸承故障及故障類型，能較好刻畫軸承的運行狀態(tài)；但該統(tǒng)計量對軸承故障發(fā)展的靈敏性較差，而威布爾分布形狀參數(shù)值對軸承運行狀態(tài)變化較敏感，且能隨軸承故障的發(fā)展呈明顯上升趨勢，故用威布爾分布形狀參數(shù)作為性能退化指標(biāo)。

趨勢預(yù)測另一重要步驟為在確定性能退化指標(biāo)后建立可靠的預(yù)測模型。文獻(xiàn)［5，10－11］分別采用神經(jīng)網(wǎng)絡(luò)及灰色理論建立預(yù)測模型，預(yù)測滾動軸承的退化趨勢。由于灰色理論缺乏平穩(wěn)性檢驗，且忽略非平穩(wěn)信息，出現(xiàn)較大預(yù)測誤差。神經(jīng)網(wǎng)絡(luò)存在局部極小、小樣本推廣能力差等難解決問題。而最小二乘支持向量機［6］在小樣本、高維、非線性等數(shù)據(jù)空間下具有較好的泛化能力，但其核參數(shù)σ及正規(guī)化參數(shù)γ對模型預(yù)測、推廣性能影響較大，用粒子群算法優(yōu)化核參數(shù)及正規(guī)化參數(shù)可明顯提高預(yù)測精度及預(yù)測結(jié)果的穩(wěn)定性。

本文將威布爾分布形狀參數(shù)與粒子群優(yōu)化的最小二乘支持向量機相結(jié)合，提出軸承退化趨勢預(yù)測新方法。該方法充分發(fā)揮威布爾分布形狀參數(shù)在早期故障敏感度及粒子群優(yōu)化的最小二乘支持向量機預(yù)測模型優(yōu)勢，預(yù)測精度較高。

1 軸承性能退化指標(biāo)

滾動軸承出現(xiàn)早期故障時其信息較易淹沒于噪聲信號中。希爾伯特變換［8］對軸承早期故障信息提取效果良好，故在提取軸承性能退化指標(biāo)前先對振動數(shù)據(jù)作希爾伯特變換。對正常軸承及不同類型故障軸承而言，振動信號的包絡(luò)值均可用威布爾分布建立模型［12］。威布爾分布概率密度函數(shù)［8］為

據(jù)牛頓迭代法求參數(shù)β，η。由于極大似然方法估計參數(shù)時存在誤差，采用修正方法［13］，修正形狀參數(shù)β可減少估計誤差。修正公式為

式中：N為數(shù)據(jù)長度；βU為修正值；β為估計值。

由Cincinnati大學(xué)實測的滾動軸承全壽命實驗數(shù)據(jù)計算威布爾分布的兩模型參數(shù)、極大似然函數(shù)負(fù)值（nLog）［9］及威布爾分布均值（mu）與方差（Std）［8］，同時計算時域特征值（均值、峭度、波形指標(biāo)、裕度指標(biāo)）、頻域特征值（均值頻率、均方根頻率）、融合指標(biāo)（主成分分析 （PCA）第一主成分（PC1）［6］）。各特征值趨勢曲線見圖1、圖2。由兩圖看出，威布爾分布均值及方差、尺度參數(shù)η及時域特征量均值、均方根頻率不隨故障發(fā)展呈上升趨勢，不適合作軸承的性能退化指標(biāo)。而波形指標(biāo)、裕度指標(biāo)雖呈現(xiàn)一定退化上升趨勢，但指標(biāo)總體信息較嘈雜，亦不適合作軸承的衰退性能指標(biāo)。峭度、均值頻率指標(biāo)隨故障發(fā)展呈較好的上升趨勢，但在700點后才呈現(xiàn)上升趨勢，此時軸承已趨向嚴(yán)重失效狀態(tài)，不能反應(yīng)軸承的早期故障。威布爾分布形狀參數(shù)及極大似然函數(shù)負(fù)值在約500點呈上升趨勢，與PC1退化趨勢相近，能較好、全面反映軸承的退化趨勢，且對早期故障的敏感度更高。與極大似然函數(shù)負(fù)值相比，威布爾分布形狀參數(shù)在500點前基本趨于一直線，整體趨勢較純凈。故選威布爾分布的形狀參數(shù)作為性能退化指標(biāo)。

圖1 全壽命軸承信號時、頻域指標(biāo)Fig.1 The time and frequency domain character indicators of the whole life bearing test signal

圖2 全壽命軸承信號統(tǒng)計指標(biāo)Fig.2 The statistical indicators of the whole life bearing test signal

用威布爾分布形狀參數(shù)作為軸承性能退化指標(biāo)具有兩點優(yōu)勢：① 計算簡單；② 對軸承早期故障的敏感程度高，能較好反映軸承性能退化趨勢。

2 基于粒子群優(yōu)化的最小二乘支持向量機

2.1 最小二乘支持向量機理論

最小二乘支持向量機將最小二乘線性理論引入支持向量機，為標(biāo)準(zhǔn)支持向量機理論的擴展利用［6］。LSSVM定義與標(biāo)準(zhǔn)支持向量機不同的約束函數(shù)，并據(jù)拉格朗日函數(shù)及KKT條件將不等式約束化成等式約束，構(gòu)成LSSVM優(yōu)化模型為

式中：ω為權(quán)向量；γ為正則化參數(shù)或曰懲罰因子，其可決定對超出誤差樣本的懲罰程度，為支持向量機擬合程度與推廣能力的平衡參數(shù)；b為偏差向量；ξi為松弛變量，表明實際對象對逼近函數(shù)在樣本數(shù)據(jù)點的誤差期望。

由式（5）得最小二乘支持向量機的拉格朗日函數(shù)為

式中：αi為拉格朗日乘子。

對式（7）進(jìn)行優(yōu)化，即 L（ω，b，e，α）分別對 ω，ei，b，αi求偏導(dǎo)數(shù)，并令偏導(dǎo)數(shù)為零。消去變量ω，e得

用徑向基函數(shù)（RBF）作為核函數(shù)，表達(dá)式為

式中：σ為核函數(shù)寬度。

2.2 粒子群優(yōu)化的最小二乘支持向量機

LSSVM雖能較好解決小樣本、非線性、局部極小點等問題，但實際應(yīng)用中LSSVM核參數(shù)σ及正規(guī)化參數(shù)γ對其預(yù)測推廣性能影響較大，因此用LSSVM預(yù)測時需對核參數(shù)及正規(guī)化參數(shù)進(jìn)行優(yōu)化。LSSVM優(yōu)化方法較多，如交叉驗證法、梯度下降法、遺傳算法、粒子群算法等。粒子群算法具有能獲取全局最優(yōu)結(jié)果、收斂速度快、魯棒性好、尋優(yōu)能力強及編程易實現(xiàn)等優(yōu)點［14］。故選粒子群算法優(yōu)化最小二乘支持向量機的模型參數(shù)。

選均方誤差（MSE）作為PSO算法的粒子適應(yīng)度函數(shù)，表達(dá)式為

式中：n為訓(xùn)練樣本個數(shù)；yi為實際值；y＾i為預(yù)測值。

PSO優(yōu)化LSSVM算法中，將尋優(yōu)目標(biāo)參數(shù)γ、σ作為粒子，均方誤差作為適應(yīng)度函數(shù)，算法基本步驟為：① 初始化種群參數(shù)，包括種群規(guī)模、最大粒子數(shù)及終止條件等；隨機設(shè)置參數(shù)γ、σ的初始位置及速度。② 作一次LSSVM訓(xùn)練，計算訓(xùn)練樣本均方誤差作為粒子的初始適應(yīng)值。③ 更新速度及位置向量。④ 將更新的γ、σ值重新代入SVM模型，并據(jù)②重新訓(xùn)練，保存其輸出結(jié)果，計算粒子的適應(yīng)值。⑤ 將④中所得適應(yīng)值與當(dāng)前粒子適應(yīng)值比較，若優(yōu)于當(dāng)前粒子適應(yīng)值，則更新當(dāng)前粒子適應(yīng)值，并將參數(shù)γ、σ更新為與④中相對應(yīng)的γ、σ。⑥ 判斷是否滿足停止條件，滿足則終止迭代，否則返回③。⑦ 輸出最優(yōu)解，結(jié)束算法。

3 基于威布爾分布形狀參數(shù)及粒子群優(yōu)化的LSSVM軸承趨勢預(yù)測

本文用威布爾分布形狀參數(shù)作為性能退化指標(biāo)，利用LSSVM作趨勢預(yù)測。由于早期故障易淹沒于振動信號中，造成趨勢預(yù)測精度下降，故引入Hilbert變換對振動信號進(jìn)行預(yù)處理。對Hilbert處理后全壽命振動數(shù)據(jù)提取威布爾分布形狀參數(shù)，將其作為性能退化指標(biāo)，輸入LSSVM模型中預(yù)測軸承的退化趨勢；用粒子群算法優(yōu)化LSSVM的核參數(shù)σ及正規(guī)化參數(shù)γ提高預(yù)測精度。具體流程見圖3：

（1）先對振動數(shù)據(jù)作Hilbert變換，再計算其威布爾分布的形狀參數(shù)；

（2）將該形狀參數(shù)作為性能退化指標(biāo)，輸?shù)絃SSVM模型中；

（3）將粒子群優(yōu)化的核參數(shù)及正規(guī)化參數(shù)作為LSSVM模型參數(shù)建立預(yù)測模型，獲取預(yù)測結(jié)果；

（4）按步驟（1）、（2）處理軸承全壽命實驗數(shù)據(jù)，將結(jié)果輸入預(yù)測模型中進(jìn)行預(yù)測。

圖3 滾動軸承退化趨勢預(yù)測流程Fig.3 The flowchart of the rolling bearing degradation trend forecasting process

4 基于軸承全壽命實驗驗證

用實測的滾動軸承全壽命數(shù)據(jù)［15］進(jìn)行驗證，全壽命實驗裝置及采集儀器布置見圖4。軸承實驗臺轉(zhuǎn)軸上安裝4個Rexnord公司ZA－2115雙列滾子軸承，交流電機通過帶傳動2 000 r／min恒定轉(zhuǎn)速帶動轉(zhuǎn)軸旋轉(zhuǎn)，實驗時為軸承施加6 000 lbs徑向載荷。每個軸承的X，Y向各安裝一PCB 353B33加速度傳感器，采樣頻率20 kHz，采樣間隔10 min，采樣長度20 480點。軸承持續(xù)運行7天，直到失效。

圖4 軸承全壽命實驗裝置及采集儀器布置圖［15］Fig.4 The whole life bearing test equipment

選威布爾分布形狀參數(shù)作為軸承性能退化指標(biāo)，該指標(biāo)部分?jǐn)?shù)據(jù)用于訓(xùn)練LSSVM模型并完成預(yù)測，將所得預(yù)測值與實際性能退化指標(biāo)進(jìn)行對比分析。確定性能退化指標(biāo)后利用粒子群算法優(yōu)化最小二乘支持向量機的核參數(shù)σ及正規(guī)化參數(shù)γ。粒子群數(shù)目設(shè)為20，最大迭代次數(shù)設(shè)為200，加速因子 c1，c2均設(shè)為1.5，粒子群算法優(yōu)化最小二乘支持向量機核參數(shù)σ為24、正規(guī)化參數(shù)γ為94.5，用此兩模型參數(shù)構(gòu)建最小二乘支持向量機模型進(jìn)行訓(xùn)練及預(yù)測。在700點之前軸承基本處于正常運行狀態(tài)。因此用700～900點進(jìn)行訓(xùn)練，900點后的數(shù)據(jù)作預(yù)測。為評價預(yù)測結(jié)果的準(zhǔn)確性，用平均相對誤差作為預(yù)測效果［1］的評價指標(biāo)。

式中：n為測試數(shù)據(jù)個數(shù)；yi為實際值；y＾i為預(yù)測值。

預(yù)測中預(yù)測模型保持不變，將預(yù)測結(jié)果作為輸入進(jìn)行下一步預(yù)測，獲得預(yù)測結(jié)果后再次作為輸入進(jìn)行下一步預(yù)測，如此循環(huán)迭代實現(xiàn)預(yù)測。該預(yù)測方法優(yōu)勢在于預(yù)測結(jié)果可作為輸入用于預(yù)測之后結(jié)果，可用于新數(shù)據(jù)較難獲取情況。預(yù)測數(shù)據(jù)部分結(jié)果見圖5。由圖5看出，20步預(yù)測范圍內(nèi)預(yù)測曲線與實際性能退化曲線較接近，平均相對誤差為0.158 2。與實測數(shù)據(jù)作為輸入預(yù)測方法對比知，平均相對誤差為0.114 9。通常軸承的全壽命實驗數(shù)據(jù)較難獲取，用預(yù)測數(shù)據(jù)作為下一步輸入可使預(yù)測更具實用性、推廣性。

為進(jìn)一步驗證形狀參數(shù)作為退化指標(biāo)方法的準(zhǔn)確性及整體算法各環(huán)節(jié)的必要性，進(jìn)行對比驗證，比較預(yù)測結(jié)果的平均相對誤差值：① 峭度作為性能退化指標(biāo)，PSO優(yōu)化LSSVM建立預(yù)測模型；② 均方根值作為性能退化指標(biāo)，PSO優(yōu)化LSSVM建立預(yù)測模型；③ 裕度指標(biāo)作為性能退化指標(biāo)，PSO優(yōu)化LSSVM建立預(yù)測模型；④ PC1作為性能退化指標(biāo)，PSO優(yōu)化LSSVM建立預(yù)測模型；⑤ 形狀參數(shù)作性能退化指標(biāo)，人為固定SVM模型參數(shù)建立預(yù)測模型；⑥ 形狀參數(shù)作性能退化指標(biāo)，粒子群算法優(yōu)化SVM模型參數(shù)建立預(yù)測模型；⑦ 形狀參數(shù)作性能退化指標(biāo)，固定LSSVM的模型參數(shù)建立預(yù)測模型；⑧ 形狀參數(shù)作性能退化指標(biāo)，采用神經(jīng)網(wǎng)絡(luò)建立預(yù)測模型。對比結(jié)果見圖6及表1。

圖5 滾動軸承趨勢預(yù)測曲線與實際狀態(tài)趨勢曲線Fig.5 The actual state trend curve and the forecast trend curve

圖6 其它方法趨勢預(yù)測曲線與實際狀態(tài)趨勢曲線Fig.6 The forecast trend curve obtained by other methods and the actual state trend curve

由圖5、圖6、表1看出，形狀參數(shù)作為性能退化指標(biāo)能較好反應(yīng)軸承的退化趨勢，且預(yù)測效果優(yōu)于峭度、裕度指標(biāo)及均方根值作為性能退化指標(biāo)的預(yù)測效果。而基于峭度、裕度指標(biāo)及均方根等特征指標(biāo)由于對早期故障不敏感，導(dǎo)致預(yù)測結(jié)果出現(xiàn)較大誤差。形狀參數(shù)作為性能退化指標(biāo)的預(yù)測效果與PC1預(yù)測效果相近，且無需多特征信息融合，可減少計算復(fù)雜程度及運算量。采用PSO優(yōu)化LSSVM建立預(yù)測模型預(yù)測效果優(yōu)于PSO優(yōu)化SVM預(yù)測效果，采用PSO優(yōu)化LSSVM模型參數(shù)預(yù)測效果優(yōu)于固定LSSVM模型參數(shù)預(yù)測效果。表明減少模型參數(shù)個數(shù)，合理選擇模型參數(shù)可減少預(yù)測誤差、提高預(yù)測精度。因此，威布爾分布形狀參數(shù)作為性能退化指標(biāo)能有效反映軸承性能退化趨勢，預(yù)測精度較高。

表1 各方法預(yù)測結(jié)果與實際性能退化指標(biāo)間誤差比較Tab.1 The prediction error comparison of different methods

工程實際中軸承趨勢預(yù)測通常用一個軸承訓(xùn)練建模預(yù)測另一同型號、同工況軸承的退化趨勢。試驗中2、3號為同型號、工況類似的兩滾動軸承，用2號軸承訓(xùn)練模型預(yù)測3號軸承退化趨勢，形狀參數(shù)作為性能退化指標(biāo)所得3號軸承預(yù)測結(jié)果的平均相對誤差為0.235 8，而PC1作性能退化指標(biāo)所得平均相對誤差為0.522 9。反之，用3號軸承訓(xùn)練模型預(yù)測2號軸承退化趨勢，形狀參數(shù)作為性能退化指標(biāo)所得2號軸承的平均相對誤差為0.1496，而PC1作性能退化指標(biāo)所得平均相對誤差為0.271 4。實驗表明，威布爾分布的形狀參數(shù)作性能退化指標(biāo)預(yù)測結(jié)果最優(yōu)。

5 結(jié) 論

（1）本文通過分析常用時頻域特征作為性能退化指標(biāo)存在問題，提出采用威布爾分布的形狀參數(shù)作為性能退化指標(biāo)，在識別早期故障的同時亦能較好預(yù)測軸承的性能退化趨勢。

（2）為獲得更準(zhǔn)確預(yù)測結(jié)果，選粒子群算法優(yōu)化LSSVM模型參數(shù)，利用粒子群算法的尋優(yōu)能力提高LSSVM預(yù)測精度。

（3）通過對比分析全壽命振動數(shù)據(jù)，表明本文所提趨勢預(yù)測方法的有效性及各環(huán)節(jié)間關(guān)聯(lián)性，能充分發(fā)揮各部分優(yōu)勢，實現(xiàn)較高精度的退化趨勢預(yù)測。

［1］申中杰，陳雪峰，何正嘉，等.基于相對特征和多變量支持向量機的滾動軸承剩余壽命預(yù)測［J］.機械工程學(xué)報，2013，49（2）：183－189.SHEN Zhongjie，CHEN Xuefeng，HE Zhengjia，et al.Remaining life predictions of rolling bearing based on relative features and multivariable support vector machine［J］.Chinese Journal of Mechanical Engineering，2013，49（2）：183－189.

［2］Gebraeel N，Lawley M，Liu R，et al.Residual life predictions from vibrationbased degradation signals：a neural network approach［J］.Industrial Electronics，IEEE Transactions on，2004，51（3）：694－700.

［3］Shao Y，Nezu K.Prognosis of remaining bearing life using neural networks［J］.Proceedings of the Institution of Mechanical Engineers，Part I：Journal of Systems and Control Engineering，2000，214（3）：217－230.

［4］Janjarasjitt S，Ocak H，Loparo K A.Bearing condition diagnosis and prognosis using applied nonlinear dynamical analysis of machine vibration signal［J］.Journal of Sound and Vibration，2008，317（1／2）：112－126.

［5］奚立峰，黃潤青，李興林，等.基于神經(jīng)網(wǎng)絡(luò)的球軸承剩余壽命預(yù)測［J］.機械工程學(xué)報，2007，43（10）：137－143.XI Lifeng，HUANG Runqing，LI Xinglin，et al.Residual life predictions for ball bearing based on neural networks［J］.Chinese Journal of Mechanical Engineering，2007，43（10）：137－143.

［6］劉寶英，楊仁剛.基于主成分分析的最小二乘支持向量機短期負(fù)荷預(yù)測模型［J］.電力自動化設(shè)備，2008，28（11）：13－17.LIU Baoying，YANG Rengang.Shortterm load forecasting model based on LSSVM with PCA［J］.Electric Power Automation Equipment，2008，28（11）：13－17.

［7］Liao H，Zhao W，Guo H.Predicting remaining useful life of an individual unit using proportional hazards model and logistic regression model［C］.Reliability and Maintainability Symposium，2006：127－132.

［8］Jiménez G A，Munoz A O，DuarteMermoud M A.Fault detection in induction motors using Hilbert and wavelet transforms［J］.Electrical Engineering，2007，89（3）：205－220.

［9］Abbasion S，Rafsanjani A，F(xiàn)arshidianfar A，et al.Rolling element bearings multifault classification based on the wavelet denoising and support vector machine［J］.Mechanical Systems and Signal Processing，2007，21（7）：2933－2945.

［10］Yan J，Guo C，Wang X.A dynamic multiscale Markov model based methodology for remaining life prediction［J］.Mechanical Systems and Signal Processing，2011，25（4）：1364－1376.

［11］徐東，徐永成，陳循，等.基于 EMD的灰色模型的疲勞剩余壽命預(yù)測方法研究［J］.振動工程學(xué)報，2011，24（1）：104－110.XU Dong，XU Yongcheng，CHEN Xun，et al.Residual fatigue life prediction based on grey model and EMD［J］.Journal of Vibration Engineering，2011，24（1）：104－110.

［12］龍凱，李剛成.基于 Hilbert變換與威布爾分布的軸承早期故障診斷［J］.科學(xué)技術(shù)與工程，2010，10（17）：4163－4167.LONG Kai，LI Gangcheng.A method of incipient fault diagnosis for rolling bearing based on hilbert transform and weibull distribution ［J］. Science Technology and Engineering，2010，10（17）：4163－4167.

［13］Ross R.Formulas to describe the bias and standard deviation of the mlestimated weibull shape parameter［J］.Dielectrics and Electrical Insulation， IEEE Transactions on，1994，1（2）：247－253.

［14］Eberhart R，Kennedy J.A new optimizer using particle swarm theory［C］.Pro Sixth International Symposium on Micro Machine and Human Science，Japan，1995：39－43.

［15］Lee J，Qiu H，Yu G，et al.Rexnord technical services，“bearing data set”［R］.IMS，University of Cincinnati，NASA Ames Prognostics Data Repository，＜http：／／ti.arc.nasa. gov／project／ prognosticsdatarepository ＞， NASA Ames，Moffett Field，CA，2007.