摘要： 針對(duì)強(qiáng)噪聲下難以估計(jì)信號(hào)瞬時(shí)頻率的問題，提出了基于質(zhì)心動(dòng)態(tài)路徑規(guī)劃（DPPB）的變轉(zhuǎn)速齒輪箱時(shí)頻脊線索引算法。該算法在剖析多路徑匹配追蹤（MMP）脊線索引算法及其在強(qiáng)噪聲下失效原因的基礎(chǔ)上，通過對(duì)MMP算法得到的脊線集加窗求質(zhì)心，構(gòu)建信號(hào)的脊線質(zhì)心稀疏矩陣，并針對(duì)質(zhì)心稀疏矩陣設(shè)計(jì)動(dòng)態(tài)路徑規(guī)劃函數(shù)索引脊線上的質(zhì)心點(diǎn)，根據(jù)脊線代價(jià)函數(shù)值計(jì)算最優(yōu)時(shí)頻脊線。將相似度系數(shù)Ra和置信度σRa作為脊線提取效果的衡量指標(biāo)，通過仿真和試驗(yàn)驗(yàn)證了DPPB算法可有效提取強(qiáng)噪聲信號(hào)的時(shí)頻脊線，且不同程度噪聲下的可靠性和魯棒性均優(yōu)于峰值索引算法和MMP脊線索引算法。

關(guān)鍵詞： 故障診斷； 齒輪箱； 質(zhì)心動(dòng)態(tài)路徑規(guī)劃（DPPB）； 時(shí)頻脊線索引； 多路徑匹配追蹤（MMP）； 時(shí)頻分析

中圖分類號(hào)： TH165+.3； TH132.41 文獻(xiàn)標(biāo)志碼： A 文章編號(hào)： 1004-4523（2024）06-1077-12

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

引 言

齒輪箱在變轉(zhuǎn)速工況下會(huì)產(chǎn)生與轉(zhuǎn)速相關(guān)的非平穩(wěn)振動(dòng)信號(hào)，使用傳統(tǒng)的平穩(wěn)信號(hào)分析方法難以診斷齒輪箱故障［1］。瞬時(shí)轉(zhuǎn)速不僅是反映齒輪箱運(yùn)行狀態(tài)的重要參數(shù)，也是階次跟蹤算法、廣義解調(diào)分析等非平穩(wěn)信號(hào)分析方法的必需參量［2?4］。所以準(zhǔn)確獲得變轉(zhuǎn)速齒輪箱瞬時(shí)轉(zhuǎn)速對(duì)齒輪箱健康狀態(tài)監(jiān)測(cè)和故障診斷具有重要意義。

利用高精度轉(zhuǎn)速測(cè)量設(shè)備可獲得振動(dòng)信號(hào)對(duì)應(yīng)的轉(zhuǎn)速，但是由于實(shí)際工作環(huán)境的制約，大部分齒輪箱不具備安裝轉(zhuǎn)速測(cè)量設(shè)備的條件。齒輪箱振動(dòng)信號(hào)中的特征頻率通常與轉(zhuǎn)速相關(guān)，所以通過估計(jì)瞬時(shí)頻率（Instantaneous Frequency，IF）可達(dá)到求取瞬時(shí)轉(zhuǎn)速的目的［5］。郭瑜等［6］針對(duì)變轉(zhuǎn)速旋轉(zhuǎn)機(jī)械，利用峰值索引算法在振動(dòng)信號(hào)時(shí)頻分布中獲得一階轉(zhuǎn)速對(duì)應(yīng)的IF時(shí)頻脊線，實(shí)現(xiàn)了無轉(zhuǎn)速計(jì)條件下估計(jì)振動(dòng)信號(hào)的轉(zhuǎn)頻信息。峰值索引算法計(jì)算簡(jiǎn)單但對(duì)噪聲較為敏感，在實(shí)際工程中往往難以得到滿足要求的時(shí)頻曲線。為了提高估計(jì)算法對(duì)噪聲的魯棒性， latsenko等［7］以頻率平均值和標(biāo)準(zhǔn)差構(gòu)造自適應(yīng)代價(jià)函數(shù)代替時(shí)頻峰值進(jìn)行脊線索引，引入跳頻懲罰保證脊線的連續(xù)性，一定程度上降低了噪聲對(duì)脊線索引的影響。江星星等［8］提出了一種雙向搜索時(shí)頻脊融合算法，運(yùn)用代價(jià)函數(shù)脊線索引方法從正反兩個(gè)方希搜索同一IF分量的時(shí)頻脊線，通過脊線融合算法合并最優(yōu)脊線段從而消除噪聲引起的局部誤差。Ding等［9］則根據(jù)軸承振動(dòng)信號(hào)的頻率特征，運(yùn)用概率密度函數(shù)將基頻與其對(duì)應(yīng)的倍頻進(jìn)行融合，從而達(dá)到修正時(shí)頻脊線、提高轉(zhuǎn)速估計(jì)精度的目的。由于受到Heisenberg不確定原理的影響，時(shí)頻分布中頻率發(fā)散、多分量頻率信號(hào)混疊等現(xiàn)象嚴(yán)重影響時(shí)頻脊線的索引結(jié)果。張炎等［10］運(yùn)用同步壓縮小波變換（Synchrosqueezing Wavelet Transform，SWT）對(duì)軸承振動(dòng)信號(hào)的包絡(luò)譜進(jìn)行時(shí)頻重排，通過提高時(shí)頻包絡(luò)譜的分辨率和聚集性，減小頻譜發(fā)散對(duì)脊線提取的影響，從而準(zhǔn)確提取軸承的特征頻率曲線。上述算法一定程度上提高了時(shí)頻脊線索引的精度，但對(duì)索引初始點(diǎn)的依賴性很高，若選擇錯(cuò)誤的索引初始點(diǎn)便無法提取正確的時(shí)頻脊線。多路徑匹配追蹤（Multipath Matching Pursuit，MMP）［11］通過建立基于樹的路徑搜索模型，運(yùn)用貪心算法搜尋最優(yōu)路徑，降低初始搜索點(diǎn)選擇對(duì)搜索結(jié)果的影響。相繼有學(xué)者利用同步壓縮變換（Synchro?Squeezing Transform， SST）［12］、Multisynchro?squeezing Transform（MSST）［13］、同步提取變換（Synchro?Extracting Transform，SET）［14］、Synchro?Reassigning Transform（SRT）［15］算法獲得變轉(zhuǎn)速齒輪箱振動(dòng)信號(hào)的重排時(shí)頻分布，運(yùn)用MMP算法提取到了較為準(zhǔn)確的瞬時(shí)頻率。然而對(duì)于強(qiáng)噪聲信號(hào)，重排時(shí)頻分布會(huì)出現(xiàn)大量噪聲分量，MMP算法在干擾下依然難以識(shí)別正確的時(shí)頻脊線。

