軸承性能時間序列的概率信息求取及退化分析*

常 振，魏劍波，平曉明，曹茂來，王二化

(1.杭州軸承試驗研究中心有限公司，浙江 杭州 310022；2.機械工業軸承產品質量檢測中心(杭州)， 浙江 杭州 310022；3.常州信息技術學院 智能裝備學院，江蘇 常州 213164)

0 引 言

軸承性能時間序列是將某種性能屬性的數值按照時間測量的先后順序排列所形成的數列，根據時間序列所反映出來的發展過程、方向和趨勢，進行類推或延伸，借以預測或分析下一段時間或若干年后可達到的動態運轉水平，進而快速有效地進行產品性能故障診斷[1-3]。

時間序列的信息挖掘一直受到學術界與工程界的高度重視，特別是在航空航天、金融、經濟、天文、地質等應用十分廣泛[4-7]。而早期時間序列分析法的模型幾乎都是線性的，隨產品質量指標的要求不斷提高，研究人員發現非線性模型才能合理解釋工程實際問題[8，9]。軸承服役過程中伴隨有振動、溫度、摩擦力矩等不同屬性的非線性性能時間序列，因此，可從該類時間序列中提取系統參與演化的所有變量的大量信息，從而用于軸承某些性能信號概率信息的求取與退化特征分析。

軸承性能時間序列概率密度函數的求取是數據處理與信號分析的根基，可據此獲取某一性能信號的真值估計與區間波動。

本文基于軸承的非線性性能時間序列，首先通過自助法建立時間序列仿真數據，運用最大熵原理對大量的仿真數據進行概率密度求取，分析其區間波動情況并驗證模型的準確性；然后基于模糊集合理論[10，11]進行軸承性能時間序列的退化特征分析；兩套模型兼并使用綜合討論軸承性能信號的演變規律并識別其性能退化歷程。

1 軸承性能時間序列的概率信息求取

1.1 性能時間序列分組處理

設軸承某一性能時間序列向量表示為：

X=(x(1),x(2),…,x(n),…x(N))

(1)

式中：x(n)—軸承性能時間序列中的第n個數據，n=1,2,…,N；N—樣本個數。

取前n個數據序列作為訓練組，記為XA；后N-n個數據序列作為驗證組，記為XB，即該軸承的性能時間序列分組為訓練組XA和和驗證組XB，表示為：

XA=(x(1),x(2),…,x(n))

(2)

XB=(x(n+1),x(n+2),…,x(N))

(3)

1.2 訓練組XA灰自助再抽樣

為滿足灰預報模型GM(1,1)關于x(n)≥0的苛刻要求，在式(2)中，若有x(n)<0，則人為選取一個常數c，使得x(n)+c≥0即可。

實際分析時，X表示為：

XA=(x(n)+c;t=1,2,…,N)

(4)

從XA中取與時刻n緊鄰的前m個時刻的數據(包括時刻n的數據)，構成時刻n的動態分析子向量，即：

Xm=[xm(u)];u=n-m+1,n-m+2,…,n;n≥m

(5)

運用自助法，在時刻n，從Xm中等概率可放回地隨機抽取一個數，共抽取m次，可得到一個自助樣本Y1，它有m個數據。

按此方法重復執行B次，得到B個樣本，可表示為：

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

(6)

式中：Yb—第b個自助樣本；B—總的自助再抽樣次數，也是自助樣本的個數。

且有：

Yb=[yb(u)]

(7)

其中：u=n-m+1,n-m+2,…,n；b=1,2,…,B。

根據灰預報模型GM(1,1)[12]，可得到ω=n+1時刻，組成樣本容量為B的大數據模擬信號XGBootstrap，可表示為如下向量：

(8)

1.3 基于最大熵原理求解模擬信號XGBootstrap的概率密度函數

(9)

隨后可根據最大熵方法獲得基于模擬信號XGBootstrap的密度函數的最優估計。最大熵方法的核心原理是在所有的可行解中，滿足熵最大的解是最“無偏”的，亦即：

(10)

在具體的求解過程中，引入拉格朗日乘子，可以把概率分布求解問題轉化為對拉格朗日乘子的求解問題[13]。

1.4 積分區間的映射求解概率密度函數

為了使計算過程的收斂，筆者將序列XGBootstrap中的數據按進行遞增排序并分成β組，隨后畫出對應的直方圖。同時，可得到組中值zμ和頻數Γμ，u=2,3,…,β+1。然后將直方圖擴展成β+2組，并令Γ1=Γξ+2。將原始數據區間S映射到區間[-e,e]中，e—自然常數。

令：

