梁潔琳，丁 康，何國林，2，林慧斌，蔣 飛

（1.華南理工大學(xué)機(jī)械與汽車工程學(xué)院，廣東 廣州 510640；2.琶洲實(shí)驗(yàn)室，廣東 廣州 510335）

引言

滾動(dòng)軸承作為旋轉(zhuǎn)機(jī)械的關(guān)鍵部件，其退化狀態(tài)影響機(jī)械設(shè)備的安全可靠性能。準(zhǔn)確的軸承壽命預(yù)測有利于預(yù)防性維護(hù)系統(tǒng)，避免軸承突然失效而導(dǎo)致設(shè)備發(fā)生致命故障［1-4］。其中，運(yùn)行狀態(tài)劃分和剩余使用壽命預(yù)測是軸承壽命預(yù)測框架中的兩個(gè)重要問題。

目前常用的軸承正常運(yùn)行和損傷運(yùn)行階段劃分的方法為通過合理的閾值準(zhǔn)則選取退化初始時(shí)間（DIT）［5］。Niu 等［6］應(yīng)用國際規(guī)范的軸承閾值基準(zhǔn)線和報(bào)警系數(shù)獲取融合指標(biāo)的DIT。申彥斌等［7］提出一種基于卷積自編碼器的軸承健康指標(biāo)（Health Indicator，HI）構(gòu)建方法，根據(jù)所提HI 的變化規(guī)律人為劃分三個(gè)運(yùn)行階段；Hong 等［8］利用自組織映射網(wǎng)絡(luò)構(gòu)成全壽命退化評估曲線并手動(dòng)設(shè)定正常和早期退化的閾值界限。然而，以上方法均基于大量統(tǒng)計(jì)特性或數(shù)據(jù)經(jīng)驗(yàn)主觀選取，沒有進(jìn)一步驗(yàn)證DIT 結(jié)果的合理性，對于單個(gè)系統(tǒng)而言，所選DIT 可能會滯后或提前，導(dǎo)致檢修延遲或誤預(yù)警。

基于數(shù)據(jù)驅(qū)動(dòng)的軸承壽命預(yù)測方法主要包括機(jī)器學(xué)習(xí)和統(tǒng)計(jì)模型兩種方法。前者結(jié)合大量可用數(shù)據(jù)樣本預(yù)測退化趨勢，適用于處理復(fù)雜機(jī)械系統(tǒng)的RUL 預(yù)測問題［9-10］。但受限于成本和時(shí)間耗費(fèi)，實(shí)際應(yīng)用中難以獲取高質(zhì)量數(shù)據(jù)樣本，增加了RUL 預(yù)測難度。而統(tǒng)計(jì)模型預(yù)測方法的優(yōu)勢在于同時(shí)結(jié)合先驗(yàn)數(shù)學(xué)［11-12］或物理［13-14］模型以及最新數(shù)據(jù)樣本以更新未知參數(shù)，相對簡單易行。其中，基于指數(shù)模型的傳統(tǒng)粒子濾波算法廣泛應(yīng)用于鋰電池、飛行器等［15-16］不同對象的RUL 預(yù)測領(lǐng)域。滾動(dòng)軸承的全壽命退化過程也與指數(shù)模型相似，但由于受損傷演變特性［17］或隨機(jī)誤差的影響，表征軸承退化過程的觀測指標(biāo)會局部偏離全局指數(shù)模型，呈上下非線性波動(dòng)，導(dǎo)致傳統(tǒng)粒子濾波的預(yù)測曲線偏離實(shí)際退化趨勢，降低壽命預(yù)測精度。

