本發(fā)明涉及旋轉(zhuǎn)機(jī)械狀態(tài)監(jiān)測與故障診斷領(lǐng)域，具體涉及基于log算子和pauta檢驗(yàn)的瞬時頻率估計(jì)方法。

背景技術(shù)：

由于工作環(huán)境的復(fù)雜性，旋轉(zhuǎn)機(jī)械經(jīng)常工作在變速條件下。瞬時頻率估計(jì)是評估旋轉(zhuǎn)機(jī)械運(yùn)行狀態(tài)及進(jìn)行故障診斷的重要前提。目前常用的瞬時頻率估計(jì)方法是一步代價函數(shù)法（one-stepcostfunction）。一步代價函數(shù)法能夠在局部頻率范圍內(nèi)搜索脊點(diǎn)，但是局部頻率范圍的中心點(diǎn)依賴于上一個脊點(diǎn)的位置，這導(dǎo)致一步代價函數(shù)存在著延遲。此外，局部頻率范圍的寬度是根據(jù)經(jīng)驗(yàn)隨意設(shè)置的，且在任何時刻的寬度都是固定的，不能隨著時間發(fā)生變化，這導(dǎo)致一步代價函數(shù)缺乏足夠的自適應(yīng)性。上述缺陷導(dǎo)致一步代價函數(shù)法在估計(jì)瞬時頻率時準(zhǔn)確度和精確度較低。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明要解決的問題是針對以上不足，提出基于log算子和pauta檢驗(yàn)的瞬時頻率估計(jì)方法。與現(xiàn)有方法相比，本發(fā)明將映射后的平滑均值曲線作為目標(biāo)脊線的參考線，將映射后的置信區(qū)間作為目標(biāo)脊線的局部搜索區(qū)間，因此局部頻率搜索范圍的中心點(diǎn)不依賴上一個脊點(diǎn)的位置，沒有任何延遲，局部頻率搜索范圍能夠自動設(shè)定，搜索帶寬能夠隨著時間的變化而自動變化，瞬時頻率估計(jì)結(jié)果的準(zhǔn)確度和精確度高。

為解決以上技術(shù)問題，本發(fā)明提供基于log算子和pauta檢驗(yàn)的瞬時頻率估計(jì)方法，其特征在于，包括以下步驟：

步驟1：采用短時傅里葉變換算法將信號x（k）（k=1,2,…,n）轉(zhuǎn)換為時頻譜圖，n代表信號的長度；

步驟2：從時頻譜圖中選取一塊具有較高信噪比的局部區(qū)域，采用log檢測算法將該局部區(qū)域轉(zhuǎn)換成二值圖像，二值圖像包含多條脊帶；局部區(qū)域是指至少包含兩條脊帶，信噪比大于80db的區(qū)域；

步驟3：采用pauta檢驗(yàn)算法排除每條脊帶上下邊緣的異常值；

步驟4：將上述多條脊帶按照相互之間的運(yùn)動學(xué)比例關(guān)系疊加到其中一條輪廓最完整的脊帶上，構(gòu)建一條具有完整清晰邊緣的合成脊帶；運(yùn)動學(xué)比例關(guān)系是指脊帶所對應(yīng)的機(jī)器部件之間的傳動比；

步驟5：采用pauta檢驗(yàn)算法排除上述合成脊帶上下邊緣的異常值；

步驟6：計(jì)算上述合成脊帶的均值曲線，采用五點(diǎn)三次平滑算法對均值曲線進(jìn)行平滑處理，得到平滑均值曲線，計(jì)算該平滑均值曲線在95%置信水平上的置信區(qū)間；

步驟7：將上述平滑均值曲線及其置信區(qū)間按照平滑均值曲線與待估計(jì)目標(biāo)脊線之間的運(yùn)動學(xué)比例關(guān)系映射到目標(biāo)脊線上；

步驟8：將映射后的平滑均值曲線作為目標(biāo)脊線的參考線，將映射后的置信區(qū)間作為目標(biāo)脊線的局部搜索區(qū)間；

步驟9：采用非延遲代價函數(shù)在每個時刻所對應(yīng)的局部搜索區(qū)間內(nèi)搜索脊點(diǎn)，確定每個時刻所對應(yīng)的瞬時頻率，最后得到整個時間區(qū)間上的瞬時頻率。

進(jìn)一步地，所述步驟1中短時傅里葉變換算法包括以下步驟：

1）對信號x（k）進(jìn)行短時傅里葉變換：

，

tf(t,f)代表信號x（k）的短時傅里葉變換結(jié)果，t代表時間因子，f代表尺度因子，函數(shù)w(z)代表自變量為z的窗口函數(shù)；

