一種用于軸承故障診斷的特征選取方法
【技術領域】
[0001] 本發(fā)明涉及一種機械故障診斷方法，尤其是一種用于軸承故障診斷的特征選取方 法。
【背景技術】
[0002] 軸承是機械系統的重要部件之一，軸承故障診斷越來越受到重視，已成為保證生 產安全穩(wěn)定運行和提高產品質量的重要手段和關鍵技術。當前軸承故障特征提取方法大多 依賴于某種（類）故障的先驗知識，專門針對一種或幾種故障，而由于故障類型、工況環(huán)境 以及故障成因的多樣性，其故障特征也各有差異，故面對復雜多變的故障環(huán)境，需要多種特 征提取方法配合才有所成效，造成計算量與資源浪費，其根本原因是缺乏通用的，不隨故障 轉移的特征提取方法。
[0003] 同時軸承故障特征提取方法的偶然性導致了特征選擇的不確定性，甚至無法實施 特征選擇過程，目前特征選擇方法大多基于后續(xù)診斷的正確率效果反饋，以剔除弱故障區(qū) 分能力的特征（或選擇強故障區(qū)分能力的特征），并重新建模。故特征對診斷效果先期估計 與評判能力不足，且反復建模使得診斷時間代價較大，模型難以適應動態(tài)變化的參數和環(huán) 境變化，缺乏診斷實時性。

【發(fā)明內容】