綜上，本文提出了聯(lián)合頻域相關(guān)分析和改進(jìn)粒子濾波算法驅(qū)動(dòng)的滾動(dòng)軸承壽命預(yù)測方法?；跐L動(dòng)軸承在退化過程中振動(dòng)響應(yīng)信號頻域特征具有短期相似和長期差異這一特性，構(gòu)建平均相關(guān)系數(shù)曲線，識別合適的早期退化時(shí)間，并結(jié)合小波降噪和希爾伯特包絡(luò)解調(diào)方法驗(yàn)證所選DIT 的準(zhǔn)確性。由于觀測值的局部波動(dòng)，傳統(tǒng)粒子濾波的粒子后驗(yàn)分布也會相應(yīng)偏離真實(shí)分布，導(dǎo)致RUL 預(yù)測精度降低。因此構(gòu)建了一個(gè)同時(shí)考慮全局指數(shù)式退化趨勢與局部波動(dòng)雙重因素的改進(jìn)粒子濾波算法。所提方法通過新構(gòu)造的粒子約束器校正粒子分布，引導(dǎo)粒子的后驗(yàn)概率密度朝真實(shí)退化值更新，實(shí)現(xiàn)更高精度的壽命預(yù)測。

1 ACC-DFPF 滾動(dòng)軸承壽命預(yù)測方法

1.1 基于頻域相關(guān)分析的ACC 區(qū)間劃分

均方根（Root Mean Square，RMS）具有較好的單調(diào)性，且能有效反映振動(dòng)能量，該時(shí)域特征常作為軸承的健康指標(biāo)。然而，RMS 在早期故障時(shí)間前后的幅值變化不明顯，難以直接準(zhǔn)確識別DIT。實(shí)際上，故障特征頻率作為軸承發(fā)生故障的機(jī)理特性之一，頻域幅值和頻率倍頻會在軸承不同的退化狀態(tài)存在差異。基于此，通過構(gòu)建平均相關(guān)系數(shù)（ACC）曲線和所設(shè)閾值識別RMS 的DIT，具體步驟如圖1所示。

圖1 滾動(dòng)軸承狀態(tài)區(qū)間劃分流程Fig.1 Process of rolling bearing state division

a）對滾動(dòng)軸承的全壽命振動(dòng)響應(yīng)時(shí)域信號做滑窗采樣，并相應(yīng)進(jìn)行快速傅里葉變換得到幅值譜Pk(k=1，2，…，T)，T為幅值譜總個(gè)數(shù)；

b）運(yùn)用Pearson 相關(guān)系數(shù)評價(jià)任意兩幅值譜的相關(guān)程度，如下式：

構(gòu)建相關(guān)系數(shù)矩陣M：

式中ρst為第s和第t個(gè)幅值譜之間的相關(guān)系數(shù)，s，t=1，2，…，T；Ps，n，Pt，n分別對應(yīng)第s和第t個(gè)幅值譜的第n條譜線頻率值和分別為 每個(gè)幅 值譜的 譜線數(shù)，第s和第t個(gè)幅值譜的譜線頻率均值。

c）為直觀地反映起始健康信號與后續(xù)每段信號的相關(guān)程度變化，計(jì)算一維ACC 曲線，計(jì)算式為：

結(jié)合3σ 閾值識別DIT，最后用小波降噪和包絡(luò)解調(diào)譜驗(yàn)證結(jié)果的準(zhǔn)確性。

1.2 傳統(tǒng)粒子濾波(PF)算法框架

粒子濾波通過蒙特卡羅方法（Monte Carlo，MC）避免貝葉斯遞推過程中的高維積分計(jì)算，得到粒子后驗(yàn)概率分布和樣本均值，從而估計(jì)非線性系統(tǒng)的狀態(tài)。滾動(dòng)軸承的退化趨勢常用指數(shù)模型描述，其時(shí)變系統(tǒng)的狀態(tài)方程和觀測方程分別為：

式中uk和νk分別對應(yīng)k時(shí)刻未知模型參數(shù)ak和bk的過程誤差；yk和xk分別為k時(shí)刻的系統(tǒng)觀測指標(biāo)和模型狀態(tài)值，兩者存在觀測誤差χk；uk，νk，χk分別服從標(biāo) 準(zhǔn)差為σv，σq，σr的正態(tài)分布

傳統(tǒng)粒子濾波算法基于式（4）進(jìn)行預(yù)測和更新。

