潘崢嶸， 譙自健， 張 寧

（蘭州理工大學(xué)電氣與信息工程學(xué)院，甘肅 蘭州 730050）

基于符號(hào)序列熵的自適應(yīng)隨機(jī)共振的微弱信號(hào)檢測(cè)

潘崢嶸， 譙自健， 張 寧

（蘭州理工大學(xué)電氣與信息工程學(xué)院，甘肅 蘭州 730050）

針對(duì)現(xiàn)有隨機(jī)共振以信噪比為測(cè)度，難以進(jìn)行量化和在實(shí)際工程中應(yīng)用，提出一種符號(hào)序列改進(jìn)型香農(nóng)熵結(jié)合隨機(jī)共振的微弱信號(hào)檢測(cè)方法。介紹了隨機(jī)共振和符號(hào)序列化的基本原理，以輸出信號(hào)的符號(hào)序列熵作為判斷是否達(dá)到最佳隨機(jī)共振的測(cè)度，自適應(yīng)的調(diào)節(jié)系統(tǒng)參數(shù)a和b，使系統(tǒng)達(dá)到最佳共振并進(jìn)行頻譜分析。結(jié)果表明輸出信號(hào)的符號(hào)序列熵可有效反映共振狀況，能調(diào)節(jié)系統(tǒng)達(dá)到最佳信噪比輸出，而且易于工程實(shí)現(xiàn)，證實(shí)了該方法的有效性，進(jìn)而為衡量隨機(jī)共振提供一種新方法。

計(jì)量學(xué)；隨機(jī)共振；測(cè)度；信噪比；符號(hào)序列化；自適應(yīng)調(diào)節(jié)

1　引 言

信號(hào)檢測(cè)在機(jī)械故障診斷領(lǐng)域有著廣泛應(yīng)用。由于該領(lǐng)域的背景噪聲很強(qiáng)，而待檢信號(hào)非常微弱，檢測(cè)出有效信號(hào)非常困難［1］。傳統(tǒng)的信號(hào)檢測(cè)方法主要有相關(guān)檢測(cè)［2］、高階譜分析［3，4］、經(jīng)驗(yàn)?zāi)B(tài)分解［5］和卡爾曼濾波［6］等。這些方法主要通過抑制和去除噪聲來實(shí)現(xiàn)信號(hào)檢測(cè)，當(dāng)有用信號(hào)頻率和噪聲頻率混疊時(shí)，去除噪聲的同時(shí)對(duì)有用信號(hào)有一定損害，從而影響檢測(cè)精度。

隨機(jī)共振是Benzi［7，8］等提出的。與傳統(tǒng)方法不同，此方法是在非線性系統(tǒng)中，利用噪聲能量來增強(qiáng)有用信號(hào)，使得輸出信噪比增強(qiáng)。描述周期隨機(jī)共振的測(cè)度指標(biāo)有信噪比、信噪比增益與駐留時(shí)間分布。因?yàn)樵缙跈C(jī)械故障特征的周期性和調(diào)制性，其統(tǒng)計(jì)量呈現(xiàn)周期性變化，此時(shí)信號(hào)是循環(huán)的非平穩(wěn)的周期信號(hào)，而測(cè)度指標(biāo)信噪比和信噪比增益需要對(duì)特征信號(hào)有很好的估計(jì)。由于在實(shí)際中很難滿足，同時(shí)信噪比是以特征頻率處的幅值與同頻率的背景噪聲之比，所以存在局部性，不能從整體上衡量輸出信號(hào)。然而駐留時(shí)間分布卻很復(fù)雜，不易進(jìn)行量化，工程應(yīng)用中很難實(shí)現(xiàn)。針對(duì)以上問題，提出以二進(jìn)制符號(hào)化規(guī)則的符號(hào)序列熵作為衡量隨機(jī)共振輸出的測(cè)度，設(shè)計(jì)自適應(yīng)隨機(jī)共振系統(tǒng)進(jìn)行微弱信號(hào)檢測(cè)。

2　雙穩(wěn)態(tài)隨機(jī)共振模型

受隨機(jī)噪聲與外力驅(qū)動(dòng)的雙穩(wěn)態(tài)系統(tǒng)可以由Langevin方程描述，即