2）計(jì)算信號x（k）的時頻譜：

，

spectrogram(t,f)代表x（k）的時頻譜。

進(jìn)一步地，所述步驟2中l(wèi)og檢測算法包括以下步驟：

1)采用高斯濾波器對原始圖像i(x,y)進(jìn)行平滑濾波：

，

參數(shù)gσ(x,y)代表標(biāo)準(zhǔn)差為σ的高斯核函數(shù)，l(x,y)代表經(jīng)高斯濾波后的圖像，x代表圖像的時間點(diǎn)，y代表圖像的頻率點(diǎn)，*代表卷積計(jì)算；本發(fā)明中，σ=1；

2)對l(x,y)執(zhí)行l(wèi)aplacian運(yùn)算：

,

代表l(x,y)經(jīng)laplacian運(yùn)算后得到的結(jié)果；

3)設(shè)定一個合適的閾值，如果圖像上某一個點(diǎn)(x,y)所對應(yīng)的大于該閾值，則判定該點(diǎn)為邊緣；本發(fā)明中，閾值設(shè)置為4.63×10^-7。

進(jìn)一步地，所述步驟3中pauta檢驗(yàn)算法包括以下步驟：

1）估計(jì)信號xn（n=1,2,…,n）的標(biāo)準(zhǔn)差，

，

代表樣本均值，σ代表樣本標(biāo)準(zhǔn)差，n代表樣本長度；

2）如果，則剔除xn。

進(jìn)一步地，所述步驟9中非延遲代價函數(shù)包括以下步驟：

1）第k個時刻所對應(yīng)的局部搜索區(qū)間fbk定義為

，

fk(pmc)代表映射后的平滑均值曲線在第k個時刻的值，代表映射后的平滑均值曲線置信區(qū)間在第k個時刻寬度的一半，m代表目標(biāo)脊線的長度；

2）第k個時刻所對應(yīng)的非延遲代價函數(shù)cfk定義為：

,

，

fk(i)代表在fbk范圍內(nèi)所取的頻率值，tf(tk,fk)代表tf(t,f)在第k個時刻的值，tk代表t在第k個時刻的值，fk代表f在第k個時刻的值，ek代表權(quán)重因子。

進(jìn)一步地，相對誤差≤0.775%，平均相對誤差≤0.068%。

本發(fā)明采用以上技術(shù)方案，與現(xiàn)有技術(shù)相比，本發(fā)明具有以下優(yōu)點(diǎn)：

1)本發(fā)明具有實(shí)時性：本發(fā)明以合成脊帶平滑均值曲線映射為參考線，能夠即時確定當(dāng)前時刻局部頻率搜索范圍的中心點(diǎn)，避免了對前一個脊點(diǎn)的依賴，消除了時間延遲，具有實(shí)時性。

2)本發(fā)明具有自適應(yīng)性：本發(fā)明利用合成脊帶平滑均值曲線置信區(qū)間映射所提供的局部范圍，能夠自適應(yīng)確定每個時刻所對應(yīng)的局部頻率搜索范圍，搜索帶寬能夠隨著時間的變化而自動變化，不需要憑借經(jīng)驗(yàn)設(shè)置搜索帶寬，從而消除了由于人為原因而產(chǎn)生的誤差。

3)實(shí)驗(yàn)結(jié)果表明：由本發(fā)明得到的瞬時頻率估計(jì)值與實(shí)測值之間的最大相對誤差為0.775%，平均相對誤差為0.068%；與一步代價函數(shù)法的結(jié)果相比，最大相對誤差降低95.27%，平均相對誤差降低96.82%。

下面結(jié)合附圖和實(shí)施例對本發(fā)明做進(jìn)一步說明。

附圖說明

附圖1為本發(fā)明實(shí)施例中基于log算子和pauta檢驗(yàn)的瞬時頻率估計(jì)方法的流程圖；

附圖2為本發(fā)明實(shí)施例中行星齒輪箱振動信號；

附圖3為本發(fā)明實(shí)施例中行星齒輪箱振動信號的時頻譜；

附圖4為本發(fā)明實(shí)施例中從時頻譜圖上選取的具有高信噪比的局部區(qū)域；

附圖5為本發(fā)明實(shí)施例中由log算法檢測到的局部圖像區(qū)域的邊緣；

附圖6為本發(fā)明實(shí)施例中采用pauta檢驗(yàn)算法消除每條脊帶異常點(diǎn)后的結(jié)果；