目前脊線索引算法均是以時(shí)頻分布上的離散時(shí)頻點(diǎn)作為索引對(duì)象，強(qiáng)噪聲產(chǎn)生的局部幅值極值會(huì)對(duì)脊線索引造成很強(qiáng)的干擾；同時(shí)現(xiàn)有算法將局部最優(yōu)點(diǎn)的集合作為最終的脊線，所以其結(jié)果不一定滿足全局最優(yōu)。針對(duì)現(xiàn)有算法的不足，本文提出一種質(zhì)心動(dòng)態(tài)路徑規(guī)劃（Dynamic Path Planning for Barycenter， DPPB）時(shí)頻脊線索引算法。算法在MMP算法的基礎(chǔ)上構(gòu)建脊線質(zhì)心稀疏矩陣，將索引的對(duì)象由時(shí)頻離散點(diǎn)轉(zhuǎn)換為質(zhì)心點(diǎn)，減少噪聲對(duì)脊線索引的干擾并降低時(shí)頻分布的維度；運(yùn)用動(dòng)態(tài)路徑規(guī)劃算法索引時(shí)頻脊線，盡可能使索引結(jié)果達(dá)到全局最優(yōu)。

1 MMP脊線索引算法

1.1 MMP算法原理及實(shí)現(xiàn)

對(duì)于信號(hào)，窗函數(shù)，通過滑動(dòng)窗口可得到不同時(shí)間段信號(hào)的頻譜圖，從而獲得信號(hào)的時(shí)頻分布如下式所示：

（1）

式中 和表示時(shí)間；表示頻率；表示g的共軛。

根據(jù)時(shí)頻曲線附近局部能量突出且具有一定連續(xù)性的特點(diǎn)，建立最優(yōu)脊線路徑代價(jià)函數(shù)［16］：

（2）

式中 表示幅值對(duì)代價(jià)函數(shù)的權(quán)重系數(shù)；表示第k階導(dǎo)數(shù)對(duì)代價(jià)函數(shù)的權(quán)重系數(shù)；k為求導(dǎo)階數(shù)。

對(duì)式（2）進(jìn)行離散變換，忽略高階項(xiàng)的影響，同時(shí)設(shè)置脊線搜索頻率半徑為，得到如下式所示的脊線索引代價(jià)函數(shù)：

（3）

式中 為權(quán)重系數(shù)，用于確定幅值（時(shí)頻分布）和頻率變化量對(duì)代價(jià)函數(shù)的貢獻(xiàn)率。由于噪聲的影響，時(shí)頻分布中可能出現(xiàn)噪聲時(shí)頻點(diǎn)幅值高于脊線上的幅值，若索引初始點(diǎn)落在噪聲點(diǎn)上則無法索引到正確脊線。

MMP脊線索引算法通過選擇多個(gè)索引初始點(diǎn)的方法消除索引結(jié)果過于依賴初始點(diǎn)的問題。對(duì)于時(shí)頻分布構(gòu)建基于樹的脊線索引模型，選取n個(gè)起始點(diǎn)作為路徑搜索樹的支集起點(diǎn)，以代價(jià)函數(shù)C更新支集原子，通過貪心算法搜索每個(gè)支集上的原子，得到個(gè)支集作為時(shí)頻脊線，然后根據(jù)下式確定最優(yōu)時(shí)頻脊線：

（4）

式中 為離散信號(hào)總長(zhǎng)度。

MMP脊線索引算法具體步驟如下：

（1） 對(duì)于時(shí)頻矩陣，取個(gè)索引初始點(diǎn)將時(shí)間平均分為n段，每一段的時(shí)間間隔為，每個(gè)初始點(diǎn)對(duì)應(yīng)的時(shí)刻為，初始化；

（2） 令，找到第個(gè)初始時(shí)刻的最大幅值點(diǎn)作為脊線搜索的起始點(diǎn)；

（3） 正向搜索：令，根據(jù)式（3）以半徑搜索處的脊線點(diǎn)；

（4） 重復(fù)步驟（3），直到停止迭代，；

（5） 負(fù)向搜索：令，根據(jù)式（3）以半徑搜索處的脊線所在位置，；

（6） 重復(fù)步驟（5） ，直到停止迭代，得到時(shí)頻脊線；

（7） 重復(fù)步驟（2）～（6） ，直到停止迭代；根據(jù)式（4）得到最終的脊線。

1.2 強(qiáng)噪聲下MMP算法失效機(jī)理分析

圖1（a）為信噪比為-10 dB的強(qiáng)噪聲下利用MMP算法得到的時(shí)頻脊線，其中黑色虛線代表理論時(shí)頻曲線，紅色實(shí)線為MMP算法提取的結(jié)果。從圖1（a）中可以看出，時(shí)間在0～7.5 s時(shí)與基本重合，提取效果較好；在7.5 s時(shí)與發(fā)生偏離，然后便沿著偏離后的結(jié)果繼續(xù)搜尋，致使7.5 s后未能提取到有效的時(shí)頻脊線。

圖1（b）為圖1（a）中位置A的局部放大圖，紅色虛線表示式（3）中的搜索頻率半徑。在時(shí)刻，時(shí)頻點(diǎn)a雖然在附近，但在噪聲的影響下幅值并不是搜索半徑內(nèi)的最大值，存在一個(gè)遠(yuǎn)離理論曲線的時(shí)頻點(diǎn)b，其幅值大于，所以在時(shí)刻會(huì)將b作為上的點(diǎn)。當(dāng)脊線索引至?xí)r刻，附近的時(shí)頻點(diǎn)c的幅值高于遠(yuǎn)離的d處幅值，然而此刻c點(diǎn)已經(jīng)超出的范圍，算法只能在的范圍內(nèi)選擇最大幅值點(diǎn)d作為f上的點(diǎn)而無法索引到正確的時(shí)頻點(diǎn)c，使逐步偏離正確的索引方向。