圖1　不同參數(shù)的勢(shì)函數(shù)

隨著系統(tǒng)參數(shù)a、b的變化勢(shì)高也發(fā)生變化，從而控制粒子在勢(shì)阱中來回躍遷。當(dāng)系統(tǒng)參數(shù)和噪聲強(qiáng)度達(dá)到某種匹配時(shí)，粒子在兩個(gè)勢(shì)阱中來回切換，切換速度快慢使得輸出信號(hào)與微弱周期信號(hào)同步，從而使系統(tǒng)輸出中的小周期分量得到增強(qiáng)。由于雙穩(wěn)態(tài)隨機(jī)共振模型的研究相對(duì)比較成熟，所以選取雙穩(wěn)態(tài)模型作為微弱信號(hào)檢測(cè)的方法。

3　符號(hào)序列熵隨機(jī)共振檢測(cè)原理

3.1符號(hào)序列分析方法

符號(hào)序列分析（Symbolic Sequence Analysis，SSA）通過將多個(gè)可能值的數(shù)據(jù)序列變成僅有幾個(gè)互不相同值的符號(hào)序列，是一個(gè)“粗?；边^程［9，10］。這一過程能夠捕捉具有確定性結(jié)構(gòu)的信號(hào)，尤其是周期信號(hào)特征，所以適合機(jī)械故障信號(hào)分析。對(duì)時(shí)間序列進(jìn)行符號(hào)統(tǒng)計(jì)量分析時(shí)，引入相空間離散劃分并分配符號(hào)給每個(gè)劃分單元，選取合適的閾值函數(shù)將時(shí)間序列轉(zhuǎn)化為符號(hào)序列。引入劃分p=｛p1，p2，…，pq｝，則通過式（3）的閾值過濾，從而將時(shí)間序列x（n）轉(zhuǎn)化為符號(hào)序列s（n）。

當(dāng)非線性系統(tǒng)、噪聲和有用信號(hào)趨于某種最佳匹配時(shí)，輸出信號(hào)中確定性結(jié)構(gòu)和周期性逐漸增強(qiáng)，經(jīng)過序列化后其序列結(jié)構(gòu)趨于穩(wěn)定和周期化。若故障信號(hào)無直流分量影響，則其輸出x值在0上下周期變化。因此使用一個(gè)一維空間的二進(jìn)制劃分來對(duì)x時(shí)間序列進(jìn)行符號(hào)化，從而獲得輸出信號(hào)的符號(hào)序列，它表征了原時(shí)間序列的變化規(guī)律，時(shí)間序列符號(hào)化方法如圖2所示。

圖2　時(shí)間序列符號(hào)化原理圖

1）符號(hào)序列編碼。選擇合適的時(shí)間延遲τ和符號(hào)短序列長(zhǎng)度L，實(shí)現(xiàn)對(duì)數(shù)據(jù)空間劃分，按式（4）將序列s（n）分組為短序列（L）：

式中，d是狀態(tài)空間的維數(shù)，q是劃分個(gè)數(shù)。

2）符號(hào)短序列編碼。對(duì)生成的每個(gè)短序列進(jìn)行編碼，方便符號(hào)序列的數(shù)值表示，一般根據(jù)式（5）把基n的符號(hào)序列變成基10的十進(jìn)制序列，則編碼完成。

3）符號(hào)序列香農(nóng)信息熵計(jì)算。每個(gè)符號(hào)短序列的概率反映所對(duì)應(yīng)時(shí)間序列的數(shù)據(jù)特點(diǎn)。采用香農(nóng)信息熵反映系統(tǒng)的無序化程度，定義改進(jìn)型香農(nóng)熵為：

式中，i是符號(hào)序列的序數(shù)，pi，L是長(zhǎng)度為L(zhǎng)的第i個(gè)符號(hào)序列的概率，若概率為零則剔除，從而克服了出現(xiàn)非數(shù)和無窮的情況。

