基于振動時間序列的汽車輪轂軸承壽命威布爾估計

卜德天，陳龍，劉紅彬，呂文靜，張弘毅

(1.河南科技大學 機電工程學院，河南 洛陽 471003；2.洛陽LYC汽車軸承科技有限公司，河南 洛陽 471003；3.洛陽軸承研究所有限公司, 河南 洛陽 471039)

汽車輪轂軸承[1](Hub Bearing)是輪轂總成的重要組成部分，其疲勞壽命與可靠性計算的確定對汽車行駛安全性有重大影響。Lundberg和Palmgren基于赫茲接觸理論和威布爾材料強度統計理論建立了滾動軸承疲勞壽命估計方法，被廣泛接受并成為軸承行業的計算標準[2]。

Gustafsson[3]于1962年首次提出采用加速度傳感器采集軸承振動信號分析的方法，振動監測隨后廣泛應用于回轉機械的運行狀態監測。相比于基于聲信號[4]、靜電信號[5]和定子電流[6]的監測方法，振動監測不僅能夠形象、直觀地描述軸承性能，還能夠從振動信號中提取性能退化過程的特征信息[7-8]。振動監測獲取的大量狀態演變數據也為基于振動的軸承性能演變以及剩余壽命預測提供了充足的數據源。

依據振動數據測量時間的先后關系，一些學者采用振動時間序列方法開展機械設備的運行狀態預測：文獻[9]以神經網絡為基礎構建振動頻率時間序列預測方法，提高了旋轉機械設備運行狀態預測精度；文獻[10]開展潮流能水平軸水輪機振動位移時間序列的自相似性研究，為采用自相似理論研究潮流能水輪機的動態特性奠定了數學基礎，所揭示的時間序列整體與局部之間的關系及其自相似性為潮流能水輪機的故障診斷提供了技術支撐；文獻[11]進行變壓器油箱表面振動信號時間序列的混沌特性分析，證明變壓器表面振動信號的時間序列具有混沌特性，為變壓器振動信號進一步的處理及應用提供了參考；文獻[12]基于時間序列分析機車振動信號，建模和預測結果表明時間序列分析技術能夠快速有效地處理振動信號的預測問題。

由于數據量充足，統計分析方法也被廣泛應用于分析滾動軸承的振動數據：文獻[13]在研究三參數威布爾分布參數的分位數估計、擬估計和改良極大似然估計時采用隨機模擬法對其優良性進行了探究，并給出了可靠度的置信區間；文獻[14]在研究葉片發動機的故障實例時將最小二乘法和平均秩次法結合對威布爾參數進行計算；文獻[15]基于Nelson定律，采用等累積退化量和等累積失效概率的模型得到了不同應力下的等效試驗時間；文獻[16]在研究鐵路設備關鍵零部件的可靠性分析中提出了一種新的相關系數法，對收集數據進行三參數威布爾分布參數估計并建立了鐵路設備關鍵零部件可靠性分析模型；文獻[17]在對導彈武器系統的性能檢測時研究了小樣本情況下威布爾參數的計算；文獻[18]以鉸銷式碳纖維增強基復合材料環帶為研究對象，基于SWT疲勞壽命預測模型推導了服從威布爾分布的概率疲勞壽命預測方法并證明了此預測模型的正確性。

本文在研究汽車輪轂軸承疲勞壽命試驗振動數據時，將振動時間序列劃分為多個子序列，對各子序列的數據開展統計分析，以研究不同序列中威布爾分布參數的變化，總結形狀參數與尺度參數的變化規律，實現軸承疲勞失效的提前預報。

1 試驗

試驗軸承為雙列角接觸汽車輪轂軸承，型號為DAC2F11410561，內徑為34 mm、外徑為83 mm，球材料為GCr15，球直徑為12.7 mm，球組節圓直徑為57 mm，接觸角為35°。

1.1 試驗設備

采用SYJ-LG-NJ輪轂軸承耐久試驗機進行試驗，如圖1所示，試驗機由試驗主體部件、床身部件、傳動系統、加載系統、潤滑系統、電氣控制系統、計算機軟件系統等組成。試驗機采用計算機自動控制，可按預設試驗載荷、轉速和運行時間等主要參數自動運行，在試驗過程中可實時獲取軸承每個時間點的轉速、電流、振動加速度、載荷、油溫、軸溫等數據，試驗機參數見表1。

1—軸向油缸；2—徑向油缸；3—試驗軸承；4—主軸；5—電動機。圖1 試驗機結構示意圖Fig.1 Structure diagram of tester

表1 試驗機參數Tab.1 Parameters of tester

1.2 試驗工況

試驗為定速、定載加速壽命試驗，徑向載荷為8.8 kN，軸向載荷為5.4 kN，轉速為400 r/min。試驗結束后軸承單元狀態和軸承內圈溝道狀態如圖2所示。

(a) 軸承單元 (b) 軸承內圈溝道圖2 試驗結束后軸承單元和軸承內圈溝道的狀態Fig.2 Status of bearing element and bearing inner ring raceway after test

