滾動軸承振動性能序列動態預報研究*

邵立東,常 振,陸水根,曹茂來

(1.杭州職業技術學院 吉利汽車學院,浙江 杭州 310018;2.杭州軸承試驗研究中心有限公司,浙江 杭州 310022;3.機械工業軸承產品質量檢測中心,浙江 杭州 310022)

0 引 言

滾動軸承是機電行業中應用最為廣泛的基礎件與精密件之一,如何確保其在各類復雜的環境中正常運轉,并準確高效地完成其規定功能,是軸承維護過程中的重點與難點。

由于滾動軸承壽命的離散性較大,定期更換或保養并不可行;而滾動軸承振動信號是監測其性能失效的重要指標,因此根據振動性能判斷分析軸承運轉狀況變得十分必要[1,2]。

從服役開始至失效結束,軸承退化史是一個漸變的過程,在軸承發生異常行為或故障之前,其振動信號可能早有預示。及時有效地對滾動軸承振動性能進行預報和診斷,可提醒機修人員在軸承損壞之前采取相應的維護或保養措施,進而避免機組失效、生產停滯,甚至人員傷亡等惡性事故的發生[3,4]。

然而,滾動軸承振動信號成分雖然具有明顯的隨機性及不確定性,基于已知概率密度函數的傳統統計理念的分析方法難以奏效,但其中確實又隱含有軸承健康狀況的潛在信息和確定性規律。

諸多學術界與工程界研究者對如何基于振動信號進行軸承狀態監控及預報做了研究,并已獲得了一些可行性方法。

鄧四二、張文平等人[5,6]搭建了滾動軸承的動力學數學模型,對不同表面特征參數(表面粗糙度、圓度、波紋度等)、工況參數,及諧波參數下滾動軸承的振動特性進行了理論分析。陳奕雅、GRASSO M等人[7-9]利用滾動軸承故障后的振動時間序列,對故障信號進行了特征提取,并根據二階循環統計量、經驗模式分解技術,實現了對軸承故障狀況的有效識別及診斷。程金山等人[10]根據不同的潤滑脂流變參數,測定出相應工況下的軸承振動值,研究了潤滑脂的膠體結構、彈性模量、表觀黏度等變化特性對軸承振動值的影響規律。李洪儒、SINGH S等人[11,12]針對軸承振動信號隨機強、抗干擾能力差、波動劇烈等問題,提出了二元多尺度熵和多函數融合的滾動軸承退化趨勢預測方法,可快速有效地進行軸承故障/失效預報。葉亮、SUN F等人[13-15]基于滾動軸承的振動信號提取出了軸承衰退性能指標,分別建立了數據驅動的可靠性與壽命預測模型,進而實現了滾動軸承可靠性與壽命的有效預測。

上述研究成果大都是針對軸承振動性能所進行的信號提取、故障診斷、狀態分析等研究,而對其振動性能在未來狀態下的演變預報的研究少之又少。

為了對軸承振動性能序列進行動態預報,筆者將自助法與最小二乘法進行有效融合,提出一種自助-最小二乘線性擬合的動態預報模型。首先,運用自助法將緊鄰且不斷更新的10個軸承振動數據進行10 000次仿真抽樣,形成10 000個樣本含量為10的振動序列;然后,將該振動序列用最小二乘法線性擬合,分別得到各自的線性擬合解a、c;其次,憑借最大熵原理對10 000組擬合解進行概率密度求取,在給定置信水平下得到相應的擬合真值與擬合區間;最后,根據未來時間變量實現滾動軸承振動性能真值與區間的動態預報。

1 研究方法

此處筆者設滾動軸承振動性能序列X為:

X=(x1,x2,…,xn,…,xN)

(1)

式中:xn—原始序列X的第n個振動數據;N—原始數據個數。

1.1 自助-最小二乘法