a）預(yù)測過程。已知第1 至k-1 時(shí)刻的系統(tǒng)觀測值y1∶k-1和參數(shù)向量Θk-1，采用一組粒子數(shù)為Np的粒子描述k時(shí)刻參數(shù)的先驗(yàn)概率分布P(Θk|Θk-1，y1∶k-1）。基于式（4），第j個(gè)粒子的參數(shù)預(yù)測向量和模型狀態(tài)值分別為：

b）更新過程。第k時(shí)刻獲取觀測值yk時(shí)，采用下式更新和歸一化每個(gè)預(yù)測粒子的權(quán)重

此時(shí)輸出k時(shí)刻的指標(biāo)估計(jì)值

1.3 改進(jìn)粒子濾波驅(qū)動(dòng)的軸承RUL 預(yù)測

由式（6）可知，PF 算法的粒子權(quán)重更新主要依賴最新觀測值。當(dāng)輸入的最后觀測值偏離全局退化值時(shí)，PF 算法后續(xù)預(yù)測的軸承退化軌跡會偏離真實(shí)趨勢，如圖2 所示。A1和分別為實(shí)際觀測值和全局模型值，當(dāng)A1＞時(shí)，PF 預(yù)測曲線比真實(shí)趨勢更早到達(dá)失效閾值，高估了實(shí)際故障嚴(yán)重程度；反之預(yù)測值比真實(shí)趨勢更平緩，以上兩種局部波動(dòng)觀測值均降低了傳統(tǒng)PF 算法的軸承RUL 預(yù)測精度。

圖2 傳統(tǒng)粒子濾波的預(yù)測趨勢Fig.2 Trends predicted by conventional particle filter

因此，在原有PF 方法的更新過程中，新構(gòu)造一個(gè)同時(shí)考慮全局模型參數(shù)和局部觀測值的約束器，形成改進(jìn)粒子濾波。所提約束器引導(dǎo)粒子分布朝向全局退化值更新，抑制局部波動(dòng)觀測值對軸承HI 預(yù)測結(jié)果的干擾。具體步驟為：

1）構(gòu)建k時(shí)刻的粒子損失函數(shù)：

2）k時(shí)刻的權(quán)重系數(shù)λk∈[0，1]表示為：

式中c和δ均為常 數(shù)因子。令x0，k=a0exp(b0k)為k時(shí)刻的全局模型估計(jì)值。

3）第j個(gè)粒子的參數(shù)預(yù)測向量校正為：

利用Θ0和yk計(jì)算損失函數(shù)，目的是通過平滑的全局模型趨勢和輸入觀測值同時(shí)約束粒子分布，減弱輸入樣本對后續(xù)預(yù)測曲線的影響。此外，依據(jù)yk和x0，k的差異程度，權(quán)重系 數(shù)λk可自適 應(yīng)調(diào)整式（10）的兩項(xiàng)權(quán)重。

綜上，所提ACC-DFPF 滾動(dòng)軸承壽命預(yù)測方法流程如圖3 所示，具體說明如下：

圖3 ACC-DFPF 滾動(dòng)軸承壽命預(yù)測方法流程Fig.3 Flowchart of ACC-DFPF method for rolling bearing life prediction

a）采集軸承數(shù)據(jù)，對各段頻譜進(jìn)行相關(guān)分析，構(gòu)建ACC 曲線，結(jié)合所設(shè)閾值識別DIT，完成HI 指標(biāo)的正常和損傷狀態(tài)劃分，并依據(jù)小波降噪包絡(luò)解調(diào)譜的故障頻率成分驗(yàn)證DIT 的準(zhǔn)確性。

b）采用最小二乘法擬合軸承損傷區(qū)間的HI 指標(biāo)，計(jì)算全局模型經(jīng)驗(yàn)參數(shù)Θ0=[a0b0]Τ，Θ0也作為DFPF 算法的初始化參數(shù)輸入DFPF 算法。當(dāng)k時(shí)刻獲取yk時(shí)，采用式（10）～（12）所提約束器校正粒子預(yù)測過程的先驗(yàn)分布，隨后進(jìn)行粒子權(quán)重更新和重采樣，直至再無觀測值輸入時(shí)，粒子校正停止。

c）繼續(xù)采用式（5）和式（8）計(jì)算，輸出軸承從最后觀測值時(shí)刻kh至預(yù)測失效時(shí)刻的HI 指標(biāo)估計(jì)值和估計(jì) 壽命值

2 試驗(yàn)驗(yàn)證

2.1 ACC 劃分區(qū)間方法的分析對比

選擇XJTU-SY 滾動(dòng)軸承全壽命試驗(yàn)數(shù)據(jù)集［19］驗(yàn)證所提方法的有效性。測試平臺的采樣頻率為2.56 kHz，每間隔1 min 采樣1.28 s，在3 類工況下各重復(fù)5 次試驗(yàn)，共計(jì)15 組軸承全壽命時(shí)間序列樣本。其中，bearing3-1 發(fā)生外圈磨損。

以軸承bearing3-1 為例，所得ACC 退化曲線和DIT 點(diǎn)如圖4 所示。初始運(yùn)行時(shí)系統(tǒng)不穩(wěn)定，對應(yīng)ACC 值產(chǎn)生幅度波動(dòng)。健康狀態(tài)階段的軸承平均相關(guān)系數(shù)相對平穩(wěn)，隨后在初始故障位置開始減小，并在損傷階段隨退化程度加深而整體減小。圖4中，ACC 值在500～1500 min 內(nèi)相對穩(wěn)定，由此得到3σ 閾值下限為0.6106，在 第2371 min 首次超過閾值，表明此時(shí)該軸承發(fā)生初始故障。

圖4 ACC 識別起始故障點(diǎn)結(jié)果Fig.4 The DIT identification result with ACC model

進(jìn)一步采用小波降噪和包絡(luò)解調(diào)驗(yàn)證所選DIT的準(zhǔn)確性?？紤]早期軸承故障特征微弱，先用小波降噪減弱信號樣本的噪聲干擾。由圖5 降噪包絡(luò)譜可知，第2369 min 未出現(xiàn)故障特征頻率成分，而第2371 和第2373 min 均出現(xiàn)了與外圈故障特征頻率理論值123.32 Hz 接近的頻譜成分，且第2530 min軸承故障程度嚴(yán)重，故障頻率與倍頻的幅值都增大，驗(yàn)證了軸承在第2371 min 存在早期故障。

