滾動軸承故障診斷的最優加權路圖GFT沖擊提取方法

王好將 于德介 高藝源

摘要：滾動軸承發生局部故障時將產生由共振頻率調制的周期瞬態沖擊，有效提取沖擊特征是診斷滾動軸承故障的關鍵。圖信號處理方法（Graph Signal Processing，GSP）是基于圖譜理論發展起來的新研究領域，將振動信號轉換為圖信號進行分析，能有效揭示振動信號特性。對高斯函數加權下的路圖拉普拉斯矩陣進行特征分解，發現代數連通度（Algebraic Connectivity）以內的特征向量存在明顯的沖擊，因此提出利用代數連通度以內的特征向量結合逆圖傅里葉變換（Graph Fourier Transform，GFT）重構故障信號中沖擊分量的方法。高斯加權函數中的熱核寬度決定沖擊特征向量的分布，直接影響重構結果，為解決熱核寬度的選擇問題，提出結合粒子群算法（ParticleSwarm Optimization，PSO）確定最優熱核寬度;然后利用最優特征向量組重構沖擊信號，并進行包絡解調;最后實現滾動軸承故障的有效診斷。算法仿真和應用實例表明，基于最優加權的路圖GFT方法能有效地重構滾動軸承故障沖擊特征，診斷故障類型。

關鍵詞：故障診斷;滾動軸承;沖擊提取;路圖;高斯加權函數

中圖分類號：TH165+。3;TH133.33;7N911.7文獻標志碼：A 文章編號：1004-4523（2020）03-0604-10

DOI：10.16385/j.cnki.issn.1004-4523.2020.03.020

引言

滾動軸承作為旋轉機械重要的承載零部件，其發生故障時直接影響整臺機械的性能，為避免重大故障的發生，需要監測其運行狀態。軸承內圈、外圈和滾動體存在局部故障時，伴隨著軸承的周期性旋轉，會以一定故障通過頻率產生由共振頻率調制的周期瞬態沖擊。因此，滾動軸承的故障信息存在于周期性沖擊成分中，從軸承故障振動信號中提取周期性沖擊特征是有效診斷軸承故障的關鍵。但在實際工程中，受到工作環境、振動傳輸路徑和多振動源相互耦合等因素影響，導致軸承振動信號中往往存在嚴重的背景噪聲干擾及工頻振動、諧波振動等強故障特征信號，此時直接對故障信號進行解調分析無法發現弱故障特征，容易造成故障的漏診。因此，尋找行之有效的沖擊特征提取方法一直是滾動軸承故障診斷研究的重點。

共振解調方法是提取滾動軸承故障沖擊特征的有效方法。陳祥龍等提出了基于平方包絡譜相關峭度指標的最優共振解調方法，穩健地確定最優頻帶中心頻率和帶寬，準確提取循環瞬態沖擊信號。Antoni提出了基于FIR帶通濾波器的快速譜峭度算法（SK），通過設計二叉樹結構的帶通濾波器對信號進行濾波，對各層濾波結果計算峭度值得到快速譜峭度圖，根據峭度最大原則選取最優解調分析頻帶，提取沖擊特征信號。目前譜峭度方法已廣泛應用于滾動軸承故障診斷。

近年來，以圖信號為研究對象的圖信號處理方法開拓了信號處理方法新領域，逐漸引起國內外學者的關注。圖譜理論（Graph Spectrum Theo-ry）是圖信號處理方法的基礎，其主要內容是結合代數方法分析圖信號。目前，圖譜理論已經廣泛應用于機器視覺、網絡學習、多傳感器融合等領域，不但在故障診斷領域中的應用研究則相對較少。圖傅里葉變換方法將圖信號正交投影于拉普拉斯矩陣的特征向量上，為圖結構數據的分析提供了類似“頻率”的概念，從而將圖信號研究從頂點域拓展到圖譜域，并為圖信號頂點域中卷積、平移及調制運算的實現提供了可能。為了將圖信號處理方法引入到振動信號分析中，需要將振動信號時間序列轉換為合適的圖信號。路圖僅考慮相鄰時間點間存在邊連接關系，與時間序列存在一一對應關系。研究發現，無權路圖的拉普拉斯矩陣特征向量與離散余弦變換的基函數相同，路圖傅里葉變換也與經典FT結果有對應關系。