自助-最小二乘法是自助法與最小二乘法的融合。自助法可將一組數據等概率、可放回地抽樣多次,進而構成多組數據;最小二乘法可將一組數據通過最小化誤差平方和,尋找其函數的最佳匹配;二者的有效融合,可將一組數據進行多次擬合,進而獲得多組數據的最小二乘解。

具體實施時,取緊鄰的10個振動信號xi,xi+1,xi+2,…,xi+9,其中,i=1,2,…,N-9。以1,2,…,10為自變量,x為因變量,利用自助-最小二乘法進行線性擬合。

運用自助法,對緊鄰的10個滾動軸承原始振動信號等概率、可放回地隨機抽取1個數,共抽取q=10次,得到一個自助樣本Y1,其有q=10個數據。按此方法重復執行B次,得到B個樣本,可表示為:

YBootstrap=(Y1,Y2,…,Yb,…,YB)

(2)

式中:Yb—軸承振動信號的第b個自助樣本,b=1,2,…,B;B—總的自助再抽樣次數,也是自助樣本的個數。

且有:

Yb=(y1,y2,…,yl,…,yq)

(3)

利用最小二乘法對Yb進行線性擬合:

Yb=abI+cb

(4)

且有:

I=(1,2,3,…,10)=(i1,i2,…,il,…,iq)

(5)

其中:i1～iq與數值1～10一一對應。

其最小二乘解為:

(6)

(7)

由于線性擬合共進行B次,則獲得的B個最小二乘解，可表示為:

a=(a1,a2,…,ab,…,aB)

(8)

c=(c1,c2,…,cb,…,cB)

(9)

為了獲得軸承振動信號自助樣本的最小二乘解的估計真值a0、c0,以及上、下區間[aL,aU]和[cL,cU],筆者運用最大熵原理,分別對最小二乘解a和c進行概率度求取。

1.2 最大熵原理

1.2.1 概率密度求取

此處筆者以最小二乘解系數a為例,進行概率密度求取,將式(8)的B個系數a連續化,定義最大熵的表達式為:

(10)

式中:p(a)—連續化后的數據序列a的概率密度函數。

最大熵的主要思想是在所有可行解中,滿足熵最大的解是最“無偏”的。最大熵方法能夠對未知的概率分布做出主觀偏見為最小的最佳估計。

通過調整p(a)可以使熵達到最大值,拉格朗日乘子法的解可表示為[16]:

(11)

式中:θ0,θ1,…,θβ—拉格朗日乘子;a—最小二乘解系數a的隨機變量。

根據該概率密度函數p(a)可實現該數據序列的估計真值a0與上、下區間[aL,aU]的預報。

1.2.2 參數估計

由隨機變量a的概率密度函數p(a),可得序列a的估計真值a0為:

(12)

對于雙側分位數,有概率:

(13)

(14)

式中:aU,aL—最小二乘解系數a的上界值和下界值;δ—實數,δ∈(0,1);[aL,aU]—δ水平下的置信區間[17]。

同理,可得到最小二乘解系數c的估計真值c0與上、下界[cL,cU]。所以,在自助-最小二乘法線性擬合時,根據最大熵原理,可實現最小二乘解的最優估計及上、下區間預報。將最小二乘解代入擬合方程(4),即可獲得滾動軸承振動信號的最優估計值x0和上、下區間[xL,xU]。

1.3 步驟分析

具體的步驟如下:

(1)采用滾動軸承振動信號傳感器,采集滾動軸承振動性能序列X;

(2)緊鄰的10個軸承原始振動數據運用自助-最小二乘法線性擬合,獲得樣本含量為B的最小二乘解向量序列a和c;

(3)利用最大熵原理對序列a和c建立各自的概率密度函數,在給定置信水平下,求取最小二乘解的最優估計值a0、c0和上下限[aL,aU]、[cL,cU];

(4)再將最優估計值a0、c0和上下限[aL,aU]、[cL,cU]代入式(4)的線性擬合方程,設定自變量為11～20,進而預報出下一時間段(第11～20)的滾動軸承振動性能的真值x0及上、下區間[xL,xU];