由于MMP算法在索引處的脊線點(diǎn)時(shí)僅根據(jù)確定的索引范圍，得到關(guān)于的局部最優(yōu)值；若在噪聲的影響下偏離，的索引范圍也隨之發(fā)生偏離；當(dāng)?shù)钠屏吭谒饕秶鷥?nèi)，索引結(jié)果僅在附近波動(dòng)；當(dāng)?shù)睦塾?jì)偏移量超出索引范圍，偏離了并且偏離量可能越來越大，得到如圖1（a）中的結(jié)果。如果根據(jù)已搜索到的所有脊線點(diǎn)確定的索引范圍，可以有效抑制局部最優(yōu)導(dǎo)致的脊線偏移，但是在高采樣頻率下數(shù)據(jù)點(diǎn)多，影響計(jì)算效率。所以本文通過構(gòu)造脊線質(zhì)心稀疏矩陣重構(gòu)時(shí)頻分布，減少噪聲干擾并降低脊線索引的數(shù)據(jù)量，進(jìn)而利用動(dòng)態(tài)路徑規(guī)劃算法解決脊線索引時(shí)容易陷入局部最優(yōu)的問題。

2 DPPB時(shí)頻脊線索引算法

2.1 構(gòu)建脊線質(zhì)心稀疏矩陣

由第1節(jié)可知，MMP算法得到最終脊線之前會(huì)生成一個(gè)脊線集，如圖2所示。圖2中彩色細(xì)實(shí)線為算法索引的所有脊線，紅色粗實(shí)線為最終的提取結(jié)果，紫色虛線區(qū)域?yàn)槔碚摃r(shí)頻曲線所在的搜索范圍。

雖然噪聲導(dǎo)致脊線集中不一定存在一條脊線能夠較好地描述信號(hào)真實(shí)的時(shí)頻分布，但是對(duì)于，存在和集合，使得：

（5）

即時(shí)頻脊線上存在某一段正好落在內(nèi)，所有滿足式（5）的脊線點(diǎn)集合便構(gòu)成目標(biāo)脊線。

質(zhì)心可以反映質(zhì)點(diǎn)集的整體分布情況，將式（2）中的視為脊線的質(zhì)量，設(shè)窗函數(shù)為，時(shí)窗半徑為，窗口平移步長(zhǎng)為λ，那么第k條脊線的質(zhì)心可由下式求得：

（6）

針對(duì)離散的時(shí)頻分布，將式（6）改寫為離散形式：

（7）

式中 和分別為質(zhì)心對(duì)應(yīng)的時(shí)間和頻率；為質(zhì)心質(zhì)量；m為時(shí)間序列；為質(zhì)心序列，其中，表示向上取整運(yùn)算，為采樣頻率。

截取一段脊線加窗求質(zhì)心的結(jié)果如圖3所示，其中不同顏色代表加窗后的脊線段，“○”為該段脊線的質(zhì)心，藍(lán)色虛線表示理論曲線。從圖3中可以看出，加窗質(zhì)心運(yùn)算弱化了異常脊線點(diǎn)對(duì)時(shí)頻分布的影響，從而達(dá)到消除噪聲干擾的目的。同時(shí)，質(zhì)心點(diǎn)反映了窗口內(nèi)脊線點(diǎn)的整體分布，所以通過質(zhì)心點(diǎn)描述時(shí)頻脊線可以降低時(shí)頻脊線的維度。

利用式（7）對(duì)圖2中的n段脊線加窗求質(zhì)心，得到一個(gè)由質(zhì)心點(diǎn)構(gòu)成的稀疏矩陣V如圖4所示，圖4（b）為4（a）中區(qū)域B的局部放大圖。圖中“·”為質(zhì)心點(diǎn)的位置，顏色深淺代表該質(zhì)心質(zhì)量的大小。與圖1（a）對(duì)比，圖4中的時(shí)頻脊線更加清晰。

采用二維數(shù)組存儲(chǔ)稀疏矩陣會(huì)浪費(fèi)存儲(chǔ)單元存放零元素，在運(yùn)算中需要花費(fèi)大量時(shí)間對(duì)零元素進(jìn)行無效運(yùn)算。三元組存儲(chǔ)格式（Coordinate，COO）是一種直觀、簡(jiǎn)單的稀疏矩陣存儲(chǔ)格式，分別將二維數(shù)組中非零元素的行、列和數(shù)值存在三個(gè)一維矩陣中，極大地壓縮了原始數(shù)據(jù)量［17］。但是COO存儲(chǔ)數(shù)據(jù)不考慮存儲(chǔ)順序，在數(shù)據(jù)讀取時(shí)需要遍歷整個(gè)矩陣索引需要的數(shù)據(jù)，一定程度上降低運(yùn)算效率。針對(duì)時(shí)頻脊線索引的特點(diǎn)，按照時(shí)間順序運(yùn)用COO方法存儲(chǔ)稀疏矩陣，減少索引時(shí)的數(shù)據(jù)讀取量。首先通過COO方法將b時(shí)刻的質(zhì)心儲(chǔ)存在三個(gè)一維矩陣中，即

（8）

然后將所有時(shí)刻的存儲(chǔ)矩陣組成質(zhì)心稀疏矩陣V：

（9）

通過構(gòu)建質(zhì)心稀疏矩陣將時(shí)頻脊線索引對(duì)象由時(shí)頻分布轉(zhuǎn)化為稀疏矩陣V，保留了時(shí)頻分布中的關(guān)鍵信息，極大地壓縮了數(shù)據(jù)量；在索引某一時(shí)刻的脊線點(diǎn)時(shí)，只需在該時(shí)刻的COO矩陣中查找目標(biāo)元素，從而避免遍歷稀疏矩陣的全部數(shù)據(jù)，達(dá)到降低運(yùn)算量的目的。

2.2 構(gòu)造動(dòng)態(tài)路徑規(guī)劃脊線索引函數(shù)

