摘要： 針對(duì)最大二階循環(huán)平穩(wěn)盲解卷積（CYCBD）算法在軸承故障特征提取中的有效性及計(jì)算效率受濾波器長(zhǎng)度影響的問題，提出諧波譜峰因子（HSC）作為評(píng)價(jià)指標(biāo)自適應(yīng)確定CYCBD的濾波器長(zhǎng)度，通過編碼器時(shí)間序列重構(gòu)的方法平衡優(yōu)化過程的計(jì)算效率。根據(jù)滾動(dòng)軸承固有參數(shù)計(jì)算軸承故障階次，并根據(jù)其設(shè)置循環(huán)頻率；根據(jù)故障階次確定時(shí)間序列重構(gòu)的脈沖數(shù)；用中心差分法計(jì)算重構(gòu)后信號(hào)的瞬時(shí)角速度；采用等步長(zhǎng)搜索策略以諧波譜峰因子作為評(píng)價(jià)指標(biāo)自適應(yīng)確定CYCBD的濾波器長(zhǎng)度；根據(jù)諧波譜峰因子最大時(shí)對(duì)應(yīng)的階次譜揭示滾動(dòng)軸承故障特征。仿真和試驗(yàn)數(shù)據(jù)分析結(jié)果表明，所提方法能自適應(yīng)確定濾波器長(zhǎng)度，對(duì)提高CYCBD算法計(jì)算效率有明顯效果，適用于滾動(dòng)軸承故障特征提取。

關(guān)鍵詞： 故障診斷； 滾動(dòng)軸承； 編碼器；最大二階循環(huán)平穩(wěn)盲解卷積； 諧波譜峰因子

中圖分類號(hào)： TH165+.3； TH133.33 文獻(xiàn)標(biāo)志碼： A 文章編號(hào)： 1004-4523（2024）09-1616-09

DOI：10.163GnWSZR/CcFPf5dsFEc4jXA==85/j.cnki.issn.1004-4523.2024.09.018

引 言

目前，基于振動(dòng)信號(hào)的滾動(dòng)軸承故障檢測(cè)技術(shù)得到了廣泛應(yīng)用和發(fā)展［1］。然而，部分應(yīng)用中仍受傳輸路徑的影響以及存在不適合傳感器外部安裝等［2］限制。因此，研究基于其他類型傳感器信號(hào)的滾動(dòng)軸承檢測(cè)方法具有重要意義。

編碼器作為運(yùn)動(dòng)和位置的控制與反饋單元被廣泛裝備于高速主軸、伺服電機(jī)、數(shù)控機(jī)床、機(jī)器人、風(fēng)力渦輪機(jī)等旋轉(zhuǎn)機(jī)械系統(tǒng)中［3］。當(dāng)滾動(dòng)軸承發(fā)生局部故障時(shí)，滾動(dòng)體通過滾道的速度較無故障狀態(tài)時(shí)會(huì)發(fā)生突變，即瞬時(shí)角速度（Instantaneous Angular Speed，IAS）發(fā)生規(guī)律性波動(dòng)［4］。因此，可通過提取編碼器信號(hào)估計(jì)出的IAS周期性抖動(dòng)實(shí)現(xiàn)滾動(dòng)軸承故障檢測(cè)。在實(shí)際中，中心差分法（Central Difference Method，CDM）被廣泛應(yīng)用于IAS估計(jì)中，但由于CDM受采樣頻率和測(cè)量誤差干擾［5］，導(dǎo)致從CDM估計(jì)的IAS中難以直接識(shí)別故障引起的較弱抖動(dòng)。此外，由于軸承在旋轉(zhuǎn)機(jī)械系統(tǒng)中主要起支撐作用，對(duì)扭矩的傳遞不敏感，導(dǎo)致故障引起的IAS波動(dòng)往往被齒輪嚙合和隨機(jī)噪聲所淹沒。目前，基于編碼器信號(hào)的軸承故障特征提取方法的研究仍處于起步階段，如何降低估計(jì)IAS的誤差并提取IAS信號(hào)中故障引起的周期性沖擊面臨較大挑戰(zhàn)。