(5)不斷更新緊鄰的10個軸承原始振動數據,重復步驟2、3、4,從而實現滾動軸承振動性能真值與上、下區間的動態預報。

2 實驗及結果分析

由于這是一個超精密滾動軸承的壽命強化試驗,此處筆者采用型號為ABLT-1A的試驗機。該試驗機主要由試驗頭、試驗頭座、傳動系統、加載系統、潤滑系統和計算機控制系統組成。

軸承壽命強化試驗臺架如圖1所示。

圖1 軸承壽命強化試驗臺架

試驗樣品為超精密滾動軸承H7008C,且該類軸承有一個過渡等級,尚可達到國標P2級的精度要求。試驗在電機轉速為4 950 r/min,室溫為26 ℃,濕度為53%的環境條件下進行,所施加徑向載荷為5 kN,軸向載荷2 kN。

試驗頭原理圖如圖2所示。

圖2 軸承壽命試驗頭原理圖

2.1 振動性能序列X采集

試驗時間及軸承振動信息由計算機控制系統自動累積顯示,每間隔10 min由傳感器采集一次振動加速度信號,單位為ms-2。

筆者分別采集該型號軸承初期、中期、臨尾3個服役時間狀態下,振動加速度數據作為試驗研究的案例。

初期振動時間序列記為XA,如圖3所示。

圖3 軸承初期振動性能時間序列XA

中期振動時間序列XB,如圖4所示。

圖4 軸承中期振動性能時間序列XB

臨尾階段的振動數據序列為XC,如圖5所示。

圖5 軸承臨尾振動性能時間序列XC

由圖(3～5)的軸承振動原始數據不難看出:振動信號表現出明顯的隨機性及不確定性,難以用準確的公式或相應的概率分布函數對其進行描述;且不同服役階段振動值的大小、波動區間、變化趨勢各不相同,這對軸承振動信號的預報十分不利。

為有效地克服以上問題,筆者提出一種軸承振動信號預報模型,然后進行具體的數據分析處理。

2.2 自助-最小二乘法線性擬合

以初期振動性能時間序列XA為例,筆者將樣本含量為100的XA序列等分成10組,即每組10個樣本數據(x1-x10,x11-x20,…,x91-x100)。首先,以x1-x10為訓練值進行有效預報,預報步長為10步:x1-x10內的10個數據運用自助法等概率、可放回地抽樣10次,并連續重復10 000次,可構成樣本含量為10 000的自助序列Ybootstrap=(Y1,Y2,…,Yb,…,Y10 000),其中,Yb的樣本含量為10;再利用最小二乘法,對Yb內的10個抽樣數據進行線性擬合,可分別得到擬合參數ab和cb;最后構成樣本含量為10 000的最小二乘解序列a=(a1,a2,…,ab,…,a10 000)和c=(c1,c2,…,cb,…,c10 000)。

最小二乘解序列a結果如圖6所示。

圖6 時間序列XA中x1-x10的自助-最小二乘解序列a

最小二乘解序列c結果如圖7所示。

圖7 時間序列XA中x1-x10的自助-最小二乘解序列c

2.3 最大熵概率密度求取

采用自助-最小二乘法進行線性擬合的方式,可由圖(6～7)得到時間序列XA中x1-x10這10個緊鄰的振動數據的解集a和c;再將B個解集連續化,得到其對應的最大熵表達式,即式(10);根據拉格朗日乘子法,可得到該軸承當前序列段的擬合系數a和c的概率密度函數,即式(11)。

在時間段XA中,軸承x1-x10的擬合系數a的概率密度函數,如圖8所示。

圖8 時間序列XA中x1-x10的擬合系數a的概率密度圖

在時間段XA中,軸承x1-x10的擬合系數c的概率密度函數,如圖9所示。

圖9 時間序列XA中x1-x10的擬合系數c的概率密度圖