Ou等研究發現，在歐氏距離加權下，路圖傅里葉變換將沖擊成分投影在高階次區域，并利用這一特點實現了沖擊特征提取。然而，歐氏距離加權的形式與一般加權的定義相違背，即距離越遠權值越小，另外沖擊成分出現在高階次的現象也未從理論上解釋，限制了該方法的進一步應用。本文研究發現，在高斯函數加權下，與路圖拉普拉斯矩陣零特征值對應的特征向量呈現沖擊現象，這一特點與圖上連通分量相關聯，對應的路圖傅里葉變換將沖擊成分投影在低階次區域，尤其是代數連通度（Alge-braic Connectivity）以內的低階次區域。針對上述特點，可利用代數連通度以內的拉普拉斯矩陣特征向量通過路圖傅里葉逆變換重構沖擊信號。此外，熱核寬度是高斯函數中的重要參數，影響零特征值數目及沖擊特征在階次域的能量分布，不合理的熱核寬度將導致無法重構沖擊特征。因此，選擇合適的熱核寬度值是有效重構沖擊信號分量的關鍵。

基于高斯函數加權下路圖傅里葉變換的特點，本文提出一種結合粒子群優化算法的最優加權路圖傅里葉變換沖擊提取方法。該方法先將振動信號轉換為路圖;再以重構沖擊信號包絡譜故障頻率處的幅值為適應度函數，通過粒子群算法尋找最優熱核寬度值，在熱核寬度最優時，利用代數連通度以內的拉普拉斯矩陣特征向量重構沖擊信號分量;最后對重構沖擊信號進行包絡解調分析，診斷滾動軸承故障類型。滾動軸承數值仿真故障信號和應用實例表明，與文獻方法及譜峭度方法比較，該方法能更有效地提取軸承故障信號中的沖擊特征，實現背景噪聲干擾情況下的軸承故障的有效診斷。

1 GFT與高斯函數加權路圖

路圖僅考慮時間序列中相鄰時間點的連接關系，與時間序列的連續發展特點一致。因此，本文將振動信號轉換為路圖。不同加權方式下，圖拉普拉斯矩陣的特征向量代替了傳統三角函數基函數，為基于圖傅里葉變換的特殊信號分量提取提供了可能。本節首先介紹圖信號處理的基礎概念，基于圖給出GFT的定義;再結合典型路圖拉普拉斯矩陣特征分解的特點，解釋可以用于沖擊信號提取的原因;最后對高斯加權熱核寬度參數對沖擊信號重構的影響作了分析。

1.1 圖與圖傅里葉變換

圖由頂點通過頂點間連邊構成，一般的，無向、加權圖可表示為G=（V，E，A），其中V表示圖頂點的有限集合（頂點數量|V|=N），E為連邊的集合，A表示鄰接矩陣，其中元素a_ij為頂點v_i與v_j間連邊e_ij的加權值。特別地，對于無向圖，A為對稱矩陣。假定x ∈ R^N為定義在圖G上的信號，X_n為圖G中第，n個頂點的信號值。則頂點間權值通過頂點間距離d_ij=x_i-x_j進行構造，本文采用高斯函數加權

圖傅里葉變換（GFT）定義為圖信號關于拉普拉斯矩陣特征向量函數的展開，其主要思想是將圖信號正交投影于拉普拉斯矩陣的特征向量上，從而實現圖信號的“頻率”分析，將圖信號研究從頂點域拓展到圖譜域。

按式（2）構建拉普拉斯矩陣并進行正交分解

Lu_l=λ_lu_l，l=0，1，2…，N-1（3）

式中 l為階次;λ_l為第l階特征值;u_l為第l階特征向量。記特征向量矩陣U=[u₀u₁…u_N-1;特征值對角矩陣A=diag Lλ₀λ₁…λ_N-I]。

結合拉普拉斯矩陣的結構可知其特征分解有以下特點：

（1）L為對稱半正定矩陣，因此，特征向量矩陣為單位正交矩陣，U^TU=i.

（2）L的最小特征值為0，其數量為k（k≥1），將其按升序排列，0=λ₀…=λ_k-1<λ_k≤…≤λ_n，其中λ_k為第一個大于零的特征值，稱其為代數連通度（Algebraic Connectivity）。

（3）零特征值的數量k等于圖中連通分量的個數。連通分量定義為圖的一個子圖，且子圖中任意兩個頂點連通。

基于拉普拉斯矩陣及其特征向量可定義GFT，與傳統傅里葉變換定義類似，將圖拉普拉斯矩陣特

1.2 路圖拉普拉斯矩陣特征分解與沖擊特征

路圖是時間序列最為吻合的構圖對象，僅相鄰時間點間存在邊連接。考慮如圖1所示包含9個頂點的簡單路圖，對應時間序列為

x={8，1，1，8，1，1，8，1，1} （7）