(11)

(12)

(13)

b=e-azξ+2

(14)

式中：e—自然常數，e=2.718 282。

因此，所求的最大熵原理構建的軸承性能時間序列概率密度函數變換為：

(15)

式中：j—原點矩階數，常用j=5；c0—首個拉格朗日乘子；ci—第i+1個拉格朗日乘子，i=1,2,…,j。

對于實數α雙側分位數，有概率：

(16)

(17)

式中：XU，XL—軸承性能時間序列估計區間的上界值和下界值；[XL,XU]—估計區間。

1.5 軸承性能時間序列概率信息求取方法的驗證

已知該軸承性能時間序列的前n個數據序列作為訓練組XA，后N-n個數據序列作為驗證組XB。其中利用訓練組XA獲取上述性能時間序列的概率密度函數和對應的區間估計，則用驗證組XB來證實上述方案的可行性：

假設在N-n個驗證組XB的數據中有y個數據在估計區間[XL,XU]之內，則上述模型的可行度或準確度表示為：

(18)

此處所提準確度也是指驗證組XB中的數據落入估計區間的頻率，實現了所提模型的自我驗證；同時也獲取了該軸承性能時間序列的區間上、下限，可以為工程中及時掌控軸承的服役工況提供理論支持。

2 軸承性能時間序列的退化特征分析

2.1 性能時間序列向模糊關系的轉換

軸承性能退化特征分析時，使用的是模糊等價關系轉換。然而，在工程實際中往往得到的是模糊相似關系，必須將這種模糊相似關系轉變為模糊等價關系。傳遞閉包法可以進行這種關系的轉換，以獲得軸承性能時間序列的模糊等價關系，具體求解步驟與方法如下：

將軸承性能時間序列分組處理：平均分為p組，設每組有q個數據，即一小組可組成數據序列：

Zl=(Zl1,Zl2,…,Zlk,…,Zlq)k=1,2,…,q,l=1,2,…,p

(19)

式中：p—樣本數；q—每個樣本的樣本量；Zlk—第l個樣本的第k個數據。

則該軸承性能時間序列構成一個新的集合，即：

(20)

對于任意的模糊關系R，如若存在如下關系：

T(R)=Rh-1=Rh=…h=1,2,3,…

(21)

式中：T(R)—模糊數學中的模糊等價關系。

T(R)可以按照上述式子逐步求出，具體過程如下：

第一步，求出R2=R°R；

第二步，求出R4=R2°R2；

……

第u步，直到求出R2u=Ru為止。

其中：運算符“°”—矩陣的模糊運算M(∧,∨)；“∨”—“或”運算取最大；“∧”—“與”運算取最小；Ru—所求的該軸承性能時間序列模糊等價關系T(R)，也是模糊集合理論的傳遞閉包。

即有：

T(R)=Ru

(22)

有：

(23)

其中：0≤αil≤1，且滿足：

(24)

式中：αil—性能時間序列分組后第i個樣本Zi和第l個樣本Zl兩者之間的模糊等價關系，亦即Zi樣本信息和Zl樣本信息的符合程度。

2.2 軸承性能退化依據

在工學、文學、哲學等各大領域中，0和1都可以表示事物的真假、有無、強弱等兩個極端狀態。然而任何事物都有著或多或少的關聯，很難用0和1的確定性規則判定兩個事物的必然規律。

可以在模糊等價關系αil中設定λ閾值來診斷兩者關系密切的大小或者產生特征退變存在的顯著性：

(1)若αil>λ，則兩段軸承的性能時間序列Zi和Zl在λ閾值水平下彼此之間的相似度很高，關聯程度較大，即表明兩個系統間不存在性能退化；

(2)若αil≤λ，則兩段軸承的性能時間序列Zi和Zl在λ閾值水平下彼此之間的相似度很低，關聯程度較小，即表明兩個系統間變化大、演變強烈，存在有性能退化趨勢。

根據模糊集合理論，λ為研究對象從一個極端屬性過渡到另一極端屬性的邊界，也是判定兩段軸承性能時間序列是否有退變產生的顯著性。在工程實際中，當λ=0.5時兩系統之間的模糊關系達到最大，介于最難分辨的真假、有無、強弱之間；當λ>0.5時兩系統間關系趨于清晰，相似度較高且可判定兩者未發生退變[14，15]。

所以在應用過程中，取λ=0.5來界定軸承性能是否發生特征退化。

2.3 軸承性能退化的模糊表征

定義退化系數集合為：

U=(u1,u2,…,uξ,…,uw-1)

(25)

其中，有：

(26)