改進(jìn)型香農(nóng)熵將隨符號(hào)序列長(zhǎng)度L的變化而變化。在實(shí)際中，使L從1開始增加，當(dāng)改進(jìn)型香農(nóng)熵Hs（L）有最小值時(shí)，此時(shí)確定的短符號(hào)序列長(zhǎng)度L為最佳；延遲步長(zhǎng)τ的選取過大導(dǎo)致相鄰符號(hào)間距增大，從而對(duì)有用信號(hào)的特征表達(dá)不夠充分；若延遲步長(zhǎng)τ的選取過小，信息表達(dá)過于冗余導(dǎo)致計(jì)算量增加，權(quán)衡二者適當(dāng)選取延遲τ即可。

對(duì)純?cè)肼曅盘?hào)而言，通過仿真試驗(yàn)發(fā)現(xiàn)其符號(hào)序列直方圖（如圖3（a）所示）中每個(gè)符號(hào)出現(xiàn)的概率幾乎相等，沒有明顯的主導(dǎo)概率，且根據(jù)式（6）計(jì)算其香農(nóng)熵為0.994 4，接近于1，說明無序化程度很高。

圖3　噪聲與周期信號(hào)的符號(hào)序列直方圖

相反對(duì)于周期的純凈信號(hào)，通過試驗(yàn)發(fā)現(xiàn)其符號(hào)序列直方圖（如圖3（b）所示）中每個(gè)符號(hào)出現(xiàn)的概率有很大差距，同時(shí)有很強(qiáng)的主導(dǎo)概率，其香農(nóng)熵為0.249 6，同時(shí)發(fā)現(xiàn)其編碼后的序列圖呈現(xiàn)周期性，說明其信號(hào)周期性很強(qiáng)，存在周期性結(jié)構(gòu)。最后通過在同一純凈信號(hào)中加入不同強(qiáng)度的噪聲得到的香農(nóng)熵值隨著噪聲強(qiáng)度的遞減而遞減，最終趨近于純凈信號(hào)的熵值。通過以上研究發(fā)現(xiàn)，可以將符號(hào)序列熵作為隨機(jī)共振輸出的測(cè)度，來調(diào)節(jié)隨機(jī)共振達(dá)到最佳輸出。

3.2自適應(yīng)隨機(jī)共振微弱信號(hào)檢測(cè)

定義符號(hào)序列熵差值如式（7）所示：

式中，Hs（i）是系統(tǒng)參數(shù)第i次迭代輸出信號(hào)的序列熵，H＾s是待檢信號(hào)的估計(jì)序列熵。隨著系統(tǒng)參數(shù)迭代，隨機(jī)共振趨于達(dá)到最佳的同時(shí)，其符號(hào)序列熵的差值ΔHs趨于零，通過判斷該差值來自適應(yīng)地調(diào)節(jié)系統(tǒng)參數(shù)，獲得最佳輸出信噪比。采用序列熵的隨機(jī)共振方法檢測(cè)微弱信號(hào)的流程如圖4所示，基本思路是：

圖4　微弱信號(hào)檢測(cè)流程圖

1）初始化。根據(jù)待測(cè)信號(hào)和背景噪聲強(qiáng)度估計(jì)參數(shù)a、b和H＾s，并對(duì)a、b置初值，盡量使勢(shì)壘高度較小，以便后續(xù)調(diào)節(jié)；設(shè)定符號(hào)序列的劃分線p，本文設(shè)定p為0。

2）時(shí)間序列的符號(hào)化。將采樣數(shù)據(jù)加入雙穩(wěn)態(tài)系統(tǒng)中，采用四階龍庫塔法求解朗之萬方程得到輸出信號(hào)，計(jì)算步長(zhǎng)取信號(hào)采樣周期的倒數(shù)，獲得時(shí)間序列｛xi｝并將其離散化為序列｛x（n）｝；設(shè)定合適的延遲步長(zhǎng)τ，并根據(jù)設(shè)定閾值p，通過閾值函數(shù)過濾｛x（n）｝，將序列轉(zhuǎn)化為符號(hào)序列s（n）。

3）計(jì)算改進(jìn)型香農(nóng)信息熵。確定符號(hào)短序列長(zhǎng)度L，工程上一般取3～8，本文取L=5，τ=1生成序列｛（L）｝；對(duì)符號(hào)序列｛（L）｝進(jìn)行序列編碼，形成基10的十進(jìn)制序列s（L）；計(jì)算符號(hào)序列訪問每個(gè)單元的概率pi，L，并根據(jù)式（6）計(jì)算香農(nóng)信息熵Hs。