針對(duì)上述基于IAS信號(hào)的滾動(dòng)軸承特征提取問題，本文結(jié)合編碼器信號(hào)的特性提出編碼器時(shí)間序列重構(gòu)的方法以降低CDM估計(jì)IAS的誤差。此外，結(jié)合最大二階循環(huán)平穩(wěn)盲解卷積（Maximum Second?order Cyclostationarity Blind Deconvolution，CYCBD）［6］對(duì)沖擊的恢復(fù)效果，采用其對(duì)IAS信號(hào)進(jìn)行增強(qiáng)，從而實(shí)現(xiàn)軸承故障特征提取。然而，CYCBD的效果取決于循環(huán)頻率和濾波器長(zhǎng)度。循環(huán)頻率越接近實(shí)際故障頻率，脈沖增強(qiáng)效果越好。在一定范圍內(nèi)，濾波器長(zhǎng)度越長(zhǎng)，效果越好，但濾波器太長(zhǎng)會(huì)導(dǎo)致計(jì)算量增大［7?8］，參數(shù)的確定及計(jì)算量限制了CYCBD的應(yīng)用。值得一提的是，與振動(dòng)信號(hào)相比，編碼器信號(hào)無需通過額外安裝轉(zhuǎn)速計(jì)來計(jì)算理論故障頻率，可根據(jù)軸承參數(shù)直接計(jì)算故障階次、設(shè)置循環(huán)頻率。但如何自適應(yīng)確定濾波器長(zhǎng)度仍然是CYCBD應(yīng)用的關(guān)鍵。為了降低計(jì)算量，研究中結(jié)合編碼器信號(hào)特性提出時(shí)間序列重構(gòu)的方法對(duì)原信號(hào)進(jìn)行重采樣，在減小CDM估計(jì)IAS誤差的同時(shí)減小CYCBD的計(jì)算量，用該方法平衡優(yōu)化過程的計(jì)算效率。故障特征提取的目標(biāo)是通過信號(hào)處理的方式增強(qiáng)較弱的故障成分，從頻譜中識(shí)別出故障頻率及倍頻成分，故障頻率及倍頻分量越明顯，說明檢測(cè)效果越好［9］。張龍等［9］用包絡(luò)譜中最大值與均方根值的比值來評(píng)價(jià)滾動(dòng)軸承故障特征頻率在包絡(luò)譜中的占比情況，該值越高，表示故障特征提取的效果越好。根據(jù)該理論，本文提出諧波譜峰因子（Crest of Harmonics Spectrum，HSC）作為自適應(yīng)確定濾波器長(zhǎng)度的依據(jù)，考慮到搜索的范圍與計(jì)算效率，本文采用等步長(zhǎng)搜索策略對(duì)濾波器長(zhǎng)度進(jìn)行選擇。

綜上所述，本文提出了一種基于編碼器時(shí)間序列重構(gòu)和自適應(yīng)確定CYCBD濾波器長(zhǎng)度的滾動(dòng)軸承故障特征提取方法。首先計(jì)算軸承故障階次并根據(jù)其設(shè)置循環(huán)頻率；然后對(duì)編碼器輸出的時(shí)間序列進(jìn)行重構(gòu)，對(duì)重構(gòu)后的編碼器信號(hào)采用CDM估計(jì)IAS；在特定的濾波器長(zhǎng)度范圍內(nèi)等步長(zhǎng)執(zhí)行CYCBD算法；最后通過HSC最大時(shí)對(duì)應(yīng)的階次譜實(shí)現(xiàn)滾動(dòng)軸承故障特征提取。

1 理論基礎(chǔ)

1.1 編碼器信號(hào)的原理及測(cè)量

旋轉(zhuǎn)編碼器根據(jù)輸出信號(hào)的形式通常分成正弦（模擬）和TTL（數(shù)字）兩種。與模擬編碼器相比，數(shù)字編碼器具有容易進(jìn)行數(shù)字處理、控制可靠、價(jià)格低等優(yōu)點(diǎn)，因此被廣泛應(yīng)用于工業(yè)機(jī)器人伺服電機(jī)控制中［10］。光編碼器主要由LED光源、帶有N個(gè)等間距縫隙的光柵盤和光電二極管陣列組成。當(dāng)編碼器隨軸轉(zhuǎn)動(dòng)時(shí)，軸旋轉(zhuǎn)一圈對(duì)應(yīng)輸出N個(gè)周期的方波，通過高數(shù)計(jì)數(shù)器拾取編碼器脈沖，從而獲得一系列離散的時(shí)間序列［t1，t2，…，ti，…，tL］及其對(duì)應(yīng)的角度位置序列［φ1，φ2，…，φi，…，φL］，時(shí)間和位置序列的對(duì)應(yīng)關(guān)系如圖1所示，圖中Δt1，Δt2，…，ΔtL分別表示每個(gè)脈沖對(duì)應(yīng)的時(shí)間，Δθ表示每個(gè)脈沖對(duì)應(yīng)的角度。

由圖1可知，時(shí)間序列ti可表示為：

（1）

式中 fs為高速計(jì)數(shù)器的頻率；li為編碼器起始點(diǎn)（l0）到第i個(gè)脈沖的長(zhǎng)度。

ti對(duì)應(yīng)的位置序列可表示為：

（2）

式中 N為編碼器的線數(shù)。

1.2 瞬時(shí)角速度估計(jì)及其誤差分析

編碼器輸出的時(shí)間序列ti和位置序列φi通常使用CDM估計(jì)IAS［5］，其計(jì)算式為：

