聯合頻域相關分析和改進粒子濾波的滾動軸承壽命預測方法

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

（1.華南理工大學機械與汽車工程學院，廣東 廣州 510640；2.琶洲實驗室，廣東 廣州 510335）

引言

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

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

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

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

1 ACC-DFPF 滾動軸承壽命預測方法

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

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

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

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

b）運用Pearson 相關系數評價任意兩幅值譜的相關程度，如下式：

構建相關系數矩陣M：

式中ρst為第s和第t個幅值譜之間的相關系數，s，t=1，2，…，T；Ps，n，Pt，n分別對應第s和第t個幅值譜的第n條譜線頻率值和分別為 每個幅 值譜的 譜線數，第s和第t個幅值譜的譜線頻率均值。

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

結合3σ 閾值識別DIT，最后用小波降噪和包絡解調譜驗證結果的準確性。

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

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

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

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

a）預測過程。已知第1 至k-1 時刻的系統(tǒng)觀測值y1∶k-1和參數向量Θk-1，采用一組粒子數為Np的粒子描述k時刻參數的先驗概率分布P(Θk|Θk-1，y1∶k-1）?；谑剑?），第j個粒子的參數預測向量和模型狀態(tài)值分別為：

b）更新過程。第k時刻獲取觀測值yk時，采用下式更新和歸一化每個預測粒子的權重

此時輸出k時刻的指標估計值

1.3 改進粒子濾波驅動的軸承RUL 預測

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

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

因此，在原有PF 方法的更新過程中，新構造一個同時考慮全局模型參數和局部觀測值的約束器，形成改進粒子濾波。所提約束器引導粒子分布朝向全局退化值更新，抑制局部波動觀測值對軸承HI 預測結果的干擾。具體步驟為：

1）構建k時刻的粒子損失函數：

2）k時刻的權重系數λk∈[0，1]表示為：

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

3）第j個粒子的參數預測向量校正為：

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

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

圖3 ACC-DFPF 滾動軸承壽命預測方法流程Fig.3 Flowchart of ACC-DFPF method for rolling bearing life prediction

a）采集軸承數據，對各段頻譜進行相關分析，構建ACC 曲線，結合所設閾值識別DIT，完成HI 指標的正常和損傷狀態(tài)劃分，并依據小波降噪包絡解調譜的故障頻率成分驗證DIT 的準確性。

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

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

2 試驗驗證

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

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

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

圖4 ACC 識別起始故障點結果Fig.4 The DIT identification result with ACC model

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

圖5 不同時刻的軸承包絡譜Fig.5 Envelope spectrums of bearing at different moments

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

圖6 RMS 識別起始故障點結果Fig.6 The DIT identification result with RMS model

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

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

圖8 Kurtosis 識別起始故障點結果Fig.8 The DIT identification result with kurtosis model

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

圖9 FC 識別起始故障點結果Fig.9 The DIT identification result with FC model

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

2.2 DFPF 方法預測軸承RUL 的分析對比

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

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

基于ACC 方法和γ確定表1 驗證數據集的DIT和實際失效時刻kT，采用最小二乘法擬合對應HI 指標值得到全局經驗模型參數a0，b0以及標準差σq，σr和σv。其中，初始權重系數取λ0=1 表示沒有觀測值輸入時，粒子的校正完全依賴所得全局經驗模型參數。DFPF 算法的主要結構參數如表2 所示。

表1 驗證數據集的DIT 和失效時刻Tab.1 DIT and actual degradation time of the testing data

表2 DFPF 算法主要結構參數Tab.2 The main structural parameters of DFPF algorithm

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

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

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

表3 特定時間段下軸承RUL 預測誤差結果Tab.3 Results of bearing RUL prediction error under a specific time period

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

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

此外，分別用均方根誤差Ermse、平均絕對百分比誤差Emape和平均絕對誤差Emae指標評價兩種算法的RUL 預測精度，具體計算公式為：

以上指標越小代表預測精度越高，具體結果對比如表4 所示。

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

當輸入觀測值比例為較少的5%時，軸承RUL預測結果對應輸出圖12（a）～（c）的左起第一個點，可直觀看出其接近實際值。其原因可能是預測結果受初始模型參數影響，而以上軸承的初始模型參數均由對應的全局指標數據擬合得出，準確度較高，因此即使觀測樣本較少，RUL 預測結果仍會接近真實壽命。此外，對于bearing2-1 和bearing1-3，占比為10%至30%對應的預測RUL 也接近真實值，進一步分析原因可能是該比例下的實際觀測值波動較小，整體較接近全局經驗模型。以上分析可知在輸入觀測值個數不足的條件下，觀測值和初始模型參數會影響軸承RUL 預測。

3 結論

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

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

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