4）判斷待測(cè)信號(hào)中有無特定信號(hào)，若有則自適應(yīng)調(diào)節(jié)系統(tǒng)達(dá)到最佳共振。如果Hs≈1說明待測(cè)信號(hào)中無特定信號(hào)，否則含有特定信號(hào)，則根據(jù)ΔHs（i）自適應(yīng)調(diào)節(jié)系統(tǒng)參數(shù)a和b，遵循以下調(diào)整原則：若雙穩(wěn)態(tài)系統(tǒng)輸出信號(hào)的熵值較大，說明噪聲占主要成分，按步長(zhǎng)h調(diào)整系統(tǒng)參數(shù)a和b增加勢(shì)壘高度，直到系統(tǒng)輸出的熵值接近估計(jì)值為止，此時(shí)輸出為最佳共振輸出；若系統(tǒng)輸出熵值小于估計(jì)值，則調(diào)整參數(shù)a，b使勢(shì)壘降低，讓系統(tǒng)輸出熵值接近估計(jì)熵值為止。

4　仿真驗(yàn)證及分析

實(shí)驗(yàn)1：對(duì)單頻信號(hào)x（t）=s（t）+n（t），其中s（t）=2 sin（0.18πt），n（t）是均值為0、強(qiáng)度為9的高斯白噪聲。采樣頻率20 Hz，采樣時(shí)間50 s，通過估計(jì)得到信號(hào)的估計(jì)序列熵H＾s=0.2535，系統(tǒng)參數(shù)a、b置初值為0.1和1，迭代步長(zhǎng)h為0.1，b不變，進(jìn)行仿真，根據(jù)序列熵尋找的最佳信噪比輸出為a =2時(shí)，輸出信噪比為1.094 2 dB，序列熵Hs= 0.2604，則ΔHs=Hs-H＾s＜ε=0.01滿足檢測(cè)要求，跳出迭代循環(huán)。此時(shí)雙穩(wěn)態(tài)系統(tǒng)輸出時(shí)頻圖如圖5所示。

圖5　系統(tǒng)輸出信號(hào)時(shí)頻圖

從圖5（a）可以看出，以序列熵為測(cè)度自適應(yīng)檢測(cè)微弱信號(hào)的方法能很好地保留信號(hào)的趨勢(shì)，圖5（b）中有效提取出特征頻率0.09 Hz，為周期信號(hào)的檢測(cè)開創(chuàng)了一條新途徑。

實(shí)驗(yàn)2：多頻信號(hào)x（t）=x1（t）+x2（t）+n（t），其中x1（t）=0.2 sin（2π×0.05t），x2（t）= 0.3 sin（2π×0.1t），n（t）是均值為0、強(qiáng)度為5的高斯白噪聲。采樣頻率20 Hz，采樣時(shí)間50 s，通過經(jīng)驗(yàn)估計(jì)得到特征信號(hào)的估計(jì)序列熵H＾s=0.260 4，系統(tǒng)參數(shù)a和b置初值0.01和1，迭代步長(zhǎng)h為0.01，b不變，當(dāng)a=1.01時(shí)，序列熵Hs=0.277 1，ΔHs=Hs-H＾s=0.016 7＜ε（這里選擇ε=0.2，由于粗迭代使得下一個(gè)熵值小于估計(jì)熵值，信息出現(xiàn)丟失）則滿足檢測(cè)要求，跳出迭代循環(huán)。此時(shí)輸出信噪比-8.022 3 dB，雙穩(wěn)態(tài)系統(tǒng)輸出時(shí)域波形如圖6（a）所示。從圖中可以看出，序列熵為測(cè)度能有效自適應(yīng)地調(diào)節(jié)系統(tǒng)參數(shù)，獲得最佳系統(tǒng)輸出。本文只通過調(diào)節(jié)系統(tǒng)參數(shù)a來獲取最佳輸出，同樣也可以調(diào)節(jié)系統(tǒng)參數(shù)b或者同時(shí)調(diào)節(jié)二者來快速地得到最佳輸出。圖6（b）中可以清晰看出特征頻率0.05 Hz和0.1 Hz，文中采用了頻率小于1 Hz的信號(hào)進(jìn)行仿真，針對(duì)大頻率可以采用尺度變換結(jié)合熵測(cè)度來實(shí)現(xiàn)，由于篇幅所限不再贅述。