（3）

式中 ωi表示第i個(gè)脈沖時(shí)刻對(duì)應(yīng)的IAS。

由式（3）可知，測(cè)量的編碼器脈沖長(zhǎng)度（li+1和li-1）受量化誤差和測(cè)量噪聲的干擾［11］。此外，IAS的精度也受fs的影響。由于編碼器的角度分辨率為2π/N，因此每個(gè)位置序列對(duì)應(yīng)的量化誤差為π/N，根據(jù)式（3）獲得的IAS的測(cè)量誤差為［11］：

（4）

式中 Δl為兩相鄰脈沖之間的長(zhǎng)度；表示編碼器的量化誤差；Δfs為采樣誤差，由于采集時(shí)鐘具有較高的準(zhǔn)確性，Δfs可近似為0。

由式（4）可知，IAS的誤差與轉(zhuǎn)速無關(guān)，與編碼器的線數(shù)、高數(shù)計(jì)數(shù)器的采樣頻率和兩脈沖之間的長(zhǎng)度相關(guān)。然而，實(shí)際中編碼器往往嵌入在設(shè)備中，通過增加其線數(shù)來減小誤差，與現(xiàn)實(shí)不符。此外，采樣頻率至少滿足采樣定律，通過降低采樣頻率難以達(dá)到降低誤差的目的。綜上分析，通過增加實(shí)現(xiàn)減小誤差的目的是較為可行的方式。結(jié)合圖1可知，增加等價(jià)于將多個(gè)脈沖對(duì)應(yīng)的時(shí)間確定為新的時(shí)間序列。

1.3 CYCBD簡(jiǎn)述

CYCBD以最大二階循環(huán)平穩(wěn)指標(biāo)（ICS2）為依據(jù)，采用迭代特征值分解算法求解最優(yōu)逆濾波器h，從而獲得目標(biāo)源信號(hào)so，其目的是從噪聲觀測(cè)信號(hào)x中恢復(fù)源信號(hào)，其過程表示為［6］：

（5）

式中 “*”表示卷積運(yùn)算。

以ICS2為依據(jù)，計(jì)算出最優(yōu)逆濾波器h，ICS2可表示為：

（6）

式中 上標(biāo)“H”表示矩陣的共軛轉(zhuǎn)置；表示噪聲觀測(cè)信號(hào)x的矩陣形式；RXX和RXWX分別表示相關(guān)矩陣和加權(quán)相關(guān)矩陣；W為加權(quán)矩陣，表達(dá)式為：

（7）

式中 P［·］表示觀測(cè)到的信號(hào)s的周期成分；Ls和Nh分別為信號(hào)s和逆濾波器h的長(zhǎng)度；ls表示求和函數(shù)的初值，且初值為。

最優(yōu)逆濾波器h等價(jià)于式（6）中ICS2最大，然后通過求解一個(gè)廣義特征值獲得ICS2的最大值，最大特征值λ可通過下式計(jì)算：

（8）

其中，最大特征值λ對(duì)應(yīng)最優(yōu)的ICS2值。由于加權(quán)矩陣需要通過假設(shè)逆濾波器h來獲得，所以ICS2的最大值需要通過迭代過程實(shí)現(xiàn)。CYCBD算法實(shí)現(xiàn)的具體步驟如下：

（1）利用自回歸模型對(duì)逆濾波器h進(jìn)行初始化，獲得濾波器系數(shù)；

（2）根據(jù)觀測(cè)信號(hào)x和濾波器h計(jì)算加權(quán)相關(guān)矩陣RXWX；

（3）利用式（8）計(jì)算特征值λ及其對(duì)應(yīng)的h；

（4）返回執(zhí)行步驟（2）和（3）直至收斂，從而獲得目標(biāo)源信號(hào)。

更多細(xì)節(jié)可參閱文獻(xiàn)［6］。

2 時(shí)間序列重構(gòu)及濾波器長(zhǎng)度確定

綜上分析，CDM估計(jì)編碼器信號(hào)的IAS時(shí)，故障引起的IAS抖動(dòng)因受誤差干擾難以提取，而采用CYCBD增強(qiáng)故障引起IAS抖動(dòng)的效果取決于濾波器長(zhǎng)度，濾波器長(zhǎng)度過大會(huì)導(dǎo)致計(jì)算量增加和信號(hào)失真［8］。對(duì)此，研究提出編碼器時(shí)間序列重構(gòu)的方法以減小CDM估計(jì)IAS的誤差，同時(shí)平衡自適應(yīng)確定濾波器長(zhǎng)度的計(jì)算成本。

2.1 編碼器時(shí)間序列重構(gòu)

