時保吉，栗永非

(新鄉(xiāng)職業(yè)技術(shù)學(xué)院，河南 新鄉(xiāng) 453006)

隨著科學(xué)技術(shù)的快速發(fā)展，為確保系統(tǒng)安全可靠運行，現(xiàn)代機(jī)械對滾動軸承性能提出了更加嚴(yán)格的要求。其旋轉(zhuǎn)精度、振動、摩擦力矩波動性以及壽命等直接影響主機(jī)的可靠性和使用壽命[1-2]。對滾動軸承摩擦力矩變異過程的實時評估與預(yù)報，可以為滾動軸承摩擦磨損壽命和可靠性分析提供依據(jù)。滾動軸承摩擦力矩的時間序列屬于乏信息過程[3-4]。經(jīng)典統(tǒng)計理論依賴于已知的概率分布與大樣本數(shù)據(jù)，難以有效地解決滾動軸承摩擦力矩不確定性的動態(tài)預(yù)報問題。

目前，關(guān)于滾動軸承的摩擦力矩的研究已取得一些成果。文獻(xiàn)[5]提出了最大熵方法，對HTKB軸承摩擦力矩分布特征參數(shù)進(jìn)行了靜態(tài)評估與預(yù)測，并由試驗數(shù)據(jù)驗證了模型的正確性；文獻(xiàn)[6]以模糊集合理論為基礎(chǔ)，從小樣本數(shù)據(jù)入手，建立了航天軸承摩擦力矩參數(shù)的經(jīng)驗概率密度函數(shù)及其模糊預(yù)報模型，并對摩擦力矩的均值上界和最大值上界進(jìn)行預(yù)報，試驗研究表明，預(yù)報結(jié)果可以滿足航天工程的要求；文獻(xiàn)[7]提出了理論模型和試驗方法來表征改進(jìn)的推力球軸承在混合和全油膜潤滑條件下的摩擦力矩，并研究了潤滑劑黏度對摩擦力矩的影響；文獻(xiàn)[8]提出了轉(zhuǎn)盤軸承啟動力矩變異系數(shù)的自助最大熵評估方法，但以上成果尚未涉及摩擦力矩不確定性的動態(tài)預(yù)報。

下文融合了GM(1,1)模型完善的預(yù)報機(jī)制和自助法的模擬再抽樣能力的優(yōu)點，在乏信息條件下，提出了滾動軸承摩擦力矩不確定性的動態(tài)灰自助GBM(1,1)預(yù)報方法，通過摩擦力矩不確定性試驗，證實了該方法的正確性與有效性。

1 數(shù)學(xué)模型

假設(shè)所研究的滾動軸承摩擦力矩為隨機(jī)變量x。在滾動軸承試驗期間，對其摩擦力矩進(jìn)行定期采樣分析，假設(shè)獲得了該性能的R個時間單元的數(shù)據(jù)。令Xr為第r個時間單元采集的數(shù)據(jù)，并構(gòu)成第r個時間序列Xr

Xr={xr(n)};n=1,2,…,r=1,2,…,R,

(1)

式中：xr(n)為Xr中的第n個原始數(shù)據(jù)；n為數(shù)據(jù)序號，即時刻。

借助與時刻n緊鄰時刻的前m個數(shù)據(jù)(包括時刻n的數(shù)據(jù))評估時刻n的不確定度。在時刻n只考慮m個數(shù)據(jù)的原因是：在動態(tài)測量中，不確定度隨時間而變化，m越大，包含的舊信息越多，估計的不確定度的誤差越大?；疑到y(tǒng)理論的灰預(yù)測模型GM(1,1)和自助法要求m的最小值為4。

在時刻n，從Xr中抽取m個數(shù)據(jù)，形成一個小樣本數(shù)據(jù)子序列Xrm

Xrm={xrm(u)+qr};

u=n-m+1,n-m+2,…,n，n≥m，

(2)

式中：u為時間變量；qr為常數(shù)，在灰預(yù)測GM(1,1)中，若xrm(u)<0，qr應(yīng)滿足xrm(u)+qr≥0；若xrm(u)≥0，則取qr=0。

在時刻n，從Xrm中等概率可放回地隨機(jī)抽取一個數(shù)據(jù)，抽取m次，得到一個大小為m的自助樣本[9]。連續(xù)重復(fù)B步，得到B個自助再抽樣樣本，用向量表示為