接下來,需要對參數a和c進行估計。根據圖(8,9)所得的自助-最小二乘解系數a和c的概率密度函數,結合式(12),便可得到其各自的估計真值a0和c0;然后給定顯著水平δ=0.1(即置信水平為90%),根據式(13,14)得到其對應的區間上下限[aL,aU]和[cL,cU]。

所以,時間序列XA中x1-x10的擬合系數的預報結果,即系數a和c的估計真值與區間如表1所示。

表1 系數a和c的估計真值與區間

表1中的結果便是XA中,采用自助-最小二乘法進行線性擬合x1-x10這10個緊鄰的振動數據的解,再分別將其代入線性擬合式(4),便可得到x1-x10這10個緊鄰的軸承振動性能數據的擬合方程;令自變量等于11～20,便可獲得第11～20個時間段(即未來時間段)振動性能的預報真值及上、下限。

在進行第21～30個時間段的振動性能預報時,采取舊數據舍棄、新數據更替的原則,將原來的x1-x10原始數據舍棄,添加10個新的原始振動性能值x11-x20;同時,將x11→新x1,x12→新x2,x13→新x3,…,x20→新x10,即將x11-x20這10個緊鄰的原始振動數據轉變為新的x1-x1010個振動性能數據;然后,重復以上步驟,對其進行自助-最小二乘法線性擬合;同樣,令自變量等于11～20,便可獲得第21～30個時間段振動性能的預報真值及上下限。

同理,依次類推,可分別得到第31～40,41～50,51～60,…,等時間段的線性擬合結果。

在各時間段內,時間序列XA的自助-最小二乘法擬合系數,如表2所示。

表2 時間序列XA的擬合系數

在各時間段內,時間序列XB的自助-最小二乘法擬合系數,如表3所示。

表3 時間序列XB的擬合系數

在各時間段內,時間序列XC的自助-最小二乘法擬合系數,如表4所示。

表4 時間序列XC的擬合系數

2.4 振動性能的動態預報

根據表(2～4),運用自助-最小二乘法線性擬合,可分別得到軸承各服役時間段,10個緊鄰的振動性能的系數真值及上、下限;將其代入線性擬合方程,并設置自變量為11～20,可分別得到時間序列XA、XB、XC第11～20,第21～30,第31～40,…,等時間段振動性能的預報真值x0及上下限[xL,xU],即x0=a0×l+c0,xL=aL×l+cL,xU=aU×l+cU,其中:l=11,12,13,…,20。

上述過程反復運行,便可實現各個時間段振動性能的動態預報。

時間序列XA第11～100步的預報結果如圖10所示。

圖10 時間序列XA各時間點振動性能的預報結果

由圖10可得:時間序列XA各時間段振動性能的預報真值與實際值相差極小,第43次(即第53時間點)的預報真值與實際值相差最大,但僅為0.129 ms-2;振動性能的預報區間可將其實際值全部包羅,且上、下區間差值小,預報精度高。

以上結果表明,預報模型應用于軸承初期服役階段的振動性能預報時是準確可行性的。

時間序列XB第11～100步的預報結果如圖11所示。

圖11 時間序列XB各時間點振動性能的預報結果

由圖11可得:時間序列XB各時間段振動性能的預報真值有逐漸上升的總趨勢,與實際值的變化趨勢保持良好的一致性;第28次(即第38時間點)的預報真值與實際值相差最大,但僅為0.188 ms-2;同樣,振動性能的預報區間可將其實際值全部包羅,且上下區間差值小,預報精度高。

以上結果表明,預報模型應用于軸承中期服役階段的振動性能預報時,同樣是準確可行性的。

時間序列XC第11～100步的預報結果如圖12所示。

圖12 時間序列XC各時間點振動性能的預報結果

由圖12可得:時間序列XC各時間段振動性能的預報真值在實際值的均值附近波動,且波動范圍極小;第85次(即第95時間點)的預報真值與實際值相差最大,但僅為1.379 ms-2;同樣,振動性能的預報區間可將其實際值全部包羅,且上下區間差值小,預報精度高。