由式（4）可知，增加有利于減小誤差，即增加每個(gè)時(shí)間序列對(duì)應(yīng)的脈沖數(shù)，研究中將該步驟稱為時(shí)間序列重構(gòu)。該步驟能降低中心差分法估計(jì)IAS信號(hào)的誤差，同時(shí)減少每個(gè)周期的數(shù)據(jù)長(zhǎng)度，有利于CYCBD算法增強(qiáng)故障引起的IAS抖動(dòng)及減少算法的計(jì)算量。為了盡量減少CYCBD的計(jì)算量，在滿足IAS信號(hào)分辨率的前提下，重構(gòu)后每個(gè)周期包含的數(shù)據(jù)長(zhǎng)度應(yīng)盡量小。在故障檢測(cè)中，往往通過故障頻率及倍頻在頻譜中的占比情況確定故障是否發(fā)生。因此，可以根據(jù)滾動(dòng)軸承理論故障階次確定重構(gòu)后每一個(gè)周期的采樣點(diǎn)數(shù)。首先根據(jù)軸承參數(shù)計(jì)算各元件的理論故障特征階次，計(jì)算式［12］為：

（9）

（10）

（11）

式中 Oo，Oi和Ob分別表示外圈、內(nèi)圈和滾動(dòng)體的故障階次；n為滾動(dòng)體個(gè)數(shù)；d為滾動(dòng)體直徑；D為節(jié)圓直徑；φ為接觸角。

黃包裕等［13］通過前三階諧波顯著性指標(biāo)判斷故障特征的顯著程度。為了保證重構(gòu)后信號(hào)的階次譜中包含三倍故障階次分量，考慮頻譜中頻率為單邊譜，顯示的頻率為采樣頻率的一半，重構(gòu)后信號(hào)采樣頻率為fs=6Ox，其中Ox表示軸承內(nèi)圈、外圈或滾動(dòng)體的故障階次。因此，若將k個(gè)脈沖作為一個(gè)整體對(duì)時(shí)間序列進(jìn)行重構(gòu)，且每個(gè)周期的數(shù)據(jù)不低于fs可滿足重構(gòu)需求，重構(gòu)后每個(gè)周期的采樣點(diǎn)數(shù)及k可通過下式計(jì)算：

（12）

式中 mod（N，k）表示計(jì)算能整除N的所有k值；f表示N被k整除時(shí)的商；find［·］表示計(jì)算不小于fs的所有因子；arg min｛·｝表示取最小參數(shù)值。

式（12）表示計(jì)算能整除N且不小于fs的最小整數(shù)kopt和重構(gòu)后每個(gè)周期的采樣點(diǎn)數(shù)fopt。根據(jù)式（12）計(jì)算的k值對(duì)原始時(shí)間序列ti進(jìn)行重構(gòu)，其表達(dá)式為：

（13）

式中 表示第i×k個(gè)脈沖對(duì)應(yīng)的時(shí)間，其對(duì)應(yīng)的位置序列為：

（14）

式中 表示重構(gòu)后的位置序列。

2.2 濾波器長(zhǎng)度自適應(yīng)確定

CYCBD算法的效果主要取決于濾波器長(zhǎng)度，在保證計(jì)算效率的前提下，自適應(yīng)確定濾波器長(zhǎng)度對(duì)拓展CYCBD在故障特征提取方面的應(yīng)用具有重要作用。自適應(yīng)確定濾波器長(zhǎng)度的核心在于選擇合適的指標(biāo)評(píng)價(jià)濾波效果。在階次譜中故障頻率及倍頻分量越明顯，說明檢測(cè)效果越好。根據(jù)文獻(xiàn)［9］提出的包絡(luò)譜譜峰因子的思想，本文提出以HSC作為自適應(yīng)確定濾波器長(zhǎng)度的依據(jù)，其表達(dá)式如下：

（15）

式中 FS表示階次譜；jOx表示理論故障特征階次的第j次倍頻；max（·）和rms（·）分別表示計(jì)算FS在對(duì)應(yīng)區(qū)間內(nèi)的最大值和有效值。

需要注意的是，由于滑移導(dǎo)致理論故障特征階次與實(shí)際故障特征階次存在1%～2%的滑移誤差［13］，為了避免轉(zhuǎn)頻、嚙合頻率的干擾并保證搜索的最大值是感興趣的故障成分，式（15）中每個(gè)倍頻最大值的搜索區(qū)間為［0.98jOx，1.02jOx］。

此外，為了不過多增加計(jì)算量并保證信號(hào)不失真，文獻(xiàn)［8］提出等步長(zhǎng)（步長(zhǎng)為100）搜索策略對(duì)CYCBD的濾波器長(zhǎng)度進(jìn)行搜索，搜索次數(shù)為10次。由于文獻(xiàn)［8］中的故障特征頻率約為本文的10倍，采用等比例縮小的方法確定本文步長(zhǎng)為10。試驗(yàn)中發(fā)現(xiàn)，濾波器長(zhǎng)度的搜索范圍為［fopt∶10∶fopt+90］時(shí)具有較好的效果。此外，文獻(xiàn)［8］說明了循環(huán)頻率是一個(gè)以故障頻率為初始值和間隔，長(zhǎng)度為10倍故障頻率的數(shù)組。因此，本文的循環(huán)頻率為［Ox∶Ox∶10Ox］。