附圖7為本發(fā)明實(shí)施例中利用脊帶之間的運(yùn)動學(xué)比例關(guān)系疊加而成的合成脊帶（最下層的脊帶即為合成脊帶）；

附圖8為本發(fā)明實(shí)施例中采用pauta檢驗(yàn)算法消除合成脊帶異常點(diǎn)后的結(jié)果；

附圖9為本發(fā)明實(shí)施例中合成脊帶的均值平滑曲線及其95%置信區(qū)間；

附圖10為本發(fā)明實(shí)施例中映射得到的均值平滑曲線及其置信區(qū)間；

附圖11為本發(fā)明實(shí)施例中瞬時頻率估計(jì)值。

具體實(shí)施方式

實(shí)施例，如圖1所示，基于log算子和pauta檢驗(yàn)的瞬時頻率估計(jì)方法，包括以下步驟：

步驟1：采用短時傅里葉變換算法將信號x（k）（k=1,2,…,n）轉(zhuǎn)換為時頻譜圖，n代表信號的長度；

步驟2：從時頻譜圖中選取一塊具有較高信噪比的局部區(qū)域，采用log檢測算法將該局部區(qū)域轉(zhuǎn)換成二值圖像，二值圖像包含多條脊帶；局部區(qū)域是指至少包含兩條脊帶，信噪比大于80db的區(qū)域；

步驟3：采用pauta檢驗(yàn)算法排除每條脊帶上下邊緣的異常值；

步驟4：將上述多條脊帶按照相互之間的運(yùn)動學(xué)比例關(guān)系疊加到其中一條輪廓最完整的脊帶上，構(gòu)建一條具有完整清晰邊緣的合成脊帶；運(yùn)動學(xué)比例關(guān)系是指脊帶所對應(yīng)的機(jī)器部件之間的傳動比；

步驟5：采用pauta檢驗(yàn)算法排除上述合成脊帶上下邊緣的異常值；

步驟6：計(jì)算上述合成脊帶的均值曲線，采用五點(diǎn)三次平滑算法對均值曲線進(jìn)行平滑處理，得到平滑均值曲線，計(jì)算該平滑均值曲線在95%置信水平上的置信區(qū)間；

步驟7：將上述平滑均值曲線及其置信區(qū)間按照平滑均值曲線與待估計(jì)目標(biāo)脊線之間的運(yùn)動學(xué)比例關(guān)系映射到目標(biāo)脊線上；

步驟8：將映射后的平滑均值曲線作為目標(biāo)脊線的參考線，將映射后的置信區(qū)間作為目標(biāo)脊線的局部搜索區(qū)間；

步驟9：采用非延遲代價函數(shù)在每個時刻所對應(yīng)的局部搜索區(qū)間內(nèi)搜索脊點(diǎn)，確定每個時刻所對應(yīng)的瞬時頻率，最后得到整個時間區(qū)間上的瞬時頻率。

步驟1中短時傅里葉變換算法包括以下步驟：

1）對信號x（k）進(jìn)行短時傅里葉變換：

，

tf(t,f)代表信號x（k）的短時傅里葉變換結(jié)果，t代表時間因子，f代表尺度因子，函數(shù)w(z)代表自變量為z的窗口函數(shù)；

2）計(jì)算信號x（k）的時頻譜：

，

spectrogram(t,f)代表x（k）的時頻譜。

步驟2中l(wèi)og檢測算法包括以下步驟：

1)采用高斯濾波器對原始圖像i(x,y)進(jìn)行平滑濾波：

，

參數(shù)gσ(x,y)代表標(biāo)準(zhǔn)差為σ的高斯核函數(shù)，l(x,y)代表經(jīng)高斯濾波后的圖像，x代表圖像的時間點(diǎn)，y代表圖像的頻率點(diǎn)，*代表卷積計(jì)算；本發(fā)明中，σ=1；

2)對l(x,y)執(zhí)行l(wèi)aplacian運(yùn)算：

,

代表l(x,y)經(jīng)laplacian運(yùn)算后得到的結(jié)果；

3)設(shè)定一個合適的閾值，如果圖像上某一個點(diǎn)(x,y)所對應(yīng)的大于該閾值，則判定該點(diǎn)為邊緣；本發(fā)明中，閾值設(shè)置為4.63×10^-7。

步驟3中pauta檢驗(yàn)算法包括以下步驟：

1）估計(jì)信號xn（n=1,2,…,n）的標(biāo)準(zhǔn)差，

，

代表樣本均值，σ代表樣本標(biāo)準(zhǔn)差，n代表樣本長度；

2）如果，則剔除xn。

步驟9中非延遲代價函數(shù)包括以下步驟：