圖6　系統(tǒng)輸出信號(hào)時(shí)頻圖

5　工程應(yīng)用

采用美國(guó)凱斯西儲(chǔ)大學(xué)滾動(dòng)軸承數(shù)據(jù)中心的實(shí)驗(yàn)數(shù)據(jù)［11］，驅(qū)動(dòng)端電機(jī)軸由包含故障的JEM SKF 6025-2RS深溝球軸承支撐。采用驅(qū)動(dòng)端內(nèi)圈故障數(shù)據(jù)（編號(hào)：213.mat），在內(nèi)圈表面人為加工了直徑為0.533 4 mm的凹坑，采樣頻率48 kHz，采樣時(shí)間0.5 s，軸的轉(zhuǎn)速為1 797 r/min，軸承主要參數(shù)見網(wǎng)站。軸回轉(zhuǎn)頻率為frotate=29.16 Hz，由滾動(dòng)軸承故障特征頻率理論計(jì)算得到內(nèi)圈、外圈、滾珠故障頻率分別為finner=162.19 Hz，fouter=107.36 Hz和froller= 141.17 Hz。圖7給出了試驗(yàn)軸承內(nèi)圈故障信號(hào)的時(shí)域波形，可以看出軸承內(nèi)圈故障具有一定的沖擊特性和大量的噪聲，無規(guī)律可尋。

圖7　內(nèi)圈故障原始信號(hào)時(shí)域波形

受絕熱理論的限制，隨機(jī)共振無法檢測(cè)頻率大于1 Hz的信號(hào)，所以首先對(duì)軸承內(nèi)圈信號(hào)選擇尺度因子1 000，進(jìn)行歸一化尺度變換［12］使其滿足隨機(jī)共振的輸入條件，進(jìn)而采用本文方法進(jìn)行自適應(yīng)隨機(jī)共振檢測(cè)滾動(dòng)軸承故障特征，按迭代步長(zhǎng)調(diào)整系統(tǒng)參數(shù)使系統(tǒng)輸出序列熵值接近估計(jì)值0.519，系統(tǒng)輸出信號(hào)的幅值譜如圖8所示，此時(shí)系統(tǒng)參數(shù)a、b分別為1 000，950?？梢钥闯鲚S承回轉(zhuǎn)頻率29.3 Hz（理論值約為29.16），其倍頻信息（58.6≈2frotate，87.9≈3frotate等），外圈故障頻率106.9 Hz（理論值為107.36），滾珠故障頻率146.5 Hz（理論值約為141.17）以及內(nèi)圈故障161.1 Hz（接近162.19），除此之外還有明顯的倍頻和邊帶頻信息，與原始滾動(dòng)軸承內(nèi)圈故障的直接幅值譜（圖9所示）進(jìn)行對(duì)比，可以看出故障特征頻率在某種程度上有很大提高，同時(shí)噪聲頻率有效地被削弱，隨機(jī)共振對(duì)不同頻率信號(hào)增強(qiáng)程度是不同的，通過模擬仿真和工程應(yīng)用驗(yàn)證了該方法的有效性。

圖8　隨機(jī)共振系統(tǒng)輸出信號(hào)頻譜

圖9　滾動(dòng)軸承內(nèi)圈故障信號(hào)頻譜

6　結(jié) 論

選取迭代步長(zhǎng)時(shí)可以進(jìn)行“粗迭代”，到靠近估計(jì)熵值時(shí)，進(jìn)行“細(xì)迭代”，從而有效提高迭代速率和確保更好的精度。熵值過小表示其信號(hào)中所含信息量較少，所以在選取差值時(shí)，不可過小，同時(shí)所得到的熵值必須稍大于估計(jì)值，從而能保留更多特征信息。符號(hào)序列熵需要先驗(yàn)知識(shí)的支撐來得到估計(jì)熵值，所以在一定范圍內(nèi)存在局限性。下一步工作是研究如何有效估計(jì)熵值。

