張雨琦,鄒金慧,馬 軍
(1.昆明理工大學信息工程與自動化學院,云南 昆明 650500;2.云南省礦物管道輸送工程技術(shù)研究中心,云南 昆明 650500;3.昆明理工大學機電工程學院,云南 昆明 650500)
滾動軸承是旋轉(zhuǎn)機械中用途最為廣泛的零部件,其剩余壽命的長短與設(shè)備運行安全、運行工況間存在著直接關(guān)聯(lián),由其引發(fā)的故障是造成設(shè)備失效的重要原因。因此,對滾動軸承剩余壽命進行預測對于預防設(shè)備失效和實現(xiàn)基于狀態(tài)的設(shè)備維護具有重要意義。
目前,剩余壽命預測的方法主要可以分為機理建模和數(shù)據(jù)驅(qū)動建模兩種類型[1]。機理建模的方法主要是指根據(jù)設(shè)備的內(nèi)在運轉(zhuǎn)機理及工作原理來建立設(shè)備的退化模型從而對剩余壽命進行預測[2-3];數(shù)據(jù)驅(qū)動建模的方法主要是指通過數(shù)據(jù)擬合的方式對能反映設(shè)備退化性能的主要性能變量進行擬合,并根據(jù)其變化趨勢來預測設(shè)備的剩余壽命??紤]到滾動軸承的機理具有一定的復雜性,而據(jù)此構(gòu)建機理模型較為困難,因此以數(shù)據(jù)驅(qū)動為支撐的方法在剩余壽命預測中應(yīng)用較為廣泛。文獻[4]以最小量化誤差為衰退指標,利用自組織映射神經(jīng)網(wǎng)絡(luò)對軸承剩余壽命進行預測。文獻[5]選擇改進的EMD以全壽命周期振動信號為對象進行分解,并將不同階段的退化特征量輸入灰色模型進行訓練并對軸承剩余壽命作出預測。文獻[6]選擇飛機空氣制冷機作為研究對象,應(yīng)用譜峭度和最小二乘支持向量機的方法對其壽命趨勢進行預測。
然而,這些模型都存在一定缺陷,只考慮了一個特征參數(shù)或單獨分析了幾個特征值的狀況和趨勢,并沒有將相關(guān)的特征值作整體以及系統(tǒng)性的把握。本文針對此問題,提出了多退化變量灰色預測模型的滾動軸承剩余壽命預測方法。
灰色系統(tǒng)理論是探究某些既存在已知信息,同時也存在未知或無法確知信息的理論和方法。以此理論為基礎(chǔ),可以在離散、有限與雜亂的數(shù)據(jù)中梳理出規(guī)律,并構(gòu)建相應(yīng)的灰色模型,對數(shù)據(jù)的趨勢與發(fā)展狀況進行探究。在灰色預測分析中,最經(jīng)典的模型是單一變量灰色預測模型,即通過一個變量的一階微分方程揭示序列的變化趨勢。但對于復雜工況下的預測,需同時跟蹤多個特征參數(shù)的影響,并對他們進行融合。因此引入多退化變量灰色預測模型來根據(jù)多個變量的變化趨勢對軸承壽命進行預測[7-8]。
(1)
式(1)中,i=1,2,…,n;k=1,2,…,m,m為每個序列的數(shù)據(jù)個數(shù)。
dX(1)/dt=AX(1)+B
(2)
記ai=(ai1,ai2,…,ain,bi)T,i=1,2,…,n,則:
(3)
式(3)中,L=(L1,L2,…,Lj,…Ln,1),
則多變量灰色預測模型的計算值為:
(4)
(5)
(6)
當n=1時,多退化變量灰色預測模型退化為單一變量灰色預測模型。
在實際工程中,當滾動軸承開始工作后,隨著時間的增加,其狀態(tài)會由正常逐漸發(fā)生退化,對于滾動軸承的剩余壽命預測而言,如果從滾動軸承起始正常無故障運轉(zhuǎn)狀態(tài)進行剩余壽命預測,則會出現(xiàn)較大誤差,因此需要對早期故障點即突變點進行識別,從突變點開始進行剩余壽命預測會更為合理?;诖?引入CUSUM算法對突變點進行判斷。
CUSUM算法是通過對事件的突變狀態(tài)進行累積,將過程中的小偏移量累加起來,求其累積和進而判斷是否發(fā)生突變。該算法要求的假定條件較少,能有效反映過程變化的靈敏性,非常適合用于滾動軸承壽命退化過程中的突變檢測。其算法步驟如下所示[9]:
(7)
步驟2 令累積和Ti為:
(8)
式(8)中,i=1,2,…,n。
步驟3 找出Ti中的最大值Tmax,其對應(yīng)的點的橫坐標xmax就是早期故障點即突變發(fā)生時刻:
|Tmax|=max(|Ti|)
(9)
由于單一退化變量灰色預測模型缺乏對能夠表征滾動軸承退化過程的其他變量的分析考量,可能導致有效的信息丟失,使得預測的準確性與精度受到一定的限制,因此本文以對滾動軸承性能退化過程敏感的多個特征參數(shù)為基礎(chǔ),建立多退化變量灰色預測模型對滾動軸承剩余壽命進行預測。其預測流程如下:
1) 選取滾動軸承全壽命周期數(shù)據(jù)并提取對其退化趨勢敏感的均方根、峭度、功率譜密度均值三個指標。
2) 將1)所述的三個指標分別輸入多變量灰色預測模型進行預測,通過各種誤差精度來計算指標權(quán)重,并根據(jù)CUSUM累積和理論判斷出早期故障文件。
3) 從早期故障文件開始,選取相同間隔工作時間文件的三個特征值輸入多退化變量灰色預測模型中進行訓練并得到辨識參數(shù)值。
4) 結(jié)合辨識參數(shù)值對訓練指標進行擬合并預測下一時刻的退化性能指標。
5) 通過多退化變量灰色預測模型建立退化性能指標與壽命之間的非線性映射關(guān)系,并對下一時刻的滾動軸承剩余壽命進行預測。
算法流程圖如圖1所示。
圖1 滾動軸承剩余壽命預測流程圖Fig.1 Flowchart of theresidual life prediction in rolling bearing
筆者采用美國辛辛那提大學智能維護系統(tǒng)的滾動軸承全壽命周期數(shù)據(jù)進行實驗[10]。圖2展示了其實驗平臺狀況,主軸上裝配4個Rexnord ZA-2115雙排列軸承,主軸由直流電機以皮帶為媒介提供動力,每個軸承每排包含16個滾動體,所有軸承均采用油潤滑,轉(zhuǎn)速為2 000 r/min,通過彈性系統(tǒng)向軸和軸承上施加6 000 lb(2 721.5 kg)的徑向載荷,借助NI DAQ 6062 E數(shù)據(jù)采集卡采集數(shù)據(jù),所有數(shù)據(jù)均為加速疲勞實驗。本文選取的數(shù)據(jù)集采樣頻率為20 kHz,數(shù)據(jù)采集間隔時間為10 min,采集數(shù)據(jù)次數(shù)為984次,在每組數(shù)據(jù)中,采樣點的數(shù)量為20 480個,在實驗的末尾,軸承1在持續(xù)運轉(zhuǎn)約163.83 h時,出現(xiàn)了失效問題,具體表現(xiàn)為外圈損傷嚴重,圖3描述了其從正常運轉(zhuǎn)到外圈故障失效的全壽命過程時域波形圖。
圖2 軸承全壽命周期數(shù)據(jù)實驗平臺Fig.2 Experimental platform for lifetime data of rolling bearing
圖3 軸承全壽命周期數(shù)據(jù)時域波形圖Fig.3 Time domain waveform diagram of lifetime data
滾動軸承在運作過程中局部出現(xiàn)故障時,其振動信號中不僅包含故障引起的瞬態(tài)沖擊信息,還存在軸轉(zhuǎn)頻、倍頻等諧波成分和噪聲,以致其振動信號的能量和波形產(chǎn)生一定的變化?;诖?在訓練樣本的選擇上,不能直接以振動信號的振幅特征為準。在實踐操作中,主要的統(tǒng)計特征包括峰值、絕對平均值、均方根、波形系數(shù)、峰態(tài)因子、脈沖因子、裕度因子、偏度、峭度、功率譜密度均值等。結(jié)合現(xiàn)場實踐狀況[11-14],本文選取RMS(均方根值)、峭度、功率譜密度均值三個特征來描述滾動軸承整個壽命周期的變化趨勢:
1) 均方根值:即有效值,由于均方根值是對時間的平均,所以對存在表面裂紋且表現(xiàn)出無規(guī)則振動波形的故障表現(xiàn)出較顯著的敏感特性,針對故障測量值可做出恰當?shù)脑u價[15],描述為:
(10)
2) 峭度:峭度是反映波形尖峰度的歸一化累積量。當滾動軸承在無故障條件下運行時,其所表現(xiàn)出的振動信號幅值基本是符合正態(tài)分布規(guī)律的,峭度指標值K≈3。當滾動軸承局部表現(xiàn)出故障時,因故障帶來的沖擊,振動信號概率密度會呈現(xiàn)出顯著異常的增加,信號幅值分布不會再貼合正態(tài)分布,相應(yīng)的峭度值會顯著增大,由此可以將其作為有效反映軸承是否存在故障的指標[16]。
峭度指標可表示為:
(11)
3) 功率譜密度均值:以功率譜變化(是否有額外譜峰)為依據(jù),可以對故障是否存在作出科學判定。而功率譜密度均值則是對每段數(shù)據(jù)作采集,并依次完成功率譜密度求解,累加平均之后,其結(jié)果即定義為功率譜密度均值[17]。
將自相關(guān)函數(shù)Rx(τ)的傅里葉變換通過下式作具體描述:
(12)
Sx(f)即為x(t)的功率譜密度函數(shù)。則其均值為:
(13)
式(13)中,n=1,2,…,n表示數(shù)據(jù)分段數(shù)。
分別選取上述三個指標:均方根值、峭度值、功率譜密度均值,并計算采集的每組數(shù)據(jù)的三個指標量化到區(qū)間[0,1]的歸一化值(共984組數(shù)據(jù)文件)繪制出如圖4—圖6所示的描述全壽命周期數(shù)據(jù)變化趨勢的各個特征圖。當軸承表面剛開始出現(xiàn)故障時,形成小的剝落或裂紋,均方根、峭度和功率譜密度均值隨之逐漸變大,之后由于連續(xù)的滾動接觸,三個指標均經(jīng)歷一段平滑階段,當損傷擴展到更大范圍時,沖擊再次變大,三個指標持續(xù)攀升且攀升幅度明顯增大,此時軸承出現(xiàn)嚴重故障并使機械設(shè)備無法繼續(xù)運轉(zhuǎn)[10,18]。由圖3—圖6均可知,軸承前期并沒出現(xiàn)故障,如果根據(jù)整體周期數(shù)據(jù)進行預測并不合理,因此需要對最早出現(xiàn)故障的文件標號數(shù)進行判斷。
圖4 均方根特征歸一化趨勢圖Fig.4 Normalization trend of root mean square
圖5 峭度特征歸一化趨勢圖Fig.5 Normalization trend of kurtosis
圖6 功率譜密度均值特征歸一化趨勢圖Fig.6 Normalization trend of mean power spectral density
由于上述三個敏感指標對滾動軸承壽命數(shù)據(jù)退化情況的表征不同,依據(jù)均方誤差(MSE)、平均絕對誤差(MAE)、平均絕對百分比誤差(MAPE)、均方百分比誤差(MSPE)、均方根誤差(RMSE)這五種誤差來計算不同特征值所占權(quán)重[19],分別將三個特征值輸入多退化變量灰色預測模型之中,并根據(jù)權(quán)重占比來綜合判斷其最早出現(xiàn)故障的文件。
(14)
(15)
(16)
(17)
(18)
式(14)—式(18)中,Q0代表原始值,Qm代表預測值,N為樣本組數(shù)。
(19)
式(19)中,i為第i個特征,N代表特征總數(shù),本文取N=3。
由上式分別計算出峭度、均方根值、功率譜密度均值在綜合指標中所占權(quán)重,如表1所示。根據(jù)不同特征的權(quán)重計算出綜合指標并繪制出如圖7所示的經(jīng)平滑處理(Moving Average, MA)后的歸一化綜合指標退化趨勢圖,平滑處理如式(20)所示。
(20)
式(20)中,Ns表示去掉均值后的數(shù)據(jù)總數(shù),s表示存在于子矩陣中的數(shù)據(jù)數(shù),一般情況下為奇數(shù)且s?
Ns,dn表示原始綜合指標中的第n個指標值,man為經(jīng)MA處理后的新指標值。
表1 三個特征在綜合指標中所占權(quán)重表
由文獻[20—22]可知,故障均發(fā)生于400號文件以后,計算出前400號文件的綜合指標均值,如圖7中直線所示。根據(jù)綜合指標線與均值線相交部分的放大圖可知,突變文件大概位于400—600號文件中間,根據(jù)1.2節(jié)介紹的CUSUM算法,以600號文件為界,繪制出累積和曲線圖如圖8所示,可知當文件標號為532號時,曲線發(fā)生突變。
結(jié)合圖9繪制出的532號文件前后文件的包絡(luò)譜變化圖,可以更清楚地顯示出532號文件處故障還未明顯發(fā)生,在533號文件處出現(xiàn)早期故障,在534號文件時故障表示明顯,由此可以判定軸承從533號文件開始發(fā)生故障,即533號文件為早期故障點(對應(yīng)時間是:軸承已經(jīng)工作88.67 h時)。
圖7 歸一化綜合指標退化趨勢圖(經(jīng)平滑處理)Fig.7 Degradation trend of normalized comprehensive index (smoothed)
圖8 累積和曲線圖Fig.8 Cumulative sum curve
圖9 文件包絡(luò)譜Fig.9 Envelope spectrum
由于已判斷出全壽命周期數(shù)據(jù)是從533號文件開始發(fā)生故障,則根據(jù)故障發(fā)生規(guī)律和圖10所示,分別從早期故障點533號文件和第一次出現(xiàn)明顯故障表示(即第一個波形尖端)的704號文件進行壽命預測。從533號文件開始,每隔20個文件取一次數(shù)據(jù),共取16組數(shù)據(jù)的三個退化變量值作為訓練樣本,并將其輸入多退化變量灰色預測模型,分別獲取預測辨識參數(shù),預測之后7組數(shù)據(jù)的壽命值。從704號文件開始,每隔20個文件取一次數(shù)據(jù),共取9組數(shù)據(jù)的三個退化變量值作為訓練樣本,同樣以上述方法預測后6組的壽命值。預測結(jié)果分別如圖11、圖12所示,由圖可知,從早期故障點開始進行預測,其特征值經(jīng)歷了整體的故障變化趨勢,預測結(jié)果與實際壽命值誤差不大,而從明顯出現(xiàn)故障表示的故障點開始預測,由于其在初期的故障后,經(jīng)歷一段時間的平穩(wěn)過渡階段再進入深度故障期(即特征值出現(xiàn)下降再上升的過程),則預測的后期會出現(xiàn)較大誤差,因此本文取早期故障點的特征值為初始特征輸入多退化變量灰色預測模型中進行剩余壽命的預測。
圖10 壽命預測故障點提取圖Fig.10 Fault point extraction for life prediction
圖11 多退化變量灰色預測模型壽命預測對比圖(533號文件起)Fig.11 Life prediction comparison diagram of grey prediction model with multiple degenerate variables (from No.533)
圖12 多變量灰色預測模型壽命預測對比圖(704號文件起)Fig.12 Life prediction comparison diagram of grey prediction model with multiple degenerate variables (from No.704)
為了對比和評估預測結(jié)果,將灰色預測模型(多退化變量與單一變量)、SVR(支持向量回歸)預測模型、BP神經(jīng)網(wǎng)絡(luò)預測模型的分析結(jié)果作進一步對比,對比圖如圖11、圖13、圖14、圖15、圖16所示(由于預測已工作時間可能出現(xiàn)超過實際壽命值的情況,因此對于剩余壽命為負值的情況全部設(shè)置剩余壽命為0),其中圖14為圖13預測曲線部分(即虛線部分)的局部放大圖。分別計算出各預測壽命曲線與實際壽命曲線相應(yīng)的MSE,MAE,MAPE,MSPE,RMSE與NSE,將它們作為評價模型優(yōu)劣的指標。其中NSE指納什系數(shù),該系數(shù)是用來反映模型質(zhì)量的參數(shù)值,若納什系數(shù)越接近于1,則證明模型質(zhì)量越優(yōu)[23],納什系數(shù)表示為:
圖13 單一變量灰色預測與多變量灰色預測對比圖(533號文件起)Fig.13 Comparison diagram between grey prediction with single variable and multiple variables(from No.533)
圖14 單一變量灰色預測與多變量灰色預測對比局部放大圖Fig.14 Local enlarged diagram of comparison between grey prediction with single variable and multiple variables
(21)
圖15 BP神經(jīng)網(wǎng)絡(luò)預測模型壽命預測對比圖(533號文件起)Fig.15 Comparison diagram of BP neural network model for life prediction(from No.533)
圖16 支持向量回歸預測模型壽命預測對比圖(533號文件起)Fig.16 Comparison diagram of SVR model for life prediction(from No.533)
圖17 多退化變量灰色預測曲線95%置信區(qū)間圖Fig.17 95% confidence interval diagram of grey
圖18 BP神經(jīng)網(wǎng)絡(luò)預測曲線95%置信區(qū)間圖Fig.18 95% confidence interval diagram of BP
圖19 SVR預測曲線95%置信區(qū)間圖Fig.19 95% confidence interval diagram of SVR prediction curve with multiple degenerate variables neural network prediction curve prediction curve
預測數(shù)據(jù)組數(shù)1234567實際已使用壽命/h142.000 0145.333 3148.666 7152.000 0155.333 3158.666 7162.000 0單一變量灰色預測(峭度)預測已使用壽命/h140.994 4143.846 0146.605 6149.269 3151.833 2154.294 4156.650 4相對誤差/%0.7081.0231.3861.7972.2532.7563.302單一變量灰色預測(均方根)預測已使用壽命/h140.903 1143.689 1146.368 0148.935 5151.388 2153.723 2155.938 6相對誤差/%0.7721.1311.5462.0162.5393.1163.742單一變量灰色預測(功率譜密度均值)預測已使用壽命/h138.974 8141.231 9143.352 7145.339 4147.195 3148.923 9150.529 3相對誤差/%2.1302.8223.5744.3825.2396.1407.081多變量灰色預測預測已使用壽命/h142.471 3146.143 2149.932 8153.886 8158.063 4162.532 3163.83 相對誤差/%0.3310.5570.8511.2411.7572.4361.129
表3 三種模型壽命預測值與實際壽命值誤差對比表
圖17—圖19分別展示了三種預測模型預測結(jié)果曲線的95%置信區(qū)間。由圖中可知,雖然三種模型的預測結(jié)果均處于95%置信區(qū)間內(nèi),但BP神經(jīng)網(wǎng)絡(luò)預測模型和SVR預測模型所預測出的結(jié)果表現(xiàn)出較大的離散性,而多退化變量灰色預測模型則呈現(xiàn)出較高的精度。單一變量灰色預測模型與多變量灰色預測模型的相對誤差如表2所示。三種預測模型和灰色預測模型中的多退化變量預測與單一退化變量預測對比的各種誤差評價指標計算結(jié)果如表3所示。從表2、表3中可以看出多退化變量灰色預測模型預測結(jié)果能更好地逼近真實壽命值,各種誤差均優(yōu)于單一變量灰色預測和其他兩種預測模型,而能表征模型質(zhì)量的參數(shù)納什系數(shù)則高于其他兩種預測模型。因此,結(jié)合了多個退化變量的灰色預測在小樣本條件下預測軸承壽命更優(yōu)于單一變量灰色預測,而灰色預測模型,則要優(yōu)于SVR預測模型與BP神經(jīng)網(wǎng)絡(luò)預測模型兩類模型。
本文提出了多退化變量灰色預測模型的滾動軸承剩余壽命預測方法。該方法結(jié)合多個表征軸承退化趨勢的特征參數(shù)與早期故障突變點對滾動軸承進行早期故障識別,并利用軸承壽命與特征參數(shù)之間的映射關(guān)系建立多退化變量灰色預測模型對滾動軸承剩余壽命進行預測。仿真實驗結(jié)果表明,本文所提取的參數(shù)集能夠有效表征滾動軸承退化趨勢,多退化變量灰色預測模型的預測精度及可靠性均優(yōu)于單一退化變量灰色預測模型、BP神經(jīng)網(wǎng)絡(luò)預測模型及SVR預測模型,對滾動軸承的剩余壽命研究具有重要的指導意義。