2.3 故障特征提取流程

本文提出方法的計(jì)算流程如圖2所示，其主要實(shí)現(xiàn)步驟如下：

（1）根據(jù)式（9）～（11）計(jì)算理論故障階次，然后通過式（12）計(jì)算出重構(gòu)時(shí)間序列的脈沖數(shù)kopt及濾波器長(zhǎng)度搜索的初始值fopt；

（2）根據(jù)式（13）對(duì)時(shí)間序列進(jìn)行重構(gòu)，根據(jù)式（14）計(jì)算出重構(gòu)后的位置序列，然后求IAS信號(hào)；

（3）根據(jù)理論故障特征階次設(shè)置循環(huán)頻率［Ox∶Ox∶10Ox］，根據(jù)fopt確定濾波器長(zhǎng)度搜索范圍［fopt∶10∶fopt+90］，迭代次數(shù)和使算法收斂的相對(duì)誤差值對(duì)算法效果影響不大，根據(jù)文獻(xiàn)［6］的推薦值分別選為30和10-3，然后執(zhí)行CYCBD計(jì)算，并根據(jù)式（15）求解卷積后信號(hào)的HSC；

（4）根據(jù)濾波器設(shè)定的范圍取不同的值分別執(zhí)行步驟（3）獲得每個(gè)濾波器長(zhǎng)度對(duì)應(yīng)的HSC，選擇HSC最大時(shí)對(duì)應(yīng)的階次譜揭示滾動(dòng)軸承故障。

3 仿真分析

為了驗(yàn)證所提方法的有效性，首先采用仿真信號(hào)對(duì)算法進(jìn)行驗(yàn)證。文獻(xiàn)［5］將故障引起的IAS抖動(dòng)簡(jiǎn)化成一種正弦變化模型，IAS信號(hào)仿真模型如下：

（16）

式中 v0為恒定轉(zhuǎn)速；x1（t）表示軸承外圈故障引起的IAS變化；x2（t）表示齒輪嚙合干擾成分；n（t）為測(cè)量噪聲和非平穩(wěn)負(fù)載變化引起的速度變化；S1和S2分別表示軸承故障引起的IAS波動(dòng)幅值和嚙合干擾的幅值；β1和β2為衰減系數(shù)；fr為轉(zhuǎn)頻；T1和T2分別表示軸承故障周期和齒輪嚙合周期；τ為隨機(jī)滑移的時(shí)間。仿真參數(shù)如表1所示。

仿真中，故障周期T1對(duì)應(yīng)的軸承故障階次為5.4×，編碼器的線數(shù)為5000，隨機(jī)噪聲的幅值為0.5。根據(jù)上述仿真模型和參數(shù)獲得的IAS信號(hào)及其階次譜如圖3所示。

從圖3（b）中可較明顯地看出仿真的嚙合階次（1）及其倍頻，但軸承外圈故障階次（5.4×）及其倍頻難以識(shí)別。將仿真的故障階次（5.4×）和編碼器線數(shù)5000代入式（12）中，計(jì)算得到重構(gòu)后每個(gè)周期的采樣點(diǎn)數(shù)fopt為40，即以125個(gè)原始脈沖對(duì)應(yīng)的時(shí)間作為一個(gè)板塊重構(gòu)時(shí)間序列，重構(gòu)后的IAS信號(hào)及其階次譜如圖4所示。對(duì)比圖4（b）和圖3（b），重構(gòu)后信號(hào)的頻率成分與原信號(hào)相似，說明重構(gòu)后的信號(hào)在降低每個(gè)周期數(shù)據(jù)長(zhǎng)度的同時(shí)有效保留了故障信息。

由于齒輪嚙合頻率具有嚴(yán)格的周期性，本文采用時(shí)域同步平均技術(shù)獲取嚴(yán)格的周期信號(hào)，然后用原信號(hào)減去周期分量獲得包含軸承信號(hào)的分量，再對(duì)殘余信號(hào)進(jìn)行CYCBD處理。根據(jù)仿真的故障階次5.4×，設(shè)置循環(huán)頻率為［5.4∶5.4∶54］，濾波長(zhǎng)度的搜索范圍為［40∶10∶130］，濾波器長(zhǎng)度與HSC指標(biāo)的變化關(guān)系如圖5所示。由圖5可知，當(dāng)濾波器長(zhǎng)度為50時(shí)，指標(biāo)值最大，說明此時(shí)CYCBD算法沖擊增強(qiáng)效果較好，對(duì)應(yīng)的階次譜如圖6所示。