該時間序列可看作沖擊信號。理想情況下，設定距離較遠的頂點間不連接且權值為0，否則為1.則僅頂點2與3，5與6，8與9存在邊連接，且權值a_ij=1.顯然，圖1中存在6個連通分量，分別為{{v₁}，{v₂，v₃}，{v₄}，{v₅，v₆}，{v₇}，{v₈，v₉}}

。

按式（2）計算拉普拉斯矩陣并進行標準正交分解。特征值對角矩陣如下式所示

A=diag[000000222]（8）

可見零特征值的數量與圖1中連通分量數目相同，代數連通度λ_k=2此外，圖1中頂點v₁，v₄，v₇，呈現的沖擊特征，分別對應第1，3，5個孤立連通成分，在特征向量矩陣中體現為特征向量u₀，u₂，u₄為單位脈沖向量，且沖擊發生時刻與頂點位置一致，特征向量的沖擊特性如圖2所示。因此，原始時間序列中的沖擊特征可以由特征向量u₀，u₂，u₄進行路圖傅里葉逆變換重構。對于實際信號，由于存在噪聲與其他分量信號干擾，考慮利用零特征值對應的所有特征向量重構信號中的沖擊分量，如下式所示

1.3 熱核寬度參數對沖擊特征重構結果的影響

對實際信號，由于無法直接判定頂點間的連接關系，故采用式（1）所示的加權方式衡量頂點間連接關系與權重，熱核寬度值δ將影響沖擊特征的重構結果。因此，需要分析熱核寬度參數對沖擊特征重構結果的影響。以圖3所示仿真沖擊信號為例，沖擊特征頻率為50Hz，定義代數連通度能量占比為路圖傅里葉變換中代數連通度對應階次以內能量與所有階次能量的比值

比較在不同信噪比下（SNR=-5，0，5），零特征值數目尾、代數連通度能量占比r及重構信號包絡譜在特征頻率處的幅值隨熱核寬度（0.001≤δ≤0.06）變化的情況，結果如圖4所示。

從圖4（a）中可看出，隨著熱核寬度增加，零特征值數目逐漸減少，因此沖擊特征也逐漸集中于較低階次的有限特征向量中;從圖4（b）中可看出，隨著熱核寬度增加，代數連通度能量占比逐漸減小;圖4（c）中特征頻率處幅值存在最大值，此時熱核寬度為最優熱核寬度值。綜合比較，熱核寬度最優時，既要保證零特征值的數目較小又要保證代數連通度能量占比較大;相同沖擊信號在不同信噪比下的最優熱核寬度值不同，且信噪比越小取值越大，對應幅值也隨噪聲增加而減小。因此，針對具體信號最優熱核寬度值是一個變動的值，不同信號取值不同，本文中以重構信號包絡譜在特征頻率處的幅值為適應度函數，結合粒子群優化算法對熱核寬度參數進行尋優，以達到最優的沖擊信號重構效果。

2 基于最優加權路圖GFT的滾動軸承故障診斷

滾動軸承存在局部故障時會產生周期性的瞬態沖擊特征，通常使用包絡解調方法即可診斷出軸承故障。然而，對于滾動軸承的早期故障，振動信號中包含軸承故障信息的沖擊成分能量很小，往往淹沒在機械系統諧波振動與背景噪聲中，此時直接使用包絡解調方法難以取得好的效果。因此，需要對包含軸承故障信息的沖擊成分進行提取。

本文提出基于最優加權路圖GFT沖擊提取的滾動軸承故障診斷方法，該方法將振動信號轉換為路圖，在高斯加權下使用有限拉普拉斯矩陣特征向量重構沖擊信號;以重構沖擊信號故障特征處的包絡譜值為適應度函數，利用粒子群算法尋優最優高斯熱核寬度參數;對最優沖擊分量進行包絡解調診斷軸承故障。故障診斷流程如圖5所示，步驟如下：

（1）采集滾動軸承故障振動信號，采樣點數為N，并將其轉換為N頂點路圖。

（2）使用粒子群優化算法對高斯熱核寬度參數進行尋優。首先確定熱核寬度值的尋優范圍。本文中選取0.001≤δ≤δ_max，δ_max為信號方差;

（3）初始化粒子群。設置粒子群參數維度、粒子群種群大小及最大迭代次數;

（4）計算各個粒子的適應度。首先構建路圖對應的圖拉普拉斯矩陣;然后對拉普拉斯矩陣進行特征分解，通過代數連通度λ_k確定用于沖擊信號重構的特征向量個數k;最后，根據式（9）進行信號重構，并計算其故障特征頻率處對應的包絡譜值作為粒子的適應度值;