路徑規(guī)劃可通過現(xiàn)有的數(shù)據(jù)集預(yù)測(cè)接下來一段時(shí)間的數(shù)據(jù)分布，多項(xiàng)式擬合是常用的路徑規(guī)劃方法。但是數(shù)據(jù)量和擬合階數(shù)對(duì)多項(xiàng)式擬合的結(jié)果影響較大，決定了擬合結(jié)果的精度和泛化能力。為了提高脊線點(diǎn)預(yù)測(cè)的魯棒性且盡可能保留脊線局部的細(xì)節(jié)信息，本文采用了動(dòng)態(tài)路徑規(guī)劃的脊線點(diǎn)索引算法。

設(shè)置閾值確定擬合點(diǎn)數(shù)量和擬合階數(shù)，當(dāng)已搜索到脊線的質(zhì)心個(gè)數(shù)時(shí)選用低階數(shù)對(duì)序列的質(zhì)心集進(jìn)行擬合；當(dāng)時(shí)選用高階數(shù)對(duì)序列的質(zhì)心集進(jìn)行擬合；然后計(jì)算序列為的預(yù)測(cè)值，如下式所示：

（10）

式中 ，分別為預(yù)測(cè)點(diǎn)的時(shí)間集和頻率集；；。

定義新的代價(jià)函數(shù)如下式所示：

（11）

式中 為時(shí)刻的預(yù)測(cè)頻率；表示距離。代價(jià)函數(shù)是一個(gè)與距離和質(zhì)量有關(guān)的函數(shù)。

將式（11）得到的代價(jià)函數(shù)值取最大值時(shí)的質(zhì)心點(diǎn)作為的脊線點(diǎn)，并把此時(shí)的作為脊線點(diǎn)的幅值，如下式所示：

（12）

將每條脊線各點(diǎn)的幅值合作為脊線的判別參數(shù)，選擇最大幅值合的脊線作為最終脊線，如下式所示：

（13）

2.3 DPPB脊線索引算法實(shí)現(xiàn)

對(duì)于一段時(shí)間長(zhǎng)度為的信號(hào)，運(yùn)用DPPB算法搜索時(shí)頻脊線的具體步驟如下：

（1） 根據(jù)第1節(jié)的MMP脊線索引算法搜索時(shí)頻分布的時(shí)頻脊線，得到脊線集合。

（2） 設(shè)置窗函數(shù)，時(shí)窗半徑，窗口平移步長(zhǎng)λ，對(duì)脊線集加窗求質(zhì)心，構(gòu)建質(zhì)心稀疏矩陣，計(jì)算初始化，確定擬合階數(shù)閾值m。

（3） 令，搜索中的最大值所對(duì)應(yīng)的質(zhì)點(diǎn)作為第i條脊線的搜尋起始點(diǎn)，即。

（4） 正向搜索：，若，執(zhí)行步驟（5）；否則，令，執(zhí)行步驟（6）。

（5） 由式（10）確定預(yù)測(cè)值，由式（11）計(jì)算時(shí)刻的代價(jià)函數(shù)集，由式（12）搜索脊線點(diǎn)；若，重復(fù)步驟（4）。

（6） 反向搜索：，若，執(zhí)行步驟（7）；否則，執(zhí)行步驟（8）。

（7） 由式（10）確定預(yù)測(cè)值，由式（11）計(jì)算時(shí)刻的代價(jià)函數(shù)集，由式（12）搜索脊線點(diǎn)，重復(fù)步驟（6）。

（8） 令，重復(fù)步驟（3）～（6），直到停止。由式（13）得到最終的時(shí)頻脊線。

DPPB算法流程圖如圖5所示。

2.4 關(guān)于算法的討論

2.4.1 參數(shù)設(shè)置對(duì)結(jié)果的影響

在DPPB算法中，搜索半徑、窗寬和步長(zhǎng)是算法的三個(gè)關(guān)鍵參數(shù)。搜索半徑的大小影響MMP算法對(duì)噪聲的魯棒性。減小可以降低索引區(qū)域內(nèi)的噪聲點(diǎn)數(shù)量，避免脊線波動(dòng)，但同時(shí)也降低了索引到正確脊線點(diǎn)的概率，容易在噪聲的干擾下偏離正確索引路徑；增大可以增加索引到正確脊線點(diǎn)的概率，當(dāng)索引路徑偏離理論脊線時(shí)有更大的機(jī)會(huì)重新返回正確的索引路徑，但是也增加了噪聲點(diǎn)的數(shù)量，容易使提取結(jié)果產(chǎn)生波動(dòng)。加窗求質(zhì)心的目的是降低噪聲對(duì)脊線提取的干擾，使質(zhì)心點(diǎn)能夠準(zhǔn)確描述理論脊線。窗口大小和步長(zhǎng)會(huì)影響質(zhì)心計(jì)算的準(zhǔn)確性、抗噪聲干擾能力和復(fù)雜度。隨著增大，質(zhì)心計(jì)算的抗噪聲干擾能力增強(qiáng)，而準(zhǔn)確性則先增高后降低；增大時(shí)，質(zhì)心計(jì)算的準(zhǔn)確性提高，但是提升效果會(huì)隨的減小而降低。同時(shí)，和影響算法的復(fù)雜度，增大和會(huì)增加運(yùn)算量。

目前尚未有合適的理論確定數(shù)據(jù)的最優(yōu)搜索半徑、窗寬和步長(zhǎng)，本文僅根據(jù)經(jīng)驗(yàn)和試驗(yàn)分析獲得各參數(shù)的較優(yōu)取值范圍。對(duì)于，一般取10～35 Hz效果較好。對(duì)于，如果待提取的脊線變化較為劇烈，通常選擇0.1～0.5 s的較小窗口；如果變化較為緩慢，通常選擇0.5～1 s的較大窗口。對(duì)于，原則上越小越好，但是考慮到運(yùn)算量，一般選擇。

2.4.2 算法的計(jì)算成本

本文所提算法的計(jì)算成本主要包括MMP索引、質(zhì)心求解和動(dòng)態(tài)路徑規(guī)劃三部分。若信號(hào)時(shí)頻分布的數(shù)據(jù)量為×，為離散信號(hào)總長(zhǎng)度，為頻率維度，算法的計(jì)算成本可以表示為：

（14）