為進(jìn)一步驗(yàn)證本文算法的有效性，將直接用CYCBD算法處理的效果和VMD算法處理的效果進(jìn)行對(duì)比。直接用CYCBD增強(qiáng)信號(hào)后的階次譜如圖7所示。本文采用文獻(xiàn)［14］所提的觀察中心頻率法確定VMD算法的模態(tài)數(shù)為6，分別求出各IMF分量的階次譜，結(jié)果如圖8所示。

對(duì)比圖6，7和8，圖8中IMF1分量的轉(zhuǎn)頻階次及倍頻明顯，但故障階次無法識(shí)別，IMF4分量的階次譜中雖然可以識(shí)別到故障階次及倍頻，但受噪聲干擾嚴(yán)重。所提方法和直接用CYCBD處理后的階次譜中均能較好地識(shí)別出軸承故障階次，但兩種方法的HSC值分別為22.4和16.7，因此，說明所提方法具有更好的故障特征提取效果。

4 試驗(yàn)分析

以滾動(dòng)軸承（NU206）作為研究對(duì)象，軸承參數(shù)如下：節(jié)圓直徑D為46 mm；滾動(dòng)體直徑d為9 mm；滾動(dòng)體個(gè)數(shù)n為13；接觸角φ為0°。采用線切割加工技術(shù)在軸承內(nèi)圈、外圈和滾動(dòng)體上分別加工尺寸為0.5 mm的故障，滾動(dòng)軸承模擬故障如圖9（a）所示，將上述故障分別在圖9（b）所示的試驗(yàn)平臺(tái)上進(jìn)行試驗(yàn)，試驗(yàn)中編碼器的線數(shù)N為5000。

編碼器的局部脈沖信號(hào)如圖10（a）所示，通過式（1）～（2）計(jì)算出編碼器的時(shí)間和位置序列如圖10（b）所示。

根據(jù)式（3）估計(jì)編碼器的IAS信號(hào)，并計(jì)算其階次譜，結(jié)果如圖11所示。

圖11的階次譜中軸承故障階次的識(shí)別受干擾嚴(yán)重，難以提取特征。用本文所提方法對(duì)原始時(shí)間序列進(jìn)行重構(gòu)，將軸承參數(shù)代入式（9）中計(jì)算出軸承外圈故障階次為5.23×，通過式（12）計(jì)算出重構(gòu)時(shí)間序列的脈沖數(shù)為125，重構(gòu)后每周的長(zhǎng)度為40，重構(gòu)后的IAS信號(hào)及階次譜如圖12所示。對(duì)比圖12（b）和圖11（b）可知，時(shí)間序列重構(gòu)后信號(hào)的頻率成分與原信號(hào)基本一致。齒輪的嚙合階次及其倍頻明顯，而故障階次難以識(shí)別。為降低齒輪嚙合干擾，對(duì)信號(hào)進(jìn)行同步平均處理獲得周期性成分，然后用原信號(hào)減去平均后的信號(hào)，從而獲得包含軸承信息豐富的信號(hào)，再采用本文提出的CYCBD算法自適應(yīng)確定濾波器長(zhǎng)度，濾波器長(zhǎng)度與HSC指標(biāo)的變化關(guān)系如圖13所示。由圖13可知，當(dāng)濾波器長(zhǎng)度為110時(shí)，恢復(fù)沖擊的效果較好，處理后信號(hào)的階次譜如圖14所示。

直接用CYCBD算法對(duì)外圈故障的IAS信號(hào)進(jìn)行處理，結(jié)果如圖15所示。用VMD對(duì)信號(hào)進(jìn)行處理，根據(jù)觀察中心頻率法選擇模態(tài)數(shù)為4，分別求出各IMF分量的階次譜如圖16所示。

對(duì)比圖14，15和16，本文所提方法的階次譜中故障特征階次及倍頻較明顯，說明特征提取效果更好。為了進(jìn)一步說明本文算法的有效性，采用內(nèi)圈和滾動(dòng)體的故障試驗(yàn)進(jìn)行驗(yàn)證。內(nèi)圈故障試驗(yàn)的IAS信號(hào)及其階次譜如圖17所示。

用本文方法對(duì)內(nèi)圈故障信號(hào)進(jìn)行處理，處理后信號(hào)的階次譜如圖18所示。用VMD對(duì)信號(hào)進(jìn)行處理，選擇模態(tài)數(shù)為3，分別求出各IMF分量的階次譜如圖19所示。

與圖17（b）和圖19相比，圖18階次譜中內(nèi)圈故障階次（7.7×）及其倍頻較明顯，說明所提方法能有效提取軸承內(nèi)圈故障特征。