1.3 試驗數據

按照現有疲勞壽命試驗經驗，試驗過程中被測軸承振動加速度值達到初始值2倍時判定軸承發生疲勞失效，試驗終止。本試驗歷時8 320 min，試驗過程中每5 s提取一個振動加速度值，總計提取99 840個數據。

時間序列是指某一變量或指標的數值或觀測值按時間順序以相同時間間隔排列的一組數據。為表征軸承振動的演變狀態，將振動數據按測量時間先后等間隔排列，構成的振動時間序列如圖3所示。

圖3 DAC2F1141056軸承振動時間序列Fig.3 Vibration time series of DAC2F1141056 bearing

2 振動時間序列分析

振動時間序列包含的大量數據為基于統計理論的分析研究提供了基礎條件，但一般統計分析方法難以描述時間先后特征。為表征被測軸承試驗中的狀態演變，將整體時間序列劃分為具有時間先后關系的子序列，分析各子序列威布爾分布的形狀參數與尺度參數的變化狀態，與各子序列對應的振動時間序列比較，實現被測軸承疲勞壽命的威布爾估計。

2.1 分析方法

假設所記錄數據(軸承振動值)的總量為N，則軸承振動時間序列Q可描述為

Q=[q(1),q(2),…,q(s),…,q(N)];

s=1,2,…,N，

(1)

式中：q為試驗記錄的振動值；s為數據序號。

所有試驗數據構成總體樣本，依據總體樣本的時間先后關系將振動時間序列分為J組，每個子序列包含的樣本數量為N/J，子序列可表示為

Qj=[q(j+1),q(j+2),…,q(j+s),…,

q(j+N/J)]。

(2)

不同區間組的統計參數表征被測軸承在統計區間的性能狀態，對于第j個子序列，威布爾分布的概率密度函數為

(3)

式中：k，c分別為威布爾分布的形狀參數和尺度參數。

第j個子序列的威布爾分布的累積密度函數為

(4)

可改寫為

(5)

對等式兩邊同時取2次對數可得

(6)

y=kx-klgc。

(7)

以上分析過程中可獲取不同子序列的形狀參數與尺度參數，將不同子序列的統計參數按時間先后關系排序，則可通過參數變化觀測軸承狀態演變過程的特征。

2.2 分析結果

將99 840個數據分為20組，即J=20，如圖3所示。前19組中每組包含5 000個振動數據樣本，第20組為剩余的4 840個數據。分別計算20組數據的形狀參數和尺度參數，結果見表2：前19個子序列的形狀參數和尺度參數在一定范圍內波動，無顯著規律，說明前19個子序列內軸承狀態無明顯變化；第20個子序列內軸承狀態則發生了顯著變化。

表2 不同子序列的形狀參數和尺度參數Tab.2 Shape parameters and scale parameters of different sub-series

依據(3)式、(4)式、(7)式分別繪制各子序列的概率密度函數圖、累積分布函數圖、直線方程圖，結果如圖4所示：前19個子序列對應的概率密度函數與累積分布函數分布均較為集中，無顯著變化；第20個子序列對應函數的橫坐標跨度明顯增加，峰值右移，表征第20個子序列內軸承振動加速度增加，與試驗中振動時間序列吻合；前19個子序列的斜率與截距也無明顯變化，而第20個子序列的斜率顯著降低，截距顯著增加。

以上多個方面的結果均表明：威布爾分布的形狀參數和尺度參數能描述軸承狀態演變過程，不同區間數據的統計規律能表征軸承狀態演變的內在屬性。

(a) 概率密度函數 (b) 累積分布函數 (c) 直線方程圖4 不同子序列的威布爾分布函數Fig.4 Weibull distribution functions of different sub-series

3 疲勞預測研究

以上研究采用的是整個疲勞壽命試驗的完整壽命周期數據，包含了最終失效數據，而實際工程應用中需要提前預報軸承的疲勞失效，為提升本研究的工程應用價值，選取發生疲勞失效前一段時間內的數據進一步分析。

依據原始試驗數據，第98 407個振動加速度為初始振動加速度的2倍，為達到預測目的，截取第88 288～98 287區間內的10 000個數據開展分析，選取的數據區間為由失效值往前120個數據(換算時間單位為10 min)作為研究區間的最后一個數據，截取數據在整個振動時間序列的位置標注(圖3)。將所截取數據區間內的10 000個數據劃分為10個子序列，每個子序列包含1 000個數據樣本。依據上述處理方法得到截取時間段內不同子序列的形狀參數和尺度參數，結果見表3，各子序列的威布爾分布函數如圖5所示。

表3 截取時間段內的不同子序列的形狀參數和尺度參數Tab.3 Shape parameters and scale parameters of different sub-series in interception time period

(a) 概率密度函數 (b) 累積分布函數 (c) 直線方程圖5 截取時間段內不同子序列的威布爾分布函數Fig.5 Weibull distribution functions of different sub-series in interception time period