[0004] 為了克服現有故障診斷特征提取和選擇方法的不足，本發(fā)明提供了一種用于軸承 故障診斷的特征選取方法。所述方法利用奇異值分解（SVD)進行故障信號特征提取，并定 義特征區(qū)分度公式進行特征選擇，具體采用以下技術方案：
[0005] -種用于軸承故障診斷的特征選取方法，包括如下步驟：
[0006] 步驟（1):采集軸承不同故障類型的多個振動信號；所述不同故障類型包括不同 故障直徑和不同故障部位；根據所采集的振動信號構建各個振動信號自相關函數的矩陣， 并對所構建的矩陣進行奇異值分解，得到最大奇異值元素；
[0007] 步驟（2):逐步增加矩陣維數，直至矩陣的最大奇異值元素保持穩(wěn)定；得到此時各 個振動信號的奇異值向量和矩陣維數值；
[0008] 步驟（3):根據步驟（2)得到的各個振動信號的奇異值向量和矩陣維數值得到不 同維數下各個振動信號的奇異值向量的平均曲率；根據各個振動信號的奇異值向量的平均 曲率定義反映不同類型的軸承振動信號差異的區(qū)分度公式：
[0010] 其中，dab為第a類故障類型和第b類故障類型的特征區(qū)分度（a e {1，2，···，γ}; b e {1，2,…，r})，r為故障類型總數；kJPkb為這兩類故障特征向量曲線的平均曲率；q為 分類性能指標參數統計值；
[0011] 步驟（4):逐步減小矩陣維數，計算矩陣奇異值向量，并根據步驟（3)定義的區(qū)分 度公式計算不同類型的軸承振動信號之間的區(qū)分度；當所有軸承振動信號之間的區(qū)分度最 小值為1時，此時的矩陣維數即為特征選取維數，同時得到的該維數下的奇異值向量即為 選取的特征向量。
[0012] 所述步驟（1)采集的振動信號：
[0013] yid= {x 丄（i = 1，2,…，k)
[0014] 其中，yld為第id (id = 1,2, 3···，r)類故障的軸承振動信號序列；r為總故障類別 數；Xi= x(t J為h時刻振動加速度數據；t ;= t o+Q-l) Δ t(i = 1，2,…，k) ;t。為初始采 樣時刻，△ t = 1/f為相鄰采樣時間間隔；f為采樣頻率；k為總采樣數；i為采樣計數。
[0015] 所述步驟（2)中，逐步增加矩陣維數，直至矩陣的最大奇異值元素保持穩(wěn)定過程 如下；
[0016] 步驟（21):周期擴展振動信號yld，并計算擴展信號的自相關函數
η為擴展后的振動信號第η個序列點，m為 計算自相關函數時的平移點數；選取較小的維數構造初始自相關函數Hankel矩陣C。：
[0018] 其中，C。為初始自相關函數Hankel矩陣，M?！?。，％>0為二維矩陣的初始選取維數；
[0019] 步驟（22)，對C。進行奇異值分解：C。= U1AV1/;其中，矩陣 ^jQ - I U^,U2,......G 。為 Μ。XM。維的正父矩陣，w」，，.:，為矩陣兀素向 量；1/?二(??e兵她_為NqXNq維的正交矩陣，1??:·、, Vwq為矩陣元素向 量；D。為對角陣：
[0021] 其中，H " "，％為自相關矩陣C。的奇異值，且4 > A2 > ··· >心〇，故
[0022] 步驟（23)，令Mt= M〇+t，Ns= N〇+s，（t = 1，2, 3..·，s = 1，2, 3..·)，Qp= min{Mt，Ns}， p = min{t, s};其中，t與s為矩陣橫向和矩陣縱向的擴展維數，每次交替增加 I進行擴 展；重復步驟（21)和步驟（22)，直到得到最大奇異值元素 A1達到穩(wěn)定；此時奇異值向量 /二（H.……，A(?t)為提取的特征向量。
[0023] 所述步驟（23)中得到最大奇異值元素 λ ^勺條件為若維數Q p時最大奇異值元素 等于維數Qp+1時的最大奇異值元素，或兩者之差小于〇. 01 λ ^直，此時Q p的最大奇異值元素 為穩(wěn)定的最大奇異值元素；其中Qp是指擴展后的維數；Q P+1是指Q p再次擴展后的維數。
[0024] 本發(fā)明的有益效果是：本發(fā)明提出了一種用于軸承故障診斷的特征提取和選擇方 法，所述方法利用奇異值分解（SVD)進行故障信號特征提取，能夠有效克服依賴故障特征 先驗知識的局限，并提升軸承故障診斷通用性；所述方法定義特征區(qū)分度公式，能夠有效評 判方法的故障診斷能力；特征選取能夠有效提升后續(xù)診斷正確率，平衡診斷環(huán)境適應性與 實時性要求。本發(fā)明為軸承故障診斷提供了低維的，有明確區(qū)分的特征，在變工況負載等軸 承故障診斷應用領域具有潛在經濟和社會價值。
【附圖說明】
[0025] 圖1為本發(fā)明的用于軸承故障診斷的特征提取和選擇算法流程圖；
[0026] 圖2為本發(fā)明的目標振動數據與輔助振動數采集示意圖（a :故障部位；b :故障直 徑）；
[0027] 圖3為本發(fā)明的軸承振動信號示例圖（健康軸承，12kHz，驅動端對象）；
[0028] 圖4為本發(fā)明的最大特征值與自相關矩陣維數關系曲線（故障類型：輕微，嚴 重）；
[0029] 圖5為本發(fā)明的SVD提取奇異值向量的對數曲線圖（a :故障直徑；b :故障部位）；
[0030] 圖6為本發(fā)明的奇異值特征向量的平均曲率曲線（a :電機轉速；b :故障直徑；c : 故障部位；d :試驗對象）；
[0031] 圖7為本發(fā)明的區(qū)分度與分類性能指標參數比值曲線（故障直徑）；
[0032] 圖8為小波變換特征提取與選擇方法的特征頻率的信號能量曲線（a :145. 83Hz ; b :157. 94Hz ;c :315. 89Hz ;d :473. 83Hz)。
【具體實施方式】
[0033] 下面結合附圖及具體實施例對本發(fā)明進行詳細說明。
[0034] 參見圖1，一種用于軸承故障診斷的特征提取和選擇方法，包括如下步驟：
[0035] 步驟（1):對含不同故障類型的軸承進行振動信號（加速度）采集，用于特征提 取、選擇以及故障區(qū)分性能測試；
[0036] 步驟（1)中，不同故障類型，包括不同故障直徑與不同故障部位的信號采集，其中 軸承振動信號采集示意圖如圖2所示：
[0037] a)例如，故障直徑包括：0·000_、0· 178mm、0. 356mm、0. 533mm直徑，分別定義為健 康、輕微故障、中度故障、嚴重故障；
[0038] b)例如，故障部位包括：內圈、滾珠、外圈6:00、外圈3:00、外圈12:00(負載6:00 方向）；
[0039] c)實測中，采集頻率f = 12kHz或48kHz條件的振動信號，在本發(fā)明中，采集頻率 為軸承轉速的200倍以上。
[0040] Yid= {x il, (i = 1, 2, ···, k)
[0041] 其中，yld為第id(id= 1,2,3···,γ)類故障的含觀測噪音軸承振動數據序列；r為 總故障類別數；Xi= x(t J為1^時刻振動加速度數據；t ;= t o+Q-l) Δ t (i = 1, 2,…，k); t。為初始采樣時刻，△ t = Ι/f為相鄰采樣時間間隔；k為總采樣數；i為采樣計數。
[0042] 軸承振動信號示例圖如圖3所示，可以看出，在變化工作環(huán)境下，振動信號呈現周 期非平穩(wěn)特性，并且混有高斯噪聲，使得軸承的故障特征易被正常振動引起干擾信號以及 白噪聲掩蓋，增大提取軸承故障特征提取的復雜性和困難度。
[0043] 步驟（2):利用Hankel矩陣構建振動信號自相關函數的自相關矩陣C，增大矩陣維 數，觀察其奇異值分解（SVD)得到的最大特征值元素變化，直至保持恒定，并得到穩(wěn)定時的 振動信號特征向量（奇異值向量）與矩陣維數值；
[0044] 步驟（2)中，由于軸承故障發(fā)生時，損傷點與其它表面元件碰撞產生信號沖擊，可 能覆蓋軸承系統和傳感器的固有頻率，故自相關矩陣奇異值分解的實施目的是，克服大多 軸承故障診斷方法依賴故障特征先驗知識（故障頻率、噪音干擾等）的局限，提供不依賴先 驗知識，并具備物理意義的清晰特征，提升軸承故障診斷通用性。
[0045] 步驟⑵中，自相關矩陣C的奇異值分解步驟如下：
[0046] 步驟（21)，周期擴展振動信號yld，并計算擴展信號的自相關函數
η為擴展后的振動信號第η個序列點，m為 計算自相關函數時的平移點數；選取較小的維數構造初始自相關函數Hankel矩陣C。：
[0048] 其中，C。為初始自相關函數Hankel矩陣，M?！?。，％>0為二維矩陣的初始選取維數；
[0049] 步驟（22)，對C。進行奇異值分解：C。= UQDQVQT;其中，矩陣
為Μ。X M。維的正父矩陣，?%; ，.4?..為矩陣兀素向 量；
為Ν〇ΧΝ。維的正交矩陣，…，％。為矩陣元素向 量；D。為對角陣：
[0051 ] 其中，......，心。為自相關矩陣C。的奇異值，且4 > I2 > .1(?、，故
[0052] 步驟（23)，令Mt= M〇+t，Ns= N〇+s，（t = 1，2, 3...，s = 1，2, 3...)，Qp= min{Mt，Ns}， p = min {t, s};利用MjP N s維數重復步驟B-I和B-2,得到最大奇異值元素 λ i，若維數Qp時最大奇異值元素等于（