式中 ，，和分別為算法總成本、MMP索引成本、質(zhì)心計(jì)算成本和動(dòng)態(tài)路徑規(guī)劃成本；為脊線段的數(shù)量；為個(gè)點(diǎn)的多項(xiàng)式擬合計(jì)算成本。

所提算法以較高的計(jì)算成本來換取較好的提取性能，所以計(jì)算成本較大。特別注意的是，的系數(shù)對(duì)計(jì)算成本產(chǎn)生較大的影響，所以在保證索引精度的前提下，合理減小和、增大可以有效降低算法的計(jì)算成本。

3 仿真分析

3.1 構(gòu)建仿真模型

設(shè)仿真信號(hào)x及瞬時(shí)頻率曲線f如下式所示：

（15）

設(shè)置采樣頻率fs=800 Hz，采樣時(shí)間t=10 s，信號(hào)的時(shí)域圖和STFT變換后的時(shí)頻分布如圖6所示。

對(duì)信號(hào)加高斯白噪聲后，取，運(yùn)用第1節(jié)MMP脊線提取算法得到脊線集如圖7所示。

由式（7）構(gòu)建質(zhì)心稀疏矩陣V，其中總采樣點(diǎn)數(shù)N=fs·t=8000，窗寬度g=0.5 s，步長(zhǎng)λ=0.1 s，時(shí)間序列總長(zhǎng)度B=t/λ=100。

令式（10）中的，，，，，則仿真模型的動(dòng)態(tài)路徑規(guī)劃函數(shù)可表示為：

（16）

根據(jù)式（11）構(gòu)造仿真信號(hào)的代價(jià)函數(shù)為：

（17）

式中 為歸一化質(zhì)心質(zhì)量；為距離權(quán)重，用來平衡歸一化質(zhì)心質(zhì)量和距離對(duì)代價(jià)函數(shù)的貢獻(xiàn)度，通常的取值范圍為。

3.2 時(shí)頻脊線提取效果衡量指標(biāo)

為了更直觀地分析不同算法的優(yōu)劣，此處構(gòu)造兩個(gè)量化衡量指標(biāo)：相似度系數(shù)和置信度。

絕對(duì)值倒數(shù)法是確定兩個(gè)向量間相似度的一種方法，對(duì)于提取的第j條脊線上任意一點(diǎn)的頻率與對(duì)應(yīng)的理論值之間的相關(guān)程度可以由下式確定：

（18）

式中 c為誤差允許值，表示接受與理論值的距離不大于c的提取結(jié)果。

對(duì)于某一算法O提取的第j條脊線與理論曲線的相似度，計(jì)算公式如下：

（19）

式中 Ra∈［0，1］，Ra越接近1，說明算法提取的脊線與理論曲線相似度越高，脊線提取效果越好；反之Ra越接近0，則說明相似度越低，提取效果越差。

圖8所示為不同相似度Ra時(shí)的時(shí)頻脊線對(duì)比圖，圖中誤差上、下限之間為允許誤差c的范圍。當(dāng)Ra≥0.9時(shí)，脊線與理論時(shí)頻曲線基本吻合，提取結(jié)果誤差較??；當(dāng)0.9≤Ra≤0.8時(shí)，脊線在理論曲線附近波動(dòng)變大，且存在少量脊線點(diǎn)超出了誤差范圍，但是脊線基本落在誤差范圍以內(nèi)；當(dāng)Ra≤0.7時(shí)，脊線與理論曲線保持相同的變化趨勢(shì)，但是有部分脊線嚴(yán)重超出誤差范圍，提取效果欠佳。

運(yùn)用算法O對(duì)M個(gè)信號(hào)進(jìn)行脊線提取，獲得M條脊線，根據(jù)式（19）計(jì)算得到所有脊線與理論曲線的相關(guān)系數(shù)，算法O的平均相似度系數(shù)為：

（20）

對(duì)于任意，存在k條脊線滿足，則置信水平在上的置信度為：

（21）

平均相似系數(shù)表示提取結(jié)果與理論結(jié)果偏差的統(tǒng)計(jì)平均值，可以反映算法的穩(wěn)定程度。置信度表示將與理論結(jié)果的相似度在的提取結(jié)果視為合格結(jié)果，則在整個(gè)試驗(yàn)結(jié)果中合格的提取結(jié)果所占的比率為，值越大，代表算法的提取效果越好。

3.3 仿真結(jié)果及分析

脊線提取結(jié)果及其與理論曲線對(duì)比結(jié)果如圖9所示。圖9（a）中紅色曲線為提取到的時(shí)頻脊線，脊線索引到了理論脊線附近所有的高幅值時(shí)頻點(diǎn)；從圖9（b）中可以看到，提取到的脊線與理論曲線基本重合，表明本算法可以提取得到較為準(zhǔn)確的時(shí)頻脊線。

為了分析噪聲對(duì)脊線提取結(jié)果的影響，運(yùn)用峰值索引、MMP脊線索引、DPPB脊線索引三種算法對(duì)信噪比的信號(hào)進(jìn)行脊線提取，結(jié)果如圖10所示。

在圖10（a）中，對(duì)于信噪比SNR=-5 dB的低噪聲信號(hào)，3種算法均能較好地提取時(shí)頻脊線。在SNR=-8和-10 dB的較強(qiáng)噪聲時(shí)（圖10（b）和（c）所示），峰值索引算法得到的脊線已經(jīng)完全偏離理論曲線，MMP算法得到的脊線與理論曲線有相同的變化趨勢(shì)，但存在局部誤差而偏離理論曲線，DPPB算法提取的脊線與理論曲線基本重合。在SNR=-12 dB的強(qiáng)噪聲下（圖10（d）所示），其他兩種算法已經(jīng)完全偏離理論曲線，DPPB算法可得到與理論曲線變化趨勢(shì)相同的脊線，僅存在局部誤差。由此可以初步得出，峰值索引算法和MMP算法容易受噪聲的影響，強(qiáng)噪聲下會(huì)失去脊線提取能力，而DPPB算法可以在強(qiáng)噪聲下有效提取時(shí)頻脊線，提取效果和對(duì)噪聲的魯棒性優(yōu)于其他兩種算法。

重復(fù)進(jìn)行圖10中的仿真試驗(yàn)，試驗(yàn)次數(shù)M=10000，根據(jù)式（20）計(jì)算三種算法的平均相關(guān)系數(shù)如表1所示，其中取 Hz。