YrBootstrap=(Yr1,Yr2,…,Yrb,…,YrB)，

(3)

Yrb={yrb(u)};

u=n-m+1,n-m+2,…,n，

n≥m，b=1,2,…,B，

式中：Yrb為第r個時間單元中的第b個自助樣本：yrb(u)為Yrb中的第u個自助再抽樣數(shù)據(jù)。

根據(jù)灰預(yù)測模型GM(1,1)[10]，設(shè)Yrb的一次累加生成序列向量為

(4)

u=n-m+1,n-m+2,…,n，

n≥m，b=1,2,…,B。

灰生成模型為

(5)

式中：cr1和cr2是待估系數(shù)(cr1≠0)。

用增量代替微分，即

xrb(u+1)-xrb(u)=yrb(u+1),

(6)

式中：Δu取單位時間間隔1，設(shè)均值生成序列向量為

Zrb={zrb(u)}={0.5xrb(u)+0.5xrb(u-1)};

(7)

u=n-m+2,n-m+3,…,n,

n≥m,b=1,2,…,B。

初始條件xrb(n-m+1)=yrb(n-m+1) 下， (5) 式的最小二乘解為

(8)

D=(-Zrb,I)Τ，I=(1,1,…,1) 。

時刻w=n+1的預(yù)測值為

(9)

因此，可以獲得第r個時間單元在w=n+1時刻的B個數(shù)據(jù)，構(gòu)成如下序列向量

(10)

用統(tǒng)計理論的直方圖方法，建立第r個時間單元在時刻w的概率密度函數(shù)

frw=frw(xrm) ，

(11)

式中：xrm是描述第r個時間單元中被測數(shù)據(jù)xrm(w)的隨機(jī)變量。

假設(shè)顯著性水平為α，則置信水平為

P=(1-α)×100%，

(12)

且在時刻n，在置信水平P下，xrm的估計區(qū)間為

[XrL,XrU]=[XrL(w),XrU(w)]=

[Xrα/2,Xr1-α/2],

(13)

式中：Xrα/2為xrm對應(yīng)于概率α/2的值；Xr1-α/2為xrm對應(yīng)于概率1-α/2的值；XrU，XrL分別為估計區(qū)間的上、下邊界。

定義n時刻的擴(kuò)展不確定度，即瞬時不確定度為

Ur=Ur(w)=XrU-XrL。

(14)

顯然，灰自助預(yù)報GBM(1,1)用w=n+1時刻的預(yù)報值描述n時刻的瞬時不確定度。定義 (14) 式的不確定度為時刻n的函數(shù)，即動態(tài)不確定度，其不同于經(jīng)典統(tǒng)計方法的靜態(tài)不確定度，而是隨著時刻n變化而變化。

在動態(tài)測量過程中，假設(shè)每個時間單元中，測量數(shù)據(jù)的總數(shù)為nr=N。如果有hr個數(shù)據(jù)落在估計區(qū)間[XrL,XrU]之外，則定義參數(shù)PrB為

(15)

式中：PrB是給定置信水平P下，對第r個時間單元估計結(jié)果的可靠度，用于描述GBM(1,1)預(yù)報的可信度。一般PrB≠P。根據(jù)PrB的定義知，PrB越大越好，最好是PrB≥P。

由[XrL,XrU]和PrB的定義可知，在w時刻，P越大，則Ur越大。若P=100%，則Ur可以取最大值。必須注意的是，Ur越大，[XrL,XrU]越偏離真值，進(jìn)而估計結(jié)果越失真，因此，第r個時間單元的平均不確定度Urmean需滿足以下條件[11]

(16)

Urmean|m,B,P→min。

(17)

Urmean實際是一個統(tǒng)計量，是動態(tài)測量過程中變量不確定度的均值，可以作為動態(tài)測量過程中隨機(jī)變量波動狀態(tài)的評價指標(biāo)。第r個時間單元表示隨機(jī)變量時間歷程的第r個階段。

從 (16)～(17) 式可以看出，最合適的評估結(jié)果是在PrB=100%的條件下，Urmean得到最小值。因此，在工程實踐中，應(yīng)結(jié)合具體的研究對象，選擇滿足條件的參數(shù)m，B和P。