式中：uξ—軸承性能系數，即模糊等價關系αil的分段均值；w—樣本含量；ξ—各樣本采樣的時間先后順序，是一個時間變量。

根據退化系數uξ來判斷軸承性能時間序列的演變歷程：(1)如若uξ>λ=0.5，則表明軸承服役期間保持良好的運轉特性，性能未發生而惡性退化；(2)uξ<λ=0.5，則表明其服役期間各段性能數據相似度下降，退變趨勢明顯，性能發生明顯退化甚至軸承已發生損壞。

3 實驗及結果分析

案例一：滾動軸承振動性能時間序列試驗研究。

該試驗在ABLT-1型軸承壽命強化試驗機上進行的，試驗軸承為6 208。轉速為4 000 r/min，純徑向載荷Fr=8 kN，油潤滑。每10 min采集振動均方根值一次，單位ms-2，共采集60個數據，獲得振動性能時間數據序列X振動，如圖1所示。

圖1 軸承振動性能時間數據序列X振動

案例二：滾動軸承溫度性能時間序列試驗研究。

該試驗在ABLT-3型軸承壽命強化試驗機上進行的，試驗軸承為6 203。轉速為7 000 r/min，純徑向載荷Fr=2.5 kN，油潤滑。每10 min采集溫度值一次，單位℃，共采集60個數據，獲得溫度性能時間數據序列X溫度，如圖2所示。

圖2 軸承溫度性能時間數據序列X溫度

案例三：滾動軸承摩擦力矩性能時間序列試驗研究。

研究對象為軸承組件的穩態電流信號(單位：mA)，反作用控制箱輸出指令電壓帶動真空實驗裝置中的滾動軸承轉動，由反饋裝置取樣并轉換后將得到的電流信號反饋給控制箱。反饋得到的穩態電流信號間接得到滾動軸承運轉期間摩擦力矩時間序列。等間隔地采集軸承摩擦力矩60次，獲得摩擦力矩性能時間序列X摩擦力矩，如圖3所示。

圖3 軸承摩擦力矩性能時間數據序列X摩擦力矩

3.1 概率密度求取

案例中性能時間序列X，筆者取前10個數據作為訓練組，記為XA。從訓練組XA序列中等概率可放回地進行抽樣，每次抽樣個數v=6；連續上述步驟重復抽取B=10 000次，進而可組成樣本容量為10 000的大數據模擬信號XGBootstrap。

案例一生成的灰自助數據如圖4所示。

圖4 振動性能信號生成的灰自助數據XGBootstrap

案例二生成的灰自助數據如圖5所示。

圖5 溫度性能信號生成的灰自助數據XGBootstrap

案例三生成的灰自助數據如圖6所示。

圖6 摩擦力矩性能信號生成的灰自助數據XGBootstrap

案例一求得的概率密度函數如圖7所示。

圖7 振動性能信號概率密度函數

案例二求得的概率密度函數如圖8所示。

圖8 溫度性能信號概率密度函數

案例三求得的概率密度函數如圖9所示。

圖9 摩擦力矩性能信號概率密度函數

3.2 區間預報

根據圖(7～9)概率密度函數信息(由前n=10個數據XA獲得)，并給定α=0.002 7顯著性水平下可獲得軸承性能信號的區間估計，相應置信水平為99.73%。用剩余的60-n=50個數據來檢驗該預報的可靠性，并掌控后續信號的退變信息。

其結果如表1所示。

表1 時間序列X振動、X溫度、X摩擦力矩的區間預報分析

據上述滾動軸承性能時間序列概率密度信息的求取過程可知，數據處理過程中僅用了訓練組的10個數據，且對該組數據無任何先驗信息也無任何假設約束條件，表明該模型具有優越的客觀主動性；

隨后以驗證組來證明其區間預報的準確性，在99.73%的置信水平下，使用灰自助最大熵法區間估計的上、下邊界差值較小，其中，性能時間序列X振動的區間差值為0.42；性能時間序列X溫度的區間差值為2.34；性能時間序列X摩擦力矩的區間差值為23.42，同時表明了所提方法具有較高的預報精度。更何況滾動軸承振動、溫度、摩擦力矩性能時間序列的隨機性很強，且具有明顯的不確定度，其區間預報十分復雜，所以該方法的區間預測已相當可觀。