1）第k個時刻所對應(yīng)的局部搜索區(qū)間fbk定義為

，

fk(pmc)代表映射后的平滑均值曲線在第k個時刻的值，代表映射后的平滑均值曲線置信區(qū)間在第k個時刻寬度的一半，m代表目標(biāo)脊線的長度；

2）第k個時刻所對應(yīng)的非延遲代價函數(shù)cfk定義為：

,

，

fk(i)代表在fbk范圍內(nèi)所取的頻率值，tf(tk,fk)代表tf(t,f)在第k個時刻的值，tk代表t在第k個時刻的值，fk代表f在第k個時刻的值，ek代表權(quán)重因子。

利用風(fēng)機(jī)渦輪行星齒輪箱振動數(shù)據(jù)對本發(fā)明所述算法的性能進(jìn)行驗(yàn)證。

振動數(shù)據(jù)從靠近行星齒輪系的齒輪箱外殼上采集，數(shù)據(jù)長度n=2736825，采樣頻率fs=5000hz。

采集到的行星齒輪箱振動數(shù)據(jù)如圖2所示。

采用短時傅里葉變換算法將圖2所示的行星齒輪箱振動數(shù)據(jù)轉(zhuǎn)換為時頻譜圖，得到的時頻譜圖如圖3所示。

從圖3所示的時頻譜圖上選取具有高信噪比的局部區(qū)域，得到的局部區(qū)域如圖4所示。

采用log檢測算法對如圖4所示的局部區(qū)域進(jìn)行邊緣檢測，得到的圖像邊緣如圖5所示。

采用pauta檢驗(yàn)算法消除圖5中各條脊帶的異常點(diǎn)，得到的結(jié)果如圖6所示。

按照脊帶之間的運(yùn)動學(xué)比例關(guān)系將各條脊帶疊加到其中一條輪廓最完整的脊帶上，所構(gòu)建的合成脊帶如圖7所示（最下層的脊帶即為合成脊帶）。

采用pauta檢驗(yàn)算法消除合成脊帶的異常點(diǎn)，結(jié)果如圖8所示。

計(jì)算合成脊帶的平滑均值曲線及其95%的置信區(qū)間，結(jié)果如圖9所示。

按照平滑均值曲線與目標(biāo)脊線之間的運(yùn)動學(xué)比例關(guān)系將平滑均值曲線及其置信區(qū)間映射到目標(biāo)脊線上，結(jié)果如圖10所示。

采用非延遲代價函數(shù)搜索目標(biāo)脊線的脊點(diǎn)，得到的瞬時頻率曲線如圖11所示。

經(jīng)多次實(shí)驗(yàn)表明，由本發(fā)明得到的瞬時頻率估計(jì)值與實(shí)測值之間的最大相對誤差為0.775%，平均相對誤差為0.068%，而采用一步代價函數(shù)法獲得的瞬時頻率估計(jì)值與實(shí)測值之間的最大相對誤差為16.39%，平均相對誤差為2.14%，本發(fā)明最大相對誤差降低95.27%，平均相對誤差降低96.82%。

根據(jù)實(shí)驗(yàn)結(jié)果，分析后認(rèn)為：

1)傳統(tǒng)的一步代價函數(shù)在確定當(dāng)前搜索區(qū)間的中心點(diǎn)時需要依賴上一個脊點(diǎn)的位置，存在著時間延遲現(xiàn)象，本發(fā)明利用映射后的平滑均值曲線作為參考線，可以即時確定當(dāng)前搜索區(qū)間的中心位置，完全不依賴上一個脊點(diǎn)，因此具有實(shí)時性。

2)傳統(tǒng)的的一步代價函數(shù)法缺乏自適應(yīng)性，需要人為設(shè)置搜索區(qū)間，且搜索寬度是固定的，因而不可避免地帶來誤差，本發(fā)明利用映射后的平滑均值曲線置信區(qū)間來自動確定局部搜索區(qū)間，搜索帶寬能夠隨著時間的變化而自動變化，不需要人工參與，因此具有自適應(yīng)性。

3)與傳統(tǒng)的一步代價函數(shù)法相比，本發(fā)明精確度和準(zhǔn)確度高。

本領(lǐng)域技術(shù)人員應(yīng)該認(rèn)識到，上述的具體實(shí)施方式只是示例性的，是為了使本領(lǐng)域技術(shù)人員能夠更好的理解本發(fā)明內(nèi)容，不應(yīng)理解為是對本發(fā)明保護(hù)范圍的限制，只要是根據(jù)本發(fā)明技術(shù)方案所作的改進(jìn)，均落入本發(fā)明的保護(hù)范圍。