由表1可知，當(dāng)SNR=-5 dB時(shí)，三種算法的平均相關(guān)系數(shù)均大于0.9，提取結(jié)果比較穩(wěn)定。峰值索引算法在SNR=-8，-10和-12 dB的強(qiáng)噪聲下值下降明顯；MMP算法在三種強(qiáng)噪聲下的值雖然高于峰值索引算法，但是也有較大幅度的下降；DPPB算法在三種強(qiáng)噪聲下的值降幅較小且均在0.8以上。隨著噪聲的增強(qiáng)，三種算法的值均有不同程度的下降，但是DPPB算法的值始終高于其他兩種算法，且噪聲越強(qiáng)越明顯，說明對(duì)噪聲的魯棒性優(yōu)于其他算法。

根據(jù)圖8的仿真結(jié)果，取，由式（21）計(jì)算置信水平在［0.9，1］上的置信度，結(jié)果如表2所示。

由表2可知，在SNR=-5 dB時(shí)，MMP脊線索引和DPPB算法的置信度分別為98.01%和100%，可以準(zhǔn)確地提取到符合要求的時(shí)頻脊線；峰值索引算法的也達(dá)到80%，具備較好的提取效果。當(dāng)噪聲強(qiáng)度高于-8 dB時(shí)，峰值索引算法的均為0，已經(jīng)無法提取到符合要求的脊線。MMP算法在SNR=-8 dB時(shí)為54.75%，接近一半的提取結(jié)果不符合要求，在SNR=-10和-12 dB時(shí)的分別為13.95%和2.87%，已經(jīng)很難提取到合格的脊線。DPPB算法在SNR=-8和-10 dB的較強(qiáng)噪聲條件下分別為99.91%和93.51%，合格脊線所占的比率均在90%以上，在SNR=-12 dB的強(qiáng)噪聲下為57.12%，依然有一半以上的提取結(jié)果符合要求。雖然噪聲的增強(qiáng)降低了各算法的提取效果，但DPPB算法的提取效果始終優(yōu)于其他兩種算法，且針對(duì)強(qiáng)噪聲具有良好的脊線提取能力。

通過以上仿真結(jié)果和分析可以得出，DPPB算法可以抵抗信號(hào)中強(qiáng)噪聲的干擾，相比于峰值索引算法和MMP算法更容易從強(qiáng)噪聲信號(hào)中提取到有效的時(shí)頻脊線，對(duì)噪聲有更好的魯棒性。

4 試驗(yàn)驗(yàn)證

4.1 試驗(yàn)說明

齒輪箱振動(dòng)信號(hào)中的嚙合頻率分量是分析故障的重要特征頻率之一，本文利用SpectraQuest公司研發(fā)的動(dòng)力傳動(dòng)故障診斷綜合試驗(yàn)臺(tái)（DDS）獲取不同轉(zhuǎn)速下齒輪箱的振動(dòng)信號(hào)，通過提取不同轉(zhuǎn)速下的嚙合頻率時(shí)頻脊線驗(yàn)證所提算法的有效性，試驗(yàn)臺(tái)如圖11所示。

其中齒輪箱5為二級(jí)減速器減速齒輪箱，由4組直齒輪組合而成，主要參數(shù)如表3所示。設(shè)置加速度振動(dòng)傳感器6的采樣頻率為25.6 kHz，采集齒輪箱徑向的振動(dòng)信號(hào)。

由于DDS試驗(yàn)臺(tái)不易輸出非線性變化的轉(zhuǎn)速，所以試驗(yàn)?zāi)M三組線性變化的輸入轉(zhuǎn)速，分別為震蕩波動(dòng)曲線、三角形單波谷波動(dòng)曲線、三角形三波谷波動(dòng)曲線，如圖12所示。為了與時(shí)頻脊線對(duì)應(yīng)，下文均用“轉(zhuǎn)頻”代替“轉(zhuǎn)速”。

齒輪箱振動(dòng)信號(hào)中包含與轉(zhuǎn)頻相關(guān)的第i級(jí)傳動(dòng)齒輪嚙合頻率及其倍頻，與的關(guān)系如下式所示：

（22）

因此，獲得齒輪箱振動(dòng)信號(hào)的第i級(jí)傳統(tǒng)齒輪嚙合頻率的n倍頻，便可計(jì)算得到齒輪箱的輸入轉(zhuǎn)頻。

4.2 試驗(yàn)結(jié)果及分析

為了計(jì)算方便，試驗(yàn)通過提取1級(jí)傳動(dòng)齒輪嚙合頻率的基頻計(jì)算輸入轉(zhuǎn)頻。三種不同輸入轉(zhuǎn)速下齒輪箱的振動(dòng)時(shí)域圖和STFT變換后附近的時(shí)頻分布如圖13所示。

運(yùn)用DPPB算法提取不同轉(zhuǎn)頻曲線下嚙合頻率的時(shí)頻脊線，提取結(jié)果如圖14所示。

在圖13（b）和（d）的時(shí)頻分布中，信號(hào)的SNR較高，噪聲對(duì)瞬時(shí)頻率的影響較小，可以觀察到較為明顯的時(shí)頻脊線；在圖13（f）中，頻率波動(dòng)部分受噪聲污染嚴(yán)重，SNR較低，難以觀察到時(shí)頻脊線的變化趨勢(shì)。通過圖14的提取結(jié)果可以看出，雖然三種轉(zhuǎn)頻下的信噪比差異較大，但是DPPB算法均能提取到較為準(zhǔn)確的時(shí)頻脊線，表明該算法對(duì)強(qiáng)噪聲有良好的魯棒性。

為了驗(yàn)證DPPB算法提取轉(zhuǎn)頻的準(zhǔn)確性，將圖14得到的嚙合頻率和表3中齒輪箱參數(shù)及代入式（22）可計(jì)算三組信號(hào)的轉(zhuǎn)頻。同時(shí)，運(yùn)用MMP算法提取的時(shí)頻脊線，計(jì)算轉(zhuǎn)頻；用試驗(yàn)臺(tái)中的角度編碼器測(cè)得齒輪箱輸出軸轉(zhuǎn)頻，根據(jù)傳動(dòng)比計(jì)算得到輸入軸實(shí)測(cè)轉(zhuǎn)頻。將，作為對(duì)照組與進(jìn)行對(duì)比，結(jié)果如圖15所示。