［1］ 謝平，楊玉昕，江國(guó)乾，等.基于局部均值分解的滾動(dòng)軸承故障診斷新方法［J］.計(jì)量學(xué)報(bào)，2014，35（1）：73-77.

［2］ Dong H B，Qi K Y，Chen X F，et al.Sifting process of EMD and its application in rolling element bearing fault diagnosis［J］.JournalofMechanicalScienceand Technology，2009，23（8）：2000-2007.

［3］ 楊江天，陳家驥，曾子平.基于高階譜的旋轉(zhuǎn)機(jī)械故障征兆提?。跩］.振動(dòng)工程學(xué)報(bào)，2001，14（1）：14-17.

［4］ 劉雪霞，張琦，譚業(yè)發(fā).高階譜分析技術(shù)在軸承故障信號(hào)特征提取中的應(yīng)用［J］.機(jī)床與液壓，2008，36（7）：155-160.

［5］ 臧懷剛，李清志，韓艷龍，等.基于改進(jìn)的經(jīng)驗(yàn)?zāi)B(tài)分解的滾動(dòng)軸承故障診斷研究［J］.計(jì)量學(xué)報(bào)，2013，34（2）：101-105.

［6］ 蘇理云，馬洪，唐世福.用混沌振子和Kalman濾波檢測(cè)強(qiáng)分形噪聲中的弱信號(hào)［J］.四川大學(xué)學(xué)報(bào)（工程科學(xué)版），2007，39（3）：149-154.

［7］ Benzi R，SuteraA，VulpianiA.Themechanismof stochastic resonance［J］.JournalofPhysicsA：Mathematical and General，1981，14（11）：453-457.

［8］ Benzi R，Parisi G，Sutera A，et al.A theory of stochastic resonance in climatic change［J］.SIAM Journal on Applied Mathematics，1983，43（3）：565-578.

［9］ 王清亮.基于符號(hào)序列信息熵混沌特性的微弱信號(hào)檢測(cè)［J］.儀器儀表學(xué)報(bào)，2009，30（12）：2491-2497.

［10］ 張愛華，王平，丑永新.基于動(dòng)態(tài)差分閾值的脈搏信號(hào)峰值檢測(cè)算法［J］.吉林大學(xué)學(xué)報(bào)（工學(xué)版），2014，44（3）：.847-853.

［11］ Case Western Reserve University Bearing Data Center Website［EB/OL］.http：//www.eecs.case.edu/ laboratory/bearing/download_fan.htm

［12］ 陳敏，胡蔦慶，秦國(guó)軍，等.參數(shù)調(diào)節(jié)隨機(jī)共振在機(jī)械系統(tǒng)故障早期檢測(cè)中的應(yīng)用［J］.機(jī)械工程學(xué)報(bào)，2009，45（4）：131-135.

Weak Signal Detection of Adaptive Stochastic Resonance Based on Shannon Entropy of Symbolic Series

PAN Zheng-rong， QIAO Zi-jian， ZHANG Ning
（College of Electrical and Information Engineering，Lanzhou University of Technology，Lanzhou，Gansu 730050，China）

A method of weak features detection for mechanical fault data is proposed，which combined the stochastic resonance and the improved Shannon entropy of symbolic series.Firstly，the basic principles of SR and symbolic sequence analysis（SSA）is briefly introduced.Secondly，symbolic sequence Shannon entropy of output signal is the measure whether the state of this system of SR is the best one.Finally，according to adjust system parameters a and b adaptively，the system of SR achieves the best state and analyzes signal to noise ratio（SNR）and spectrum of output signal.Information of output signal show that Shannon entropy of symbolic series can reflect the state of SR，the best SNR of output signal can be achieved and it can be easily used in practical projects，which confirms the new method proposed is correct and useful.

metrology；SR；measure；SNR；SSA；adaptive regulation

TB93

A

1000-1158（2015）05-0496-05

10.3969/j.issn.1000-1158.2015.05.11

2013-12-23；

2014-08-28

甘肅省省自然科學(xué)研究基金（1308RJZA273）

潘崢嶸（1964-），男，山東蓬萊人，蘭州理工大學(xué)教授，主要從事計(jì)算機(jī)控制、智能檢測(cè)與控制。譙自健為本文通訊作者。zijianqiao@hotmail.com