滾動(dòng)體故障試驗(yàn)的IAS信號(hào)及其階次譜如圖20所示。用本文方法處理的結(jié)果如圖21所示，用VMD方法處理的結(jié)果如圖22所示。與圖20（b）和圖22相比，圖21中能較好地識(shí)別出滾動(dòng)體故障階次（9.8×）及其倍頻，說明采用所提方法提取滾動(dòng)體故障特征的有效性。

為驗(yàn)證所提方法在計(jì)算量方面的優(yōu)勢(shì)，選用CPU為Intel Core i7?6700、內(nèi)存為8 G的計(jì)算機(jī)，取4個(gè)周期的數(shù)據(jù)（數(shù)據(jù)長(zhǎng)度為20000），對(duì)比所提方法、直接用CYCBD處理（循環(huán)頻率和濾波器長(zhǎng)度與所提方法確定的參數(shù)一致）和使用VMD算法的計(jì)算時(shí)間。上述3種處理方式的計(jì)算效率對(duì)比結(jié)果如表2所示。

由表2可知，VMD算法的計(jì)算效率與所提方法相近，但所提方法特征提取效果更好。與直接用CYCBD算法的計(jì)算效率對(duì)比，所提方法具有明顯優(yōu)勢(shì)。值得一提的是，表2中所提方法的計(jì)算時(shí)間包括濾波器長(zhǎng)度優(yōu)化過程的計(jì)算量，而直接用CYCBD的計(jì)算時(shí)間僅為所提方法對(duì)應(yīng)參數(shù)的執(zhí)行時(shí)間，如考慮濾波器優(yōu)化過程，直接用CYCBD的方式計(jì)算量會(huì)成倍增加。因此，所提時(shí)間序列重構(gòu)的方法對(duì)提高CYCBD的計(jì)算效率具有重要作用。

5 結(jié) 論

仿真和試驗(yàn)分析驗(yàn)證了所提方法的有效性，且具有以下優(yōu)勢(shì)：

（1）與振動(dòng)信號(hào)相比，采用CYCBD算法處理編碼器信號(hào)無需額外安裝轉(zhuǎn)速計(jì)來計(jì)算循環(huán)頻率。

（2）時(shí)間序列重構(gòu)在保留故障信息的基礎(chǔ)上降低中心差分法的估計(jì)誤差；重構(gòu)后信號(hào)每個(gè)周期的數(shù)據(jù)長(zhǎng)度減小，有利于提高CYCBD的計(jì)算效率，促進(jìn)CYCBD算法的應(yīng)用。

（3）通過HSC指標(biāo)實(shí)現(xiàn)CYCBD濾波器長(zhǎng)度自適應(yīng)確定，提高了算法的魯棒性。

參考文獻(xiàn)：

［1］李志農(nóng)， 劉躍凡， 胡志峰， 等. 經(jīng)驗(yàn)小波變換-同步提取及其在滾動(dòng)軸承故障診斷中的應(yīng)用［J］. 振動(dòng)工程學(xué)報(bào)， 2021， 34（6）： 1284-1292.

Li Zhinong， Liu Yuefan， Hu Zhifeng， et al. Empirical wavelet transform-synchroextracting transform and its applications in fault diagnosis of rolling bearing［J］. Journal of Vibration Engineering， 2021， 34（6）： 1284-1292.

［2］Ding C C， Zhao M， Lin J， et al. Sparsity-based algorithm for condition assessment of rotating machinery using internal encoder data［J］. IEEE Transactions on Industrial Electronics， 2020， 67（9）： 7982-7993.

［3］Miao Y H， Zhao M， Yi Y G， et al. Application of sparsity-oriented VMD for gearbox fault diagnosis based on built-in encoder information［J］. ISA Transactions， 2020， 99： 496-504.

［4］Renaudin L， Bonnardot F， Musy O， et al. Natural roller bearing fault detection by angular measurement of true instantaneous angular speed［J］. Mechanical Systems and Signal Processing， 2010， 24（7）： 1998-2011.

［5］Zhao M， Jia X D， Lin J， et al. Instantaneous speed jitter detection via encoder signal and its application for the diagnosis of planetary gearbox［J］. Mechanical Systems and Signal Processing， 2018， 98（1）： 16-31.

［6］Buzzoni M， Antoni J， D’Elia G. Blind deconvolution based on cyclostationarity maximization and its application to fault identification［J］. Journal of Sound and Vibration， 2018， 432： 569-601.

［7］朱丹宸， 張永祥， 何偉， 等. 基于自適應(yīng)CYCBD和互相關(guān)譜的滾動(dòng)軸承復(fù)合故障診斷方法［J］. 振動(dòng)與沖擊， 2020， 39（11）： 116-122.