由圖15可以看出，在震蕩波動(dòng)轉(zhuǎn)頻和三角形單波谷轉(zhuǎn)頻下，MMP脊線索引和DPPB兩種算法提取到的時(shí)頻脊線基本與實(shí)測(cè)轉(zhuǎn)頻重合（圖15（a）和（b）所示），但由MMP脊線索引算法提取的三角形單波谷波時(shí)頻脊線存在較小的誤差。在三波谷轉(zhuǎn)頻波動(dòng)下，MMP脊線索引算法提取的脊線在恒頻部分與重合較好，但波谷處與存在較大誤差；DPPB算法得到的脊線無論在恒頻部分還是波谷處均與重合較好，只有局部存在較小的誤差（如圖15（c）所示）。由此表明DPPB算法可以在強(qiáng)噪聲下較為準(zhǔn)確地估計(jì)信號(hào)的轉(zhuǎn)頻曲線。

將實(shí)測(cè)轉(zhuǎn)頻作為理論輸入轉(zhuǎn)頻，根據(jù)式（19）計(jì)算三組信號(hào)在不同算法下提取結(jié)果的相似度Ra衡量MMP算法和DPPB算法的提取效果，如表4所示。

由表4可以看出，MMP算法和DPPB算法得到的震蕩波動(dòng)信號(hào)轉(zhuǎn)頻的Ra值均為1，三角形單波谷波動(dòng)信號(hào)轉(zhuǎn)頻的Ra值分別為0.9939和1，均有較好的提取效果。對(duì)于三角形三波谷波動(dòng)信號(hào)，MMP算法得到的轉(zhuǎn)頻的Ra值為0.7678，與實(shí)際轉(zhuǎn)頻的偏差較大，而DPPB算法得到的轉(zhuǎn)頻的Ra值為0.9669，遠(yuǎn)高于MMP算法，進(jìn)一步說明DPPB算法的有效性和對(duì)噪聲的魯棒性。

5 結(jié) 論

mhqytu3B96llJVmdX9+jYw==本文提出了一種改進(jìn)的時(shí)頻脊線索引算法來提高時(shí)頻估計(jì)的準(zhǔn)確性，主要結(jié)論如下：

（1） 信號(hào)時(shí)頻分布在強(qiáng)噪聲下會(huì)產(chǎn)生多個(gè)幅值極值，MMP算法索引時(shí)頻脊線時(shí)容易受到噪聲極值點(diǎn)的干擾陷入局部最優(yōu)，導(dǎo)致算法提取到錯(cuò)誤的時(shí)頻脊線。

（2） 通過構(gòu)建脊線質(zhì)心稀疏矩陣，將脊線索引對(duì)象由時(shí)頻離散點(diǎn)轉(zhuǎn)化為質(zhì)心點(diǎn)，降低了噪聲對(duì)脊線索引的干擾，同時(shí)壓縮了索引數(shù)據(jù)量，提高了運(yùn)算速率。

（3） 提出了一種基于質(zhì)心動(dòng)態(tài)路徑規(guī)劃（DPPB）的時(shí)頻脊線索引算法，設(shè)計(jì)動(dòng)態(tài)路徑優(yōu)化函數(shù)和脊線索引代價(jià)函數(shù)識(shí)別質(zhì)心稀疏矩陣中的時(shí)頻脊線點(diǎn)，從而避免了脊線索引過程中陷入局部最優(yōu)，提高了算法對(duì)噪聲的魯棒性。通過仿真和試驗(yàn)分析，驗(yàn)證了DPPB算法索引時(shí)頻脊線的精度和魯棒性優(yōu)于傳統(tǒng)的脊線索引算法。所提算法可用于實(shí)際工程中變轉(zhuǎn)速齒輪箱的轉(zhuǎn)速估計(jì)和故障診斷，具有一定的工程應(yīng)用價(jià)值。

參考文獻(xiàn)：

［1］ 陳雪峰， 郭艷婕， 許才彬， 等. 風(fēng)電裝備故障診斷與健康監(jiān)測(cè)研究綜述［J］. 中國(guó)機(jī)械工程， 2020， 31（2）： 175-189.

Cheng Xuefeng， Guo Yanjie， Xu Caibin， et al. Review of fault diagnosis and health monitoring for wind power equipment［J］. China Mechanical Engineering， 2020， 31（2）： 175-189.

［2］ 孫鑫暉， 董翔文， 郝木明. 基于小波脊提取的扭轉(zhuǎn)振動(dòng)測(cè)試方法研究［J］. 振動(dòng)、測(cè)試與診斷， 2020， 40（1）： 135-139.

Sun Xinhui， Dong Xiangwen， Hao Muming. Research on testing method of torsional vibration based on wavelet ridge extraction［J］. Journal of Vibration，Measurement & Diagnosis， 2020， 40（1）： 135-139.

［3］ 趙德尊， 王天楊， 褚福磊. 基于自適應(yīng)廣義解調(diào)變換的滾動(dòng)軸承時(shí)變非平穩(wěn)故障特征提?。跩］. 機(jī)械工程學(xué)報(bào)， 2020， 56（3）： 80-87.

Zhao Dezun， Wang Tianyang， Chu Fulei. Adaptive generalized demodulation transform based rolling bearing time-varying nonstationary fault feature extraction［J］. Journal of Mechanical Engineering， 2020， 56（3）： 80-87.

［4］ 王曉龍， 閆曉麗， 何玉靈. 變速工況下基于IEWT能量階次譜的風(fēng)電機(jī)組軸承故障診斷［J］. 太陽能學(xué)報(bào)， 2021， 42（4）： 479-486.

Wang Xiaolong， Yan Xiaoli， He Yuling. Fault diagnosis wind turbine bearing based on IEWT energy order spec-trum under variable speed condition［J］. Acta Energiae Solaris Sinica， 2021， 42（4）： 479-486.

［5］ 陳是扦， 彭志科， 周鵬. 信號(hào)分解及其在機(jī)械故障診斷中的應(yīng)用研究綜述［J］. 機(jī)械工程學(xué)報(bào)， 2020， 56（17）： 91-107.

