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

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

0 引 言

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

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

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

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

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

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

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

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

1 研究方法

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

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

(1)

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

1.1 自助-最小二乘法

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

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

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

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

(2)

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

且有:

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

(3)

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

Yb=abI+cb

(4)

且有:

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

(5)

其中:i1～iq與數(shù)值1～10一一對應(yīng)。

其最小二乘解為:

(6)

(7)

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

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

(8)

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

(9)

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

1.2 最大熵原理

1.2.1 概率密度求取

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

(10)

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

最大熵的主要思想是在所有可行解中,滿足熵最大的解是最“無偏”的。最大熵方法能夠?qū)ξ粗母怕史植甲龀鲋饔^偏見為最小的最佳估計。

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

(11)

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

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

1.2.2 參數(shù)估計

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

(12)

對于雙側(cè)分位數(shù),有概率:

(13)

(14)

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

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

1.3 步驟分析

具體的步驟如下:

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

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

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

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

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

2 實驗及結(jié)果分析

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

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

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

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

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

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

2.1 振動性能序列X采集

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2.3 最大熵概率密度求取

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

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

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

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

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

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

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

表1 系數(shù)a和c的估計真值與區(qū)間

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

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

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

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

表2 時間序列XA的擬合系數(shù)

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

表3 時間序列XB的擬合系數(shù)

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

表4 時間序列XC的擬合系數(shù)

2.4 振動性能的動態(tài)預(yù)報

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

上述過程反復運行,便可實現(xiàn)各個時間段振動性能的動態(tài)預(yù)報。

時間序列XA第11～100步的預(yù)報結(jié)果如圖10所示。

圖10 時間序列XA各時間點振動性能的預(yù)報結(jié)果

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

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

時間序列XB第11～100步的預(yù)報結(jié)果如圖11所示。

圖11 時間序列XB各時間點振動性能的預(yù)報結(jié)果

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

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

時間序列XC第11～100步的預(yù)報結(jié)果如圖12所示。

圖12 時間序列XC各時間點振動性能的預(yù)報結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)束語

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

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

研究結(jié)果表明:

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

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

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

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