分析可知：前8個子序列形狀參數的平均值為24.840 0且數值均在25左右波動，最大值、最小值分別為27.324 7和21.094 5，與均值差異的幅度分別為9.9%和15.1%；前8個子序列尺度參數的平均值為0.047 9且數值均在0.048左右波動，最大值、最小值分別為0.048 2和0.047 3，與均值差異的幅度分別為0.6%和1.3%；第9個子序列的形狀參數迅速衰減到15.380 6，與前8個子序列均值的變化幅度為38.1%，尺度參數也迅速變化為0.050 2，與前8個子序列均值的變化幅度為4.8%；第10個子序列的形狀參數進一步衰減到8.387 4，與前8個子序列均值的變化幅度達66.3%，尺度參數也進一步變化為0.055 9，與前8個子序列均值的變化幅度為16.7%；圖5中不同子序列威布爾分布函數的對比也清晰表明了第9，10個子序列的顯著變化。

第10個子序列最后1個數據換算為時間單位是失效前10 min，假定第10個子序列形狀參數或尺度參數的變化幅度達到判定標準，即提前10 min實現被測軸承疲勞失效的預報；第10個子序列包含1 000個數據，換算為時間單位是83.3 min，假定第9個子序列形狀參數或尺度參數變化幅度達到判定標準，則提前93.3 min實現被測軸承疲勞失效的預報。

4 驗證

為驗證以上分析方法，選取相同型號、相同試驗條件下的另1套DAC2F11410561軸承和相同試驗機上試驗的另1套不同型號的汽車輪轂軸承的試驗數據開展對比驗證。

4.1 同型號軸承的驗證

另1套DAC2F11410561軸承的疲勞壽命試驗記錄了136 260個振動加速度值，其振動時間序列如圖6所示。采用前述分析方法，以振動加速度達到初始振動加速度的2倍為節點，提前10 min作為分析區間的最后1個數據，選取該軸承100 265～110 264區間內的10 000個數據依次分為10組進行分析，用于驗證的數據在整個時間序列的位置如圖6所示。

10個子序列的形狀參數和尺度參數的變化趨勢如圖7所示：2種參數的變化規律與表3具有較好的一致性，即前8個子序列的形狀參數和尺度參數均在一定范圍內波動，第9，10個子序列則發生顯著波動，變化趨勢和規律與前述分析結果一致。

圖6 驗證用DAC2F1141056軸承振動時間序列Fig.6 Vibration time series of DAC2F1141056 bearing used for verification

(a) 形狀參數 (b) 尺度參數圖7 驗證用DAC2F1141056軸承截取時間段內不同子序列的形狀參數和尺度參數Fig.7 Shape parameters and scale parameters of different sub-series in interception time period of DAC2F1141056 bearing used for verification

4.2 其他型號軸承的驗證

為驗證以上規律是否適用于所有汽車輪轂軸承，選取新開發某型汽車輪轂軸承疲勞壽命試驗中的振動數據，構建的振動時間序列如圖8所示，包含139 318個軸承振動加速度數據，與前述截取規則相同，選取123 148～133 147區間內的10 000個數據，依次分為10個子序列進行分析，計算驗證數據10個子序列的形狀參數和尺度參數，其變化趨勢如圖9所示：2種參數的變化規律與表3、圖7均具有較好的一致性，即前8個子序列的形狀參數和尺度參數在一定范圍內波動，第9，10個子序列則發生顯著波動，說明本分析方法對汽車輪轂軸承振動時間序列的分析具有通用性。

圖8 驗證用其他型號汽車輪轂軸承的振動時間序列Fig.8 Vibration time series of automobile hub bearings of other models used for verification

(a) 形狀參數 (b) 尺度參數圖9 驗證用其他型號汽車輪轂軸承截取時間段內不同子序列的形狀參數和尺度參數Fig.9 Shape parameters and scale parameters of different sub-series in interception time period of automobile hub bearings of other models used for verification

5 結束語

振動時間序列演變的統計規律對于滾動軸承在線監測裝置的應用有良好的工程價值，相較于現有滾動軸承剩余壽命預測方法，本文提出的方法容易進行工程化，且提前預報的時間長度長于現有研究方法。但是，本研究尚在進行中，后續應在以下方面作進一步研究：

1)擴展振動時間序列數據源。從汽車輪轂軸承擴展至其他類型或其他應用領域的軸承，以進一步驗證本方法的通用性。

2)判定標準的界定。通過其他理論建立關聯性，以確定形狀參數和尺度參數變化程度的量化判定標準。

3)實現滾動預測分析。本文為了方便定性分析，采用自下而上的分析方法，即以失效節點倒推截取數據；實際應用中可采用程序化實現滾動預測，在獲取最初1 000個數據后即可計算表征軸承起始狀態的形狀參數和尺度參數，當測取到1 001個數據時去除第1個數據，采用2～1001個數據計算形狀參數和尺度參數表征更新的軸承狀態，依次循環，得到形狀參數與尺度參數的連續演變規律，最終實現工程化。