圖5 不同時(shí)刻的軸承包絡(luò)譜Fig.5 Envelope spectrums of bearing at different moments

此外，將所提方法與時(shí)域特征RMS、峭度kurtosis 和頻域特征頻譜重心（Frequency Center，F(xiàn)C）識別DIT 方法的結(jié)果進(jìn)行對比。bearing3-1 的RMS 在第2340 min 后均超過上限閾值，如圖6 所示，說明RMS 比ACC 方法提早31 min 推斷軸承發(fā)生早期故障。但圖7 中的第2340 和2350 min 降噪解調(diào)譜未有明顯的故障特征頻率譜峰。此外，RMS整體比ACC 曲線振蕩，且在第2340 min 前有多處異常RMS 超過閾值造成DIT 的錯(cuò)誤預(yù)警。

圖6 RMS 識別起始故障點(diǎn)結(jié)果Fig.6 The DIT identification result with RMS model

圖7 第2340 和2350 min 的降噪包絡(luò)譜Fig.7 Denoising envelope spectrum at 2340 and 2350 min

峭度方法識別DIT 的結(jié)果如圖8 所示，在第2523 min 后均超過上限閾值，推斷該時(shí)刻發(fā)生了軸承初始故障，比ACC 方法延遲152 min，對應(yīng)檢驗(yàn)了包絡(luò)譜有突出的故障頻率譜峰。此外，在第2371 min 前有多處超過所設(shè)閾值的劇烈波動(dòng)峭度，和RMS 方法存在類似的誤預(yù)警現(xiàn)象。

圖8 Kurtosis 識別起始故障點(diǎn)結(jié)果Fig.8 The DIT identification result with kurtosis model

頻譜重心識別DIT 如圖9 所示。分別求?。?，2000］ min 和［2400，2454］ min 上升和 下降區(qū)間的FC 值斜率，交叉點(diǎn)位于第2346 min，初步推斷該時(shí)刻發(fā)生軸承初始故障，相比ACC 方法提前25 min識別了DIT。此外，該軸承的FC 指標(biāo)相較于ACC和以上時(shí)域特征整體單調(diào)性較差。