Urmean對滾動軸承摩擦力矩變異的動態(tài)預(yù)報有重要作用。Urmean越小，摩擦力矩波動越小，軸承性能變異越輕微；否則，摩擦力矩波動越大，軸承性能變異越嚴(yán)重。據(jù)此可以實時預(yù)報軸承性能隨時間的惡化程度，為及時采取相應(yīng)措施、消除安全隱患提供決策建議。

2 試驗分析

試驗軸承型號為608，摩擦力矩測量儀型號為M9908A。摩擦力矩主軸的轉(zhuǎn)速為5 r/min，軸承外圈承受的軸向載荷為15 N。取0.005 86 s為一個采樣間隔單位，獲得了2 048個摩擦力矩數(shù)據(jù)。取前2 000個數(shù)據(jù)作為研究對象，將其分為4組，即得到R=4個時間單元，每個時間單元N=500個數(shù)據(jù)。摩擦力矩的時間序列如圖1所示。由圖可知，在給定的載荷下，摩擦力矩具有變量不確定性的特點，應(yīng)對不同時間單元摩擦力矩動態(tài)不確定度進(jìn)行預(yù)報，從而實現(xiàn)滾動軸承摩擦力矩不確定性的連續(xù)評估，保證系統(tǒng)安全平穩(wěn)運行。

圖1 摩擦力矩的時間序列X1-X4

考慮到模型中m，B和P對Urmean的影響，以第1個時間單元的數(shù)據(jù)X1為例進(jìn)行說明。

對乏信息預(yù)報而言，B太小，雖然預(yù)報速度快且占用內(nèi)存少，但自助再抽樣不充分，計算結(jié)果的波動性增大，直接影響預(yù)報效果；B太大，降低了預(yù)報速度且占用內(nèi)存過多，不利于在線預(yù)報和控制。另外，B過大并不會提高預(yù)報的準(zhǔn)確性，其有一個極限值。通過大量計算可知，當(dāng)B=500時，計算結(jié)果不再波動且計算時間最短。因此，試驗中令B=500，改變m和P，基于灰預(yù)測模型GM(1,1)對時間序列X1進(jìn)行自助預(yù)報，結(jié)果見表1。

表1 m，B和P的選取對Urmean的影響

由表可知，在P1B=100%，且m=6，P=99%時，U1mean取最小值0.022 54。因此，對試驗中其他3個時間單元的試驗數(shù)據(jù)均采取m=6，P=99%，B=500進(jìn)行分析。

分別計算4個時間單元的平均不確定度，結(jié)果如圖2所示。對比圖1和圖2可以看出，對于第r個時間單元，數(shù)據(jù)波動越劇烈，Urmean的值就越大，摩擦力矩變異越顯著。在前2個時間單元中，數(shù)據(jù)波動很小，Urmean(r=1, 2)在0.02～0.025之間，接近0，摩擦力矩變異最小；而在后2個時間單元中存在局部脈沖，第3個時間單元的脈沖波動最劇烈，U3mean最大，摩擦力矩變異最顯著；第4個時間單元的脈沖強(qiáng)度弱于第3個時間單元，U4mean介于U2mean與U3mean之間。顯然，第r個時間單元的平均不確定度Urmean較好地實現(xiàn)了摩擦力矩不確定性的實時預(yù)報，揭示了摩擦力矩在各階段時間歷程的變異程度。

圖2 4個時間序列的平均不確定度

每個時間單元的振動加速度時間序列在各個時刻的估計區(qū)間如圖3所示。由此獲得的時間序列的預(yù)報可靠度如圖4所示。

圖3 時間序列的估計區(qū)間

圖4 時間序列的預(yù)報可靠度

由圖4可知，預(yù)報可靠度的最大值為100%，最小值為99.79%，表明圖3中軸承各個時刻的摩擦力矩數(shù)據(jù)至少有99.79%包絡(luò)在估計區(qū)間內(nèi)，證明所建立的預(yù)報模型能夠有效地預(yù)測軸承摩擦力矩的變化范圍，為及時控制失效提供了有利條件。

3 結(jié)束語

基于GBM(1,1)的滾動軸承摩擦力矩不確定性的動態(tài)預(yù)報方法，解決了概率分布和變化趨勢均未知的條件下摩擦力矩變異的實時預(yù)測問題。該方法可實時追蹤各種隨機(jī)變量的波動和趨勢項的變化，實現(xiàn)了乏信息條件下動態(tài)不確定度的有效預(yù)報，補充了現(xiàn)有不確定性評估理論。