基于灰自助GBM(1,1)的滾動軸承摩擦力矩不確定性的動態預報方法

時保吉，栗永非

(新鄉職業技術學院，河南 新鄉 453006)

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

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

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

1 數學模型

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

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

(1)

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

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

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

Xrm={xrm(u)+qr};

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

(2)

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

在時刻n，從Xrm中等概率可放回地隨機抽取一個數據，抽取m次，得到一個大小為m的自助樣本[9]。連續重復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個自助再抽樣數據。

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

(4)

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

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

灰生成模型為

(5)

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

用增量代替微分，即

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

(6)

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

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的預測值為

(9)

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

(10)

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

frw=frw(xrm) ，

(11)

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

假設顯著性水平為α，則置信水平為

P=(1-α)×100%，

(12)

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

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

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

(13)

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

定義n時刻的擴展不確定度，即瞬時不確定度為

Ur=Ur(w)=XrU-XrL。

(14)

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

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

(15)

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

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

(16)

Urmean|m,B,P→min。

(17)

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

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

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

2 試驗分析

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

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

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

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

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

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

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

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

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

圖3 時間序列的估計區間

圖4 時間序列的預報可靠度

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

3 結束語

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