圖9 FC 識別起始故障點(diǎn)結(jié)果Fig.9 The DIT identification result with FC model

綜上，ACC 方法可以比時(shí)域RMS，kurtosis 和頻域FC 方法更準(zhǔn)確地識別出軸承初始故障點(diǎn)，且ACC 曲線整體更平滑，避免錯(cuò)誤預(yù)警早期故障點(diǎn)的情況，為下一步進(jìn)行軸承剩余壽命預(yù)測提供起始點(diǎn)依據(jù)。

2.2 DFPF 方法預(yù)測軸承RUL 的分析對比

選取包括同一尺寸3 個(gè)工況的bearing1-3，bearing2-1 和bearing3-1 作為算法驗(yàn)證數(shù)據(jù)集。結(jié)合3 類運(yùn)行工況共6 組歷史軸承的退化規(guī)律選取統(tǒng)一的失效閾值。對RMS 健康指標(biāo)進(jìn)行最大、最小值歸一化，避免了不同工況下因軸承HI 數(shù)值差異較大導(dǎo)致失效閾值γ難選擇的問題，如圖10 所示。圖10（b）的RMS 值在0.6 以上均變化陡峭，表明該閾值后軸承加速劣化至完全損壞，因此設(shè)定失效閾值γ=0.6。

圖10 6 組軸承數(shù)據(jù)的HI 曲線Fig.10 Six sets of bearing health indicator curves

基于ACC 方法和γ確定表1 驗(yàn)證數(shù)據(jù)集的DIT和實(shí)際失效時(shí)刻kT，采用最小二乘法擬合對應(yīng)HI 指標(biāo)值得到全局經(jīng)驗(yàn)?zāi)Ｐ蛥?shù)a0，b0以及標(biāo)準(zhǔn)差σq，σr和σv。其中，初始權(quán)重系數(shù)取λ0=1 表示沒有觀測值輸入時(shí)，粒子的校正完全依賴所得全局經(jīng)驗(yàn)?zāi)Ｐ蛥?shù)。DFPF 算法的主要結(jié)構(gòu)參數(shù)如表2 所示。

表1 驗(yàn)證數(shù)據(jù)集的DIT 和失效時(shí)刻Tab.1 DIT and actual degradation time of the testing data

表2 DFPF 算法主要結(jié)構(gòu)參數(shù)Tab.2 The main structural parameters of DFPF algorithm

分別輸入特定觀測值比例，采用所提DFPF 算法得到一系列指標(biāo)預(yù)測值，與PF 算法結(jié)果相比，如圖11 所示。

圖11 DFPF 和PF 算法所得性能退化預(yù)測曲線Fig.11 Performance degradation prediction curves obtained by DFPF and PF algorithm

對于軸承bearing2-1，50%的樣本輸入比例對應(yīng)最后輸入觀測值的時(shí)刻約為第469 min，該局部觀測數(shù)據(jù)可能受測量過程的噪聲隨機(jī)影響或是退化演變波動(dòng)特性，偏離了全局經(jīng)驗(yàn)?zāi)Ｐ汀Ｓ捎赑F算法輸出的預(yù)測值主要依賴于最后輸入的觀測值，該預(yù)測失效時(shí)刻比實(shí)際失效時(shí)刻提早12 min，而所提DFPF 方法更為有效，僅提早了5 min。原因在于所提DFPF方法有全局模型參數(shù)和局部觀測值這兩項(xiàng)共同約束粒子分布，降低了后續(xù)預(yù)測結(jié)果對局部波動(dòng)數(shù)據(jù)的敏感性，使得軸承未來狀態(tài)的預(yù)測值比傳統(tǒng)PF 算法更趨近真實(shí)觀測曲線。同理對于軸承bearing1-3 和bearing3-1，DFPF 預(yù)測曲線的95%置信區(qū)間間距更窄，RUL 相對誤差更小，具體如表3 所示。由以上分析可知，所提DFPF 算法預(yù)測軸承剩余壽命精度優(yōu)于PF 算法。

表3 特定時(shí)間段下軸承RUL 預(yù)測誤差結(jié)果Tab.3 Results of bearing RUL prediction error under a specific time period