Zhu Danchen， Zhang Yongxiang， He Wei， et al. Compound faults diagnosis of rolling element bearing using adaptive CYCBD and cross-correlation spectrum［J］. Journal of Vibration and Shock， 2020， 39（11）： 116-122.

［8］Wang Z J， Zhou J， Du W H， et al. Bearing fault diagnosis method based on adaptive maximum cyclostationarity blind deconvolution［J］. Mechanical Systems and Signal Processing， 2022， 162： 108018.

［9］張龍， 成俊良， 李興林， 等. 基于自適應(yīng)頻帶沖擊強(qiáng)度的滾動(dòng)軸承故障量化評(píng)估［J］. 振動(dòng)與沖擊， 2018， 37（19）： 30-38.

Zhang Long， Cheng Junliang， Li Xinglin， et al. Fault quantitative evaluation of rolling bearings based on shock value of selected frequency band［J］. Journal of Vibration and Shock， 2018， 37（19）： 30-38.

［10］Zhao M， Lin J. Health assessment of rotating machinery using a rotary encoder［J］. IEEE Transactions on Industrial Electronics， 2018， 65（3）： 2548-2556.

［11］趙明， 林京， 王琇峰， 等. 基于增量式編碼器的信息的機(jī)床振動(dòng)溯源技術(shù)［C］∥第九屆全國(guó)設(shè)備與維修工程學(xué)術(shù)會(huì)議暨第十五屆全國(guó)設(shè)備監(jiān)測(cè)與診斷學(xué)術(shù)會(huì)議. 北京， 2012： 610-615.

Zhao Ming， Lin Jing， Wang Xiufeng， et al. Machine tool vibration tracing technology based on incremental encoder information［C］∥ The 9th National Equipment and Maintenance Engineering Academic Conference and the 15th National Equipment Monitoring and Diagnosis Academic Conference. Beijing， 2012： 610-615.

［12］Wang D， Tse P W， Tsui K L. An enhanced Kurtogram method for fault diagnosis of rolling element bearings［J］. Mechanical Systems and Signal Processing， 2013， 35（1-2）： 176-199.

［13］黃包裕， 張永祥， 趙磊. 基于布谷鳥搜索算法和最大二階循環(huán)平穩(wěn)盲解卷積的滾動(dòng)軸承故障診斷方法［J］. 機(jī)械工程學(xué)報(bào)， 2021， 57（9）： 99-107.

Huang Baoyu， Zhang Yongxiang， Zhao Lei. Research on fault diagnosis method of rolling bearings based on Cuckoo search algorithm and maximum second order cyclostationary blind deconvolution［J］. Journal of Mechanical Engineering， 2021， 57（9）： 99-107.

［14］劉長(zhǎng)良， 武英杰， 甄成剛. 基于變分模態(tài)分解和模糊C均值聚類的滾動(dòng)軸承故障診斷［J］. 中國(guó)電機(jī)工程學(xué)報(bào)， 2015， 35（13）： 3358-3365.

Liu Changliang， Wu Yingjie， Zhen Chenggang. Rolling bearing fault diagnosis based on variational mode decomposition and fuzzy C means clustering［J］. Proceedings of the CSEE， 2015， 35（13）： 3358-3365.

Application of encoder time series reconstruction and CYCBD in fault feature extraction of rolling bearing

YANG Xin-min， GUO Yu， CHEN Xin， FAN Jia-wei

（Faculty of Mechanical and Electrical Engineering， Kunming University of Science and Technology， Kunming 650500， China）

Abstract： Aiming at the problem that the effectiveness and computation cost of the maximum second-order cyclic stationarity blind deconvolution （CYCBD） algorithm are affected by parameter setting，the crest of the harmonics spectrum （HSC） was used as an evaluation index to determine the length of the CYCBD filter adaptively， and the calculation cost of the optimization process was balanced by the method of encoder time series reconstruction. The bearing fault order is calculated and the cycle frequency is set according to it. The pulse number of time series reconstruction is determined according to the fault order. The central difference method （CDM） is used to calculate the instantaneous angular speed （IAS） of the reconstructed signal. The filter length of CYCBD was adaptively selected with the HSC as the evaluation index using the equal-step search strategy. The spectrum corresponding to the maximum HSC is calculated to achieve fault feature extraction. The simulation and experimental data analysis results show that the proposed method can adaptively select the filter length， which has an obvious effect on reducing the cost of the CYCBD algorithm， and is effective for rolling bearing fault feature extraction.

Key words： fault diagnosis； rolling bearing； encoder；maximum second-order cyclostationarity blind deconvolution；crest ofharmonics spectrum

作者簡(jiǎn)介： 楊新敏（1994―），男，博士，講師，碩士生導(dǎo)師。E-mail： yxm1434@163.com。

通訊作者： 郭 瑜（1971―），男，博士，教授，博士生導(dǎo)師。E-mail： kmgary@163.com。