此外，在驗證組的50個數據進行方法驗證過程中，試驗數據X振動、X溫度中驗證組落在預報區間之外的個數極少，分別為2個和1個，表明預報結果的準確可行性，可用于工程實踐的在線監測與診斷。然而，試驗數據X摩擦力矩中驗證組落在預報區間之外的個數很多，為33個，并非表明預報結果不準確可行，而是表明驗證組的后50個數據相對于訓練組的前10個數據有著明顯的差異性，亦即軸承服役過程中隨時間的演變，其摩擦力矩性能信號產生了嚴重的退化。由圖3的原始數據X摩擦力矩可知，摩擦力矩信號呈現上升趨勢，即表明軸承運行狀況并非穩定，內部摩擦逐漸增大，其性能信號已發生變異或退化。

基于灰自助最大熵法X振動、X溫度的區間誤報率分別為4%和2%，即表明此次性能時間序列預報的可靠度在96%～98%，所以該方法進行時間序列區間的預報精度高且可靠性強；同時也說明后50個驗證組組數據與前10訓練組數據保持著良好的一致性，退化跡象十分不明顯，即軸承運轉期間保持較高的穩定性和平穩性。

基于灰自助最大熵法X摩擦力矩的區間誤報率66%，即表明此次性能時間序列預報的可靠度在44%，結果表明后50個驗證組組數據與前10訓練組數據一致性不強，退化跡象較為明顯，即軸承運轉期間已產生明顯退化和變異。

3.3 軸承退化性能分析

案例一。

數據處理后得到的模糊等價關系矩陣T(R)為：

(27)

穩定系數集合U=[0.71 0.72 0.70 0.72 0.70 0.71 0.68 0.66 0.60]，結果如圖10所示。

案例二。

數據處理后得到的模糊等價關系矩陣T(R)為：

圖10 軸承振動性能時間序列X振動退化系數

(28)

穩定系數集合U=[0.64 0.61 0.61 0.61 0.63 0.61 0.61 0.61 0.61]，結果如圖11所示。

案例三。

數據處理后得到的模糊等價關系矩陣T(R)為：

圖11 軸承溫度性能時間序列X溫度退化系數

(29)

穩定系數集合U=[0.70 0.67 0.64 0.62 0.59 0.56 0.51 0.50 0.48]，結果如圖12所示。

圖12 軸承摩擦力矩性能時間序列X摩擦力矩退化系數

由上述軸承性能時間序列X振動、X溫度、X摩擦力矩的3個案例的退化信息處理結果可知：

案例一、二的振動性能和溫度性能表現穩定，最小退化系數分別為uξ=0.600和uξ=0.609，兩者均大于λ=0.5閾值；則表明案例一的軸承振動性能和案例二的溫度性能未發生明顯異常，前后時間序列未發生明顯變異，軸承退化特征不明顯，軸承運轉狀況穩定可靠，健康狀況較為良好；

案例三的摩擦力矩性能表現不良，最小退化系數為uξ=0.477<0.5閾值，則表明案例三軸承的摩擦力矩信號已發生明顯異常，前后時間序列有明顯變異，軸承性能信號已發生明顯退化跡象，圖12中退化系數呈現出明顯的減小趨勢，該結果與其原始數據系列逐漸增大有關，且與3.2部分區間預報分析結果保持良好的一致性，均說明了該軸承摩擦力矩性能信號退化明顯，運轉狀況不再安全可靠；為避免惡性事故的發生，應及時維護或更換軸承，并對其進一步健康監控。

綜上可知，軸承性能時間序列的概率信息的求取為更深層次挖掘其隱含信息奠定了基礎，據此進行的區間預測實現了模型準確性的自我驗證，并可以較好地反映其工作期間性能信號的波動情況；另外，軸承性能時間序列的退化系數，可以很好地識別出其特定性能的演變規律和退化趨勢。

因此，筆者所提出的自助最大熵結合法可以有效地擬合軸承性能時間序列的概率密度函數，并可預測其性能區間波動狀態；基于模糊理論的退化系數可準確地監控其性能退化動態信息。

4 結束語

針對軸承性能時間序列概率信息求取及退化分析問題，筆者對軸承振動、溫度、摩擦力矩3種性能時間序列進行了研究，結論如下：

(1)在概率信息求取過程中，X振動、X溫度區間誤報率低，概率信息求取較為準確可靠；X摩擦力矩的誤報率為66%，這主要是因為該性能信號具有明顯的退化信息，而非所提模型不可靠；

(2)在退化性能分析過程中，X振動、X溫度表現出良好的服役狀況；X摩擦力矩的退化系數uξ=0.477<0.5，具有明顯退化跡象；與概率信息求取的結果保持較好的一致性；

(3)結合軸承性能時間序列的概率信息和退化系數，可有效評估軸承性能時間序列的變化趨勢本質特征和退化演變跡象，實時掌控軸承的健康狀況，并及時發現其失效隱患。