為了更直觀地說明不同時(shí)間段下所提DFPF 算法的預(yù)測優(yōu)勢，按5%，10%，15%，…，95%共19 組不同的觀測值輸入比例預(yù)測軸承剩余使用壽命，如圖12 所示。

圖12 不同起始預(yù)測時(shí)刻的軸承RUL 曲線Fig.12 RUL curves of bearings at different initial prediction moments

此外，分別用均方根誤差Ermse、平均絕對百分比誤差Emape和平均絕對誤差Emae指標(biāo)評價(jià)兩種算法的RUL 預(yù)測精度，具體計(jì)算公式為：

以上指標(biāo)越小代表預(yù)測精度越高，具體結(jié)果對比如表4 所示。

表4 不同起始預(yù)測時(shí)刻的軸承RUL 誤差Tab.4 RUL errors for bearings at different initial prediction moments

當(dāng)輸入觀測值比例為較少的5%時(shí)，軸承RUL預(yù)測結(jié)果對應(yīng)輸出圖12（a）～（c）的左起第一個(gè)點(diǎn)，可直觀看出其接近實(shí)際值。其原因可能是預(yù)測結(jié)果受初始模型參數(shù)影響，而以上軸承的初始模型參數(shù)均由對應(yīng)的全局指標(biāo)數(shù)據(jù)擬合得出，準(zhǔn)確度較高，因此即使觀測樣本較少，RUL 預(yù)測結(jié)果仍會接近真實(shí)壽命。此外，對于bearing2-1 和bearing1-3，占比為10%至30%對應(yīng)的預(yù)測RUL 也接近真實(shí)值，進(jìn)一步分析原因可能是該比例下的實(shí)際觀測值波動(dòng)較小，整體較接近全局經(jīng)驗(yàn)?zāi)Ｐ汀Ｒ陨戏治隹芍谳斎胗^測值個(gè)數(shù)不足的條件下，觀測值和初始模型參數(shù)會影響軸承RUL 預(yù)測。

3 結(jié)論

（1）構(gòu)造了表征軸承退化性能的ACC 曲線，準(zhǔn)確地定位了軸承起始故障點(diǎn)，并采用小波降噪和希爾伯特包絡(luò)解調(diào)方法輔助驗(yàn)證所選點(diǎn)的早期故障情況。由XJTU-SY 軸承數(shù)據(jù)集的驗(yàn)證結(jié)果對比可知，所提ACC 方法相較于均方根和峭度方法能夠更準(zhǔn)確找到早期故障點(diǎn)，實(shí)現(xiàn)軸承正常與損傷兩個(gè)狀態(tài)階段的有效劃分。

（2）所提DFPF 算法構(gòu)造了含有全局模型參數(shù)和局部實(shí)際觀測值的約束器用于校正粒子分布，克服傳統(tǒng)粒子濾波算法受觀測值波動(dòng)干擾的缺陷。XJTU-SY 軸承數(shù)據(jù)集的試驗(yàn)結(jié)果表明，所提DFPF算法后續(xù)預(yù)測的健康指標(biāo)值比PF 算法更接近全局運(yùn)行狀態(tài)，置信區(qū)間寬度和相對誤差Erh均更小，實(shí)現(xiàn)了更準(zhǔn)確的軸承壽命估計(jì)。

（3）輸入不同數(shù)目觀測數(shù)據(jù)的試驗(yàn)結(jié)果表明，DFPF 算法可以有效減少軸承局部異常波動(dòng)觀測數(shù)據(jù)對軸承RUL 預(yù)測結(jié)果的影響，其RUL 的Ermse，Emape和Emae誤差均小于PF 算法，實(shí)現(xiàn)了更可靠的壽命預(yù)測。此外，受限于試驗(yàn)樣本的數(shù)量，全局經(jīng)驗(yàn)?zāi)Ｐ蛥?shù)的調(diào)整缺少更多全局退化趨勢相似的軸承離線觀測數(shù)據(jù)，而目前在實(shí)際應(yīng)用中難以獲取。今后如何根據(jù)有限的測量數(shù)據(jù)獲取該參數(shù)信息，是壽命預(yù)測研究中需要重點(diǎn)考慮的問題之一。