以上結果可以說明,預報模型應用于軸承臨尾服役階段的振動性能預報時,也是準確可行的。

為直觀看出可靠度預報值與實際值的差異,筆者分別計算出時間序列XA、XB、XC預報真值與實際值之間的相對誤差,其結果如圖13所示。

圖13 可靠度預報值與實際值之間的相對誤差

從圖13可以看出:時間序列XA的最大相對誤差出現在第37次(即第47時間點),但僅為14.73%;最小相對誤差出現在第57次(即第67時間段),為0.29%;

時間序列XB的最大相對誤差出現在第28次(第38時間段),但僅為4.57%;最小相對誤差出現在第52次(第62時間段),為0.09%;

時間序列XC的最大相對誤差出現在第85次(即第95時間段),但僅為9.91%;最小相對誤差出現在第54次(第64時間段),為0.002%。

所以,振動性能預報值與實際值的相對誤差較小,最大不超過15%,再次說明預測結果是十分真實可靠的,并可較好地應用于工程實際。

總體來看,初期服役時間序列XA的預報值與實際值差值最小,但相對誤差最大,這是由于軸承服役初期振動小,即誤差求取過程中分母基數低;中期服役時間序列XB的預報值與實際值差值較小,但相對誤差最小,這是由于軸承服役中期振動較為穩定,預報值與實際值的上下波動小;

臨尾服役時間序列XC的預報值與實際值差值最大,相對誤差較小,這是由于軸承服役臨尾振動較為劇烈,預報值與實際值的上下波動大,但誤差求取過程中分母基數高。

基于時間序列XA、XB、XC的3個不同服役階段試驗案例的自助-最小二乘法線性擬合模型,可將緊鄰且不斷更新的10個振動性能值進行有效擬合;然后,根據最大熵原理將多個最小二乘解進行概率密度求取,在給定90%的置信水平下獲得擬合解的最優估計和上下區間;最后,給定自變量并代入線性擬合方程,便可獲得下一時間段振動信號的預報真值與區間上下限。

基于該模型得到的預報結果與實際值相對誤差不超過15%,滿足工程實際的預報要求;預報區間可將其實際值全部包羅,且上下區間差值小,預報精度高。

3個案例的振動性能預報值與實際值均保持良好的一致性,差值小、誤差低、精度高,說明筆者所提出的模型具有良好的準確性及可靠性。此外,該模型不需考慮數據分布的任何信息,只針對現有數據做出最真實、客觀的判斷。該模型實現了對滾動軸承自我狀況的在線監測。

3 結束語

為了對軸承振動性能序列進行動態預報,筆者將自助法與最小二乘法進行有效融合,提出了一種基于自助-最小二乘線性擬合的軸承振動性能序列動態預報模型。

首先,筆者對自助法與最小二乘法做了有效融合,對軸承振動時間序列進行了多次擬合;運用最大熵原理描述出擬合結果的概率密度函數,給定了相應的置信水平和時間變量,并不斷地更替新舊數據,實現了滾動軸承振動性能的真值與區間的動態預報;通過3組實驗數據驗證了所提模型的準確性。

研究結果表明:

(1)自助-最小二乘線性擬合方法,可將振動時間序列多次擬合,有效反映出原始數據變化趨勢的多個側面信息;

(2)最大熵原理可準確描述出擬合結果的概率密度函數,快速求取擬合解的最優測度及上下區間;

(3)所提模型還可實現預報結果的自我驗證,時間序列XA振動預報值與實際值的最大相對誤差為14.73%;時間序列XB的最大相對誤差僅為4.57%;時間序列XC的最大相對誤差僅為9.91%;該預報結果滿足工程實際的一般要求。

該預報模型可為后續滾動軸承自我健康檢測以及在線故障診斷的研究提供理論根據;并可在工程實際中及時地發現問題,提前做好軸承早期失效的預防與監測工作。