（5）計算粒子群中的個體最優值和全局最優值;

（6）更新粒子的速度和位置，達到最大迭代次數時停止更新，返回最優熱核寬度參數值;

（7）根據最優熱核寬度δ，重新按步驟（4）重構沖擊分量，對重構信號進行包絡解調分析，并診斷滾動軸承故障類型。

3 仿真分析

為了驗證本文提出方法的有效性，利用滾動軸承數值仿真模型，模擬滾動軸承局部故障，數值仿真模型表示如下

式中h（t）為沖擊信號分量;T為仿真故障沖擊的時間間隔，則故障特征頻率f=1/T;M為仿真沖擊的個數;B_m為第m個沖擊的幅值;β為阻尼比，f_re表示共振頻率;u（t）為單位階躍函數。s（t）為兩種頻率疊加的諧波信號，頻率分別為f₁，f₂。x（t）為仿真合成信號，n（t）為添加的隨機噪聲。本文中，設置故障特征頻率f=50Hz，采樣頻率f_s=4096Hz，采樣點數取N=4096，添加噪聲后合成信號的信噪比為-10dB，其他參數設置如表1所示，仿真沖擊信號，諧波信號及合成信號的時域波形如圖6所示。

從圖6（a）中可以明顯看出周期瞬態沖擊成分及發生時刻，疊加上周期干擾分量及隨機噪聲后，圖6（c）合成信號中已無法識別出任何瞬態沖擊成分。對仿真合成信號直接作包絡譜分析，結果如圖7所示，受諧波分量的干擾，包絡解調分析方法在對信號進行解調時將兩種余弦頻率之差作為解調結果解出，因此圖中出現了頻率峰值220Hz，而故障特征頻率及其倍頻處無明顯峰值，說明軸承故障沖擊特征已完全淹沒。

應用本文提出的最優加權路圖GFT方法提取沖擊信號并進行解調分析。按照圖5中的算法流程，設置熱核寬度最大尋優值為信號方差δ_max=var（x）=0.75，則尋優范圍為[0.01，0.75]。設置粒子種群大小為10，最大迭代次數為50，適應度隨進化代數曲線如圖8（a）所示，在第19次迭代運算時取得最大值，其后保持穩定。此時計算得到的最優熱核寬度值為δ*=0.279，使用最優熱核寬度值提取沖擊信號并做包絡解調分析，包絡譜圖如圖8（b）所示。從圖中可以能夠清晰地識別仿真信號的故障特征頻率f以及其2，3倍頻，以上結果說明了該方法提取沖擊特征的可行性。

作為對比，分別使用文獻提出的歐氏距離加權下的沖擊特征信號提取方法和譜峭度（SK）方法提取沖擊信號并做包絡譜分析，結果分別如圖9和10所示。從圖9中雖然可以看出1倍及2倍故障特征頻率，但均存在54和103Hz的干擾譜線;此外，其故障特征頻率處幅值也明顯低于本文提出方法。使用SK方法計算快速Kurtogram圖，譜峭度指標最大時的中心頻率為f_c=1184Hz，層階數為6.5.從圖10可以看出，沖擊信號包絡譜僅有1倍故障頻率，且存在其他較大的干擾頻率。以上結果表明，最優加權路圖GFT方法能有效重構周期瞬態沖擊信號。

為了更為客觀地對比不同方法提取滾動軸承故障信號中沖擊特征的性能，引入包絡譜故障特征比率指標R_f對本文提出方法、歐氏距離加權GFT方法和SK方法提取的沖擊信號進行評價。該指標表示為

式中 5為信號Hilbert包絡譜;S（f），S（2f）和S（3f）分別為故障特征頻率及其2，3倍頻譜線的幅值。該指標直觀地反映了故障特征頻率在包絡譜中的幅值比率。使用三種方法提取的沖擊特征信號的R_f值如表2所示，易知本文提出方法提取的沖擊特征信號對應的R_f值取得最大值，因此更為明確地診斷了故障。

4 應用實例

為驗證基于最優加權路圖GFT方法提取沖擊信號的有效性，分別設置滾動軸承內、外圈故障進行實驗。實驗在軸承振動模擬實驗臺進行，實驗裝置如圖11所示。內、外圈故障均采用激光在6307E型滾動軸承的內圈和外圈上切割寬為0.15mm、深為0.13mm的槽，以模擬內圈和外圈故障。滾動軸承內圈故障實驗時的轉軸轉速為1200r/min，外圈故障實驗時的轉軸轉速為1500r/min.采用加速度傳感器采集軸承振動信號，采樣頻率設置為4096Hz，采樣點數取1000.通過計算得到內圈故障特征頻率為98.8Hz，外圈故障特征頻率為76.5Hz。

