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

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

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

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

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

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

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

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

1 試驗

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

1.1 試驗設(shè)備

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

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

表1 試驗機參數(shù)Tab.1 Parameters of tester

1.2 試驗工況

試驗為定速、定載加速壽命試驗，徑向載荷為8.8 kN，軸向載荷為5.4 kN，轉(zhuǎn)速為400 r/min。試驗結(jié)束后軸承單元狀態(tài)和軸承內(nèi)圈溝道狀態(tài)如圖2所示。

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

1.3 試驗數(shù)據(jù)

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

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

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

2 振動時間序列分析

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

2.1 分析方法

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

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

s=1,2,…,N，

(1)

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

所有試驗數(shù)據(jù)構(gòu)成總體樣本，依據(jù)總體樣本的時間先后關(guān)系將振動時間序列分為J組，每個子序列包含的樣本數(shù)量為N/J，子序列可表示為

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

q(j+N/J)]。

(2)

不同區(qū)間組的統(tǒng)計參數(shù)表征被測軸承在統(tǒng)計區(qū)間的性能狀態(tài)，對于第j個子序列，威布爾分布的概率密度函數(shù)為

(3)

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

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

(4)

可改寫為

(5)

對等式兩邊同時取2次對數(shù)可得

(6)

y=kx-klgc。

(7)

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

2.2 分析結(jié)果

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

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

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

以上多個方面的結(jié)果均表明：威布爾分布的形狀參數(shù)和尺度參數(shù)能描述軸承狀態(tài)演變過程，不同區(qū)間數(shù)據(jù)的統(tǒng)計規(guī)律能表征軸承狀態(tài)演變的內(nèi)在屬性。

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

3 疲勞預測研究

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

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

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

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

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

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

4 驗證

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

4.1 同型號軸承的驗證

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

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

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

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

4.2 其他型號軸承的驗證

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

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

(a) 形狀參數(shù) (b) 尺度參數(shù)圖9 驗證用其他型號汽車輪轂軸承截取時間段內(nèi)不同子序列的形狀參數(shù)和尺度參數(shù)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 結(jié)束語

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

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

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

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