陳 昌，湯寶平，呂中亮

（重慶大學(xué) 機(jī)械傳動(dòng)國家重點(diǎn)實(shí)驗(yàn)室，重慶 400030）

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

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

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

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

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

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

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

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

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

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

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

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

2.1 最小二乘支持向量機(jī)理論

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 基于軸承全壽命實(shí)驗(yàn)驗(yàn)證

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

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

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

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

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

為進(jìn)一步驗(yàn)證形狀參數(shù)作為退化指標(biāo)方法的準(zhǔn)確性及整體算法各環(huán)節(jié)的必要性，進(jìn)行對(duì)比驗(yàn)證，比較預(yù)測結(jié)果的平均相對(duì)誤差值：① 峭度作為性能退化指標(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ù)測模型。對(duì)比結(jié)果見圖6及表1。

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

圖6 其它方法趨勢預(yù)測曲線與實(shí)際狀態(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)由于對(duì)早期故障不敏感，導(dǎo)致預(yù)測結(jié)果出現(xiàn)較大誤差。形狀參數(shù)作為性能退化指標(biāo)的預(yù)測效果與PC1預(yù)測效果相近，且無需多特征信息融合，可減少計(jì)算復(fù)雜程度及運(yùn)算量。采用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ù)個(gè)數(shù)，合理選擇模型參數(shù)可減少預(yù)測誤差、提高預(yù)測精度。因此，威布爾分布形狀參數(shù)作為性能退化指標(biāo)能有效反映軸承性能退化趨勢，預(yù)測精度較高。

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

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

5 結(jié) 論

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

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

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

［1］申中杰，陳雪峰，何正嘉，等.基于相對(duì)特征和多變量支持向量機(jī)的滾動(dòng)軸承剩余壽命預(yù)測［J］.機(jī)械工程學(xué)報(bào)，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］.機(jī)械工程學(xué)報(bào)，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］劉寶英，楊仁剛.基于主成分分析的最小二乘支持向量機(jī)短期負(fù)荷預(yù)測模型［J］.電力自動(dòng)化設(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］.振動(dòng)工程學(xué)報(bào)，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.