圖12（a）為滾動軸承內圈故障振動信號的時域波形圖，雖然能夠看到沖擊特征，但沖擊信號的周期性不明顯。直接對時域信號作包絡解調分析，包絡譜如圖12（b）所示，可以看出，故障特征頻率譜線f_i被其他與軸承內圈故障無關的干擾譜線淹沒，因此不能識別故障類型。

使用本文提出方法對內圈故障振動信號中的沖擊成分進行重構。首先設置熱核寬度最大尋優值為信號方差δ_max=var（x_i）=1107，則尋優范圍取[0.01，1107];然后設置粒子群大小為10，最大迭代次數為50，如圖13（a）所示適應度值在第6次迭代運算后取得最大值，之后保持穩定，此時計算得到的最優熱核寬度值為δ*=519.3;最后使用最優熱核寬度值重新提取沖擊信號并做包絡解調分析，包絡譜圖如圖13（b）所示。從圖中可以明顯看出1倍和2倍軸承內圈故障特征頻率譜線，準確診斷出了內圈故障。

同樣作為對比，使用基于歐氏距離加權的沖擊特征提取方法提取沖擊信號并作包絡解調分析，結果如圖14所示。從圖14中可知工倍及2倍故障特征頻率譜線并不明顯，且存在較多的干擾譜線，故障特征頻率及其倍頻處的幅值也均小于本文提出方法得到的結果。

使用SK方法對內圈故障信號進行濾波后并作包絡譜，首先計算快速Kurtogram圖，譜峭度指標最大時的中心頻率為f_c=1792Hz，層階數為3.5.濾波信號包絡譜如圖15所示，易知圖中存在2倍故障頻率譜線，故障頻率處無明顯峰值。

滾動軸承外圈的故障特征頻率為f_o=76.5Hz，振動信號的時域波形如圖16（a）所示，從圖中看不到明顯沖擊特征。對外圈故障時域信號作解調分析，包絡譜如圖16（b）所示，圖中外圈故障特征頻率譜線f。被周圍的64Hz等干擾譜線淹沒，因此不能夠直接診斷出軸承故障類型。

使用本文提出方法對外圈故障振動信號中的沖擊成分進行提取。首先設置熱核寬度最大尋優值為信號方差δ_max=var（xo）=16.62，則范圍取[0.01，16.62];然后初始化粒子群尋優算法，設置粒子群大小為10，最大迭代次數為50，如圖17（a）所示，適應度值在第35次進化后取得最大值，其后全局最優值保持穩定，此時計算得到的最優熱核寬度值δ*=4.8;最后使用最優熱核寬度值重新提取沖擊信號并做包絡解調分析，包絡譜圖如圖17（b）所示。從圖中可以明顯看出軸承外圈故障特征頻率及其倍頻譜線，準確診斷出了外圈故障。

作為對比，使用基于歐氏距離的沖擊特征信號提取方法和SK方法提取沖擊信號并做包絡譜分析，結果分別如圖18和19所示。從圖18中可以看出，故障特征頻率及其倍頻譜線f。和2f。并不明顯。使用SK方法，譜峭度指標最大時的中心頻率為fo=341Hz，層階數為3.5.濾波信號包絡譜如圖19所示，圖中僅存在1倍故障頻率譜線。從上述結果可知，針對滾動軸承故障的沖擊特征，本文提出方法的提取結果更為有效。

為了進一步對比三種不同方法提取內、外圈故障信號中沖擊特征的效果，使用R_f指標進行對比分析，結果如表3所示。從表中可知，對于內、外圈故障，本文提出方法提取的沖擊特征信號的R_f均取得最大值，說明了該方法提取沖擊特征的有效性。

5 結論

（1）本文提出了基于最優加權路圖GFT的沖擊分量提取方法，將圖信號處理方法運用于沖擊信號的提取中，并結合包絡解調實現了軸承的故障診斷。該方法首先將振動信號轉換為路圖，在最優高斯加權下，利用有限圖拉普拉斯矩陣特征向量結合GFT方法重構沖擊特征信號，實現了滾動軸承故障的有效診斷。

（2）實驗結果表明，本文提出方法能夠在背景噪聲或其他強故障特征分量干擾下，有效重構故障沖擊特征，并準確診斷滾動軸承故障類型。與歐氏距離加權下的重構結果和譜峭度方法重構結果對比，驗證了最優加權路圖GFT沖擊提取方法的準確性和優越性。