Chen Shiqian， Peng Zhike， Zhou Peng. Review of signal decomposition theory and its applications in machine fault diagnosis［J］. Journal of Mechanicla Egnineering， 2020， 56（17）： 91-107.

［6］ 郭瑜， 秦樹人， 湯寶平， 等. 基于瞬時(shí)頻率估計(jì)的旋轉(zhuǎn)機(jī)械階比跟蹤［J］. 機(jī)械工程學(xué)報(bào)， 2003， 39（3）： 32-36.

Guo Yu， Qin Shuren， Tang Baoping， et al. Order tracking of rotating machinery based on instantaneous frequency estimation［J］. Journal of Mechanical Engineering， 2003， 39（3）： 32-36.

［7］ Iatsenko D， McClintock P V E， Stefanovska A. Extrac-tion of instantaneous frequencies from ridges in time–frequency representations of signals［J］. Signal Processing， 2016， 125： 290-303.

［8］ 江星星， 李舜酩， 周東旺， 等. 時(shí)頻脊融合方法及時(shí)變工況行星齒輪箱故障識(shí)別［J］. 振動(dòng)工程學(xué)報(bào)， 2017， 30（1）： 127-134.

Jiang Xingxing， Li Shunming， Zhou Dongwang， et al. Time-frequency ridge fusion method and defective identification of planetary gearbox running on time-varying condition［J］. Journal of Vibration Engineering， 2017， 30（1）： 127-134.

［9］ Ding Rongmei， Shi Juanjuan， Jiang Xingxing， et al. Multiple instantaneous frequency ridge based integration strategy for bearing fault diagnosis under variable speed operations［J］. Measurement Science and Technology， 2018， 29（11）： 115002.

［10］ 張焱， 何姝鋇， 王平， 等. 無轉(zhuǎn)速計(jì)下變工況滾動(dòng)軸承故障特征量化表征提?。跩］. 儀器儀表學(xué)報(bào)， 2021， 42（8）： 104-114.

Zhang Yan， He Shubei， Wang Ping， et al. Tacholess quanti-tative characterization of rolling bearing fault feature under varying conditions［J］. Chinese Journal of Scientific Instrument， 2021， 42（8）： 104-114.

［11］ Wu M， Wu F， Yang K， et al. A multipath matching pursuit algorithm based on improved-inner product matching criterion［C］.2020 IEEE International Conference on Signal Processing， Communications and Computing （ICSPCC）. IEEE， 2020： 1-5.

［12］ Fourer D， Auger F， Flandrin P. Recursive versions of the Levenberg-Marquardt reassigned spectrogram and of the synchrosqueezed STFT［C］.2016 IEEE International Conference on Acoustics， Speech and Signal Processing （ICASSP）. IEEE， 2016： 4880-4884.

［13］ Yu G， Wang Z H， Zhao P. Multisynchrosqueezing transform［J］. IEEE Transactions on Industrial Electronics， 2018， 66（7）： 5441-5455.

［14］ Yu G， Yu M J， Xu C Y. Synchroextracting transform［J］. IEEE Transactions on Industrial Electronics， 2017， 64（10）： 8042-8054.

［15］ Li M F， Wang T Y， Kong Y， et al. Synchro-reassigning transform for instantaneous frequency estimation and signal reconstruction［J］. IEEE Transactions on Industrial Electronics， 2021， 69（7）： 7263-7274.

［16］ 江星星. 齒輪箱關(guān)鍵部件非平穩(wěn)振動(dòng)信號(hào)分析及診斷方法研究［D］. 南京：南京航空航天大學(xué)， 2016.

Jiang Xingxing. Research on the nonstationary vibration signal of key parts of gearbox and its fault diagnosis method［D］. Nanjing：Nanjing University of Aeronautics and Astronautics， 2016.

［17］ 楊世偉， 蔣國(guó)平， 宋玉蓉， 等. 基于GPU的稀疏矩陣存儲(chǔ)格式優(yōu)化研究［J］. 計(jì)算機(jī)工程， 2019， 45（9）： 23-31.

Yang Shiwei， Jiang Guoping， Song Yurong， et al. Research on storage format optimization of sparse matrix based on GPU［J］. Computer Engineering， 2019， 45（9）： 23-31.

Time-frequency ridge index algorithm for gearbox under variable speed based on dynamic path planning of barycenter

ZHANG Bo-lin1， WAN Shu-ting1， ZHAO Xiao-yan1， ZHANG Xiong1， GU Xiao-hui2

（1.Hebei Key Laboratory of Electric Machinery Health Maintenance & Failure Prevention，North China Electric Power University， Baoding 071003， China；2.State Key Laboratory of Mechanical Behavior and System Safety of Traffic Engineering Structures，Shijiazhuang Tiedao University， Shijiazhuang 050043， China）

Abstract： The paper proposes a time-frequency ridge index algorithm for gearboxes under variable speed conditions， based on Dynamic Path Planning of Barycenter （DPPB）. This algorithm addresses the challenge of estimating the instantaneous frequency of signals in a high-noise environment. The algorithm builds upon the analysis of the Multi-Path Matching Pursuit （MMP） ridge index algorithm and its limitations under high noise. By adding windows to the ridge set obtained by the MMP algorithm， a ridge barycenter sparse matrix of the signal is constructed. A dynamic path planning function is then designed for the barycenter sparse matrix to index the barycenters on the ridge line. The optimal time-frequency ridge line is calculated based on the values of the ridge line cost function. The similarity coefficient Ra and confidence σRa are used as measures of the ridge extraction effect. Simulations and experiments indicate that the DPPB algorithm can effectively extract the time-frequency ridge of signals in high-noise environments， and it is more reliable and robust than the peak index algorithm and the MMP algorithm under various noise intensities.

Key words： fault diagnosis； gearbox；dynamic path planning barycenter （DPPB）； time-frequency ridge index； multipath matching pursuit （MMP）； time-frequency analysis

作者簡(jiǎn)介： 張伯麟（1993—），男，博士研究生。電話：（0312）7525428； E-mail：393696838@qq.com。

通訊作者： 萬書亭（1970—），男，博士，教授，博士生導(dǎo)師。電話：（0312）7525455； E-mail：52450809@ncepu.edu.cn。