指數型隨機共振微弱振動信號檢測方法

張 剛,曹 莉,賀利芳,易 甜

(1.重慶郵電大學 通信學院，重慶 400065;2.信號與信息處理重慶市重點實驗室，重慶 400065)

微弱信號指自身強度很小并且淹沒在噪聲中一類信號，微弱信號檢測即為采用各種技術手段提高信噪比將淹沒在強噪聲中的微弱信號提取出來。一種技術手段針對噪聲，通過噪聲與信號之間差異抑制噪聲，例如濾波技術、時頻分析技術、相關檢測技術，其應用涉及到實際生產的各個方面，如故障信號、地震信號、生物電信號，這些技術已經成為廣大學者的研究熱點[1]。另一種技術為增強信號，將噪聲與信號通過非線性方程進行能量轉化，將噪聲能量轉化為信號能量使得信號突出被檢測。隨機共振(Stochastic Resonance,SR)最早為1981年物理學家Benzi等[2]在對地球冰川期氣候變化的研究中首次提出；1983年Fauve等[3]在施密特觸發器實驗中觀察到了隨機共振現象。隨著各領域間隨機共振發現，人們對隨機共振現象機理研究越來越火熱，并取得豐富研究成果[4-5]。

隨機共振微弱信號檢測技術在微弱信號檢測及提取方面應用廣泛，相比較于傳統的抑制噪聲技術，在保留信號信息的基礎上，最大化將噪聲能量轉化為信號能量。常采用一維Langevin方程模型即傳統雙穩[6](Classic Bistable Stochastic Resonance,CBSR)系統，在參數滿足絕熱近似理論條件下實現信號檢測，為了更好地利用噪聲，一些新的模型陸續被提出并被研究，Qiao等[7]通過構造分段雙穩從勢函數勢阱峭度研究勢阱峭度對隨機共振輸出信號飽和度影響；周玉飛等[8]使用級聯雙穩研究能量多次轉換問題；陸思良等[9]采用兩個一維隨機共振組合而成的二維互補隨機共振應用于軸承信號檢測中，這些模型大多基于一維非線性方程。隨著研究的深入發現二階非線性方程也能產生隨機共振現象。Duffing振子方程是一類典型的二階非線性方程不僅能通過混沌理論檢測微弱信號，也能夠產生隨機共振現象。Duffing振子隨機共振最早提出在20世紀80年代[10-11]，隨后Almog等[12-13]搭建硬件平臺實現Duffing振子隨機共振信號檢測，Leng等[14]研究不同參數調節Duffing振子隨機共振對輸出信號影響；Lai等[15]研究在Duffing方程下噪聲誘導隨機共振發生；Li[16]將Duffing方程隨機共振應用在MSK信號中；Zhang等[17]提出的耦合雙穩Duffing隨機共振系統不僅能提高輸出信噪比，有效檢測正弦信號及方波信號。

Duffing振子隨機共振擴展了隨機共振理論與實踐研究，本文提出指數型雙穩作用在Duffing方程隨機共振系統，推導其逃逸率并研究阻尼參數對隨機共振輸出信噪比影響；采用衰減沖擊信號和振動諧波信號進行驗證，證明所提系統有增強信號作用；并且提出隨機共振與經驗模態分解結合下信號檢測方法，應用于軸承故障信號檢測，研究結果表明有效的檢測到故障頻率以及故障的二倍頻信號、三倍頻信號。

1 勢函數模型

1.1 簡諧勢阱模型

簡諧勢阱為線性隨機共振系統原始模型，后來演化為冪函數型勢阱。系統結構簡單參數單一，其表達式[18]為

(1)

式中，a為系統參數，不同系統參數a下勢函數圖形如圖1所示。a越大勢阱越陡峭、勢阱底部中心寬度越窄粒子運動穿過勢阱底部時間越短。調節系統參數、信號參數、噪聲使三者協同產生隨機共振。簡諧勢阱只有一個系統參數不用考慮系統參數間的相互影響，使之成為研究熱點之一。

圖1 簡諧勢函數Fig.1 Harmonic potential function

1.2 Gaussian potential勢阱模型

在核物理學中，常用GP模型來描述復核散射，徑向GP模型[19]可以表示為

(2)

式中，V表示勢阱深度，R表示勢阱寬度。圖2表示GP勢阱的勢函數，由圖可知，勢阱的兩端迅速收斂于0。由圖2(a)可知，固定V，單獨調節R可以改變勢阱的寬度，并且隨著R的減小，勢阱寬度逐漸減小，勢阱壁逐漸變得陡峭。圖2(b)顯示，GP勢的勢阱深度由V唯一表征。

(a)V=1

(b)R=1

1.3 組合模型

基于以上兩種模型，將簡諧勢阱模型和GP勢阱模型相結合，提出一種指數型雙勢阱模型，勢函數如下

(3)

圖3為結構參數為：a=1，V=1，R=0.5的指數型雙穩勢函數。勢阱函數的兩個勢阱是對稱的，每個勢阱的寬度和勢阱壁的陡峭度均可通過參數a，V，R進行調節。在簡諧勢阱模型中加入GP勢阱使其變成雙勢阱，使振蕩粒子由在單一勢阱中運動變成在兩勢阱間進行躍遷，提高噪聲的利用率從而改善輸出信噪比，得到較好的隨機共振效果。從該新型系統的勢函數來看，通過調節系統參數，可以實現模型在單穩態勢阱和雙穩態勢阱之間的切換。因此，該新型模型兼具單穩勢阱模型與雙穩勢阱模型的優點。

圖3 指數型雙穩勢函數Fig.3 The exponential bistable potential function

2 隨機共振

2.1 隨機共振機理

隨機共振模型為

(4)

粒子逃逸率反映布朗粒子運動規律，噪聲引起布朗粒子躍遷Kramers逃逸率[20]為

(5)

(6)

布朗粒子在某一勢阱中的平均駐留時間TK=1/rk與勢函數變化時間相等才使得隨機共振發生

TK=T/2

(7)

(8)

由于噪聲太大信號被淹沒,從波形圖和頻譜圖中無法識別，如圖4(a)和4(b)所示。經過隨機共振后輸出波形接近方波如圖4(c)所示，0.01 Hz信號頻譜幅度為87，如圖4(d)所示，相對于周圍噪聲幅度信號突出，說明指數型雙穩隨機共振使信號能量大幅度提高。

2.2 系統參數對雙穩隨機共振影響

傳統Duffing雙穩隨機共振在噪聲系數過大下，增大阻尼比k，達到調節粒子躍遷速度使之與勢函數周期變化趨勢同步的目的，產生共振現象。指數型Duffing方程隨機共振條件為式(8)，在滿足絕熱近似理論條件下定義函數F表達式如下

F(a,u,V,D,k)=

(a)輸入信號波形

(b)輸入頻譜

(c)輸出信號波形

(d)輸出頻譜

(9)

固定一組參數采用數值仿真，信號采用幅度0.1，頻率0.01 Hz正弦周期信號，在采樣頻率5 Hz下采樣5 000點，a=1.74，V=1，R=0.8不同k值下仿真500次取平均SNR曲線圖如圖5(a)所示。由圖5(a)可知隨著D的增大輸出信噪比出現先增大后減小趨勢，不同阻尼系數k在同一噪聲強度對應輸出信噪比值一樣大，D=0.31的為隨機共振輸出最高點，此時輸出信噪比為-9.825 dB，相對于輸入信號(信噪比-31.416 9 dB)信噪比提高21 dB。固定信號參數、采樣頻率、采樣點、噪聲強度D=0.31、V=1不變。R=0.8時參數a對輸出信噪比的影響如圖5(b)所示，隨著a的增加，輸出信噪比出現先增后減最后趨于平穩趨勢。在較小a的情況下，雙穩勢函數勢阱邊緣平坦，粒子躍遷無法通過，即使存在噪聲也無法驅動粒子完成躍遷。隨著a的增大布朗粒子運動在一個勢阱內局部運動，此時系統沒達到躍遷條件，部分信號能量轉化為噪聲能量導致輸出信噪比降低。a繼續增大，系統條件轉化為對信號有利條件，噪聲能量逐漸轉化為信號能量，當a達到某一值時達到躍遷條件隨機共振現象發生。固定a=1.75，圖5(c)為參數R對輸出信噪比影響曲線，隨著R的逐漸增加，輸出信噪比先減后增大到極大值后再減小，隨著R的增大勢壘高度逐漸降低，勢壘高度達到布朗粒子運動躍遷高度，達到隨機共振發生條件此時三者協同作用下噪聲能量最大化轉化為信號能量。由此可通過調節噪聲強度以及系統參數實現隨機共振。

2.3 隨機共振參數優化

當信號噪聲固定時，需要調節系統參數達到最佳共振效果，單一的調節會忽略參數之間影響。目前隨機共振的衡量指標主要有功率譜放大系數、信噪比增益、相關系數等等。其中，信噪比增益更能直觀地反應隨機共振效果，本文采用噪比增益作為衡量指標，其公式為[21]

(10)

式中，SNRout為輸出信噪比，SNRin為輸入信噪比。SNRgain大于1時表明隨機共振系統對輸入待測弱信號有明顯的改善和增強作用。本文采用粒子群算法以信噪比增益為適應度函數，對每一組實現數據進行參數選擇。

(a)D

(b)a

(c)R

步驟1 種群初始化。固定信號幅度，噪聲強度，及阻尼比固定為0.5，設置種群數量，參數a，V，R搜索范圍以及最大進化代數Tmax，參數最大搜索速度取最大調整步長的10%～25%，這里的最大調整步長指的是所設定的粒子位置范圍的上限值減去粒子位置范圍的下限值所得的差值。隨機初始化搜索點的位置，并計算出其相應的個體極值，記錄整個粒子群中個體極值最大的粒子序號，設置Nbest為該最大粒子的當前位置。

步驟2 評價每個粒子。適應度函數采用SNRgain公式，計算粒子適應度值。與該粒子當前個體極值進行比較，若大于后者，則更新粒子個體極值，若在該粒子的鄰域內所有粒子的個體極值中最大的大于當前的Nbest，則設置Nbest為該粒子的位置，記錄該粒子的序號，并更新Nbest的函數值。

步驟3 粒子的更新,采用粒子群算法迭代公式[22]更新粒子速度及位置。

步驟4 檢驗是否符合結束條件。判斷當前的迭代次數是否達到最大進化代數Tmax，滿足條件則停止迭代，并輸出最優參數，否則轉至步驟2。

步驟5 檢測結果。根據對a，V，R優化輸出的最優解，計算適應度值，輸出最優a，V，R。

3 經驗模態分解

經驗模態分解(Empirical Mode Decomposition，EMD)是用來實現降低噪聲的方法之一，通過數據特征時間尺度來獲得固有波動模式，從而自適應選擇函數實現信號分解。經驗模態分解后可以明顯的得到頻率大、小的信號，頻率大的信號優先被分解得出，頻率低的隨后被分解，所有的IMF分量疊加得到的信號為最原始信號。分解過程中產生分解量為模態函數(IMF)。IMF信號分量中，信號極值點和過零點數目相等或相差1，為了得到滿足條件模態函數，需要多次篩選，具體算法如下[23]

步驟1 計算信號局部極大值極小值，通過3次樣條插值法擬合出上下包絡線平均值m1(t)，并且認為h1(t)=x(t)-m1(t)為殘余分量。

步驟2 理想情況下h1為第一個IMF分量，判斷h1是否為滿足IMF分量，不滿足需要反復篩分，接下來將h1作為新信號，重復上述步驟，循環k次后，得到IMF條件h1k(t)。其中篩選次數約束滿足柯西準則[24]

(11)

式中，T為信號時間長度，ε為門限值范圍為0.2～0.3。

步驟3 得到第一階IMF的c1(t)即為h1k(t)，r1(t)=x(t)-c1(t)將c1(t)作為原始信號，反復重復上面兩個步驟，得到c2(t),c3(t),…,cn(t)，和剩余分量rn(t)，分解結束條件為rn(t)單調，由此可以將信號分解為n個模態分量。

經驗模態分解對于小波分解優點在于不需要預先設定基函數，只需根據信號自身特性進行平滑處理，最后達到降低噪聲，檢測出信號的效果。

4 仿真信號分析

為了驗證本文所提指數型雙穩隨機共振系統在微弱信號檢測中的效果，分別構造沖擊衰減信號和諧波振動信號進行模擬實際工程仿真，更加貼近實際情況。

4.1 沖擊衰減信號

機械運動軸承信號產生都是具有沖擊性，這種信號淹沒在噪聲中更加難以發現，為模擬沖擊震蕩信號構造沖擊衰減信號

sn(t)=H(t)×exp[-100(t-t/T/T)]×sin[1 200π(t-t/T/T)]+n(t)

(12)

式中，H(t)為階躍函數，T=0.128 s是沖擊周期，n(t)為均值為0方差為0.6的高斯噪聲，采樣點為2 048，采樣頻率fs=2 000 Hz，特征信號頻率為f=1/T=7.8 Hz，固定H(t)=0.2構造周期沖擊衰減信號如圖6所示。圖中顯示周期沖擊衰減信號加入噪聲后無法分辨出信號，將含噪信號送入指數型雙穩隨機共振系統中，特征信號頻率為大參數需要預處理變為小參數，進行二次采樣處理大頻率信號，二次采樣頻率為fsr=5 Hz，采樣壓縮比為R=fs/fsr=400，相當于待測信號壓縮為0.018 75 Hz滿足小參數條件。

(a)波形圖

(b)頻譜圖

經過參數為(a,b)=(1,1)經典雙穩隨機共振后特征信號7.8 Hz頻譜幅度為27.96，周圍噪聲依然很大，信號不能識別如圖7所示。由于阻尼系數對隨機共振影響不大，固定k=0.5使用3.3節尋優參數找到(a,V,R)=(1.528 3,3.512 7,0.151 8)進行指數型隨機共振，將輸出信號進行大頻恢復，特征信號7.8 Hz頻譜幅度為44.9，如圖8所示，相對于周圍噪聲信號成分突出明顯，證明此系統在振動信號中檢測有用性。

圖7 經典雙穩隨機共振輸出Fig.7 The output signal spectrum after CBSR

4.2 諧波振動信號分析

諧波振動信號是一種典型的多頻振動信號，通過分析這種振動信號，可以使我們準確了解機器的運行狀態。構造一種包含兩種諧波振動成分的含噪信號

圖8 指數型隨機共振輸出Fig.8 The output signal spectrum after exponential SR

sn(t)=A1sin(2π×ft)+A2sin(2π×2ft)+n(t)

(13)

圖9 含噪信號Fig.9 Signal with noise

將二次采樣后含噪信號送入指數型雙穩系統，固定k=0.5使用3.3節尋優參數找到(a,V,R)=(1.687 4,0.470 4,0.050 8)進行指數型隨機共振，隨機共振后信號頻譜如圖11所示，高頻噪聲降低，低頻成分升高，振動信號成分能量增強。圖11可知振動成分f頻譜幅度為203.1、2f頻譜幅度為176.8。上述兩實驗分別驗證兩種信號下指數型隨機共振在信號檢測方面優于經典雙穩系統。

含噪信號在指數型雙穩隨機共振后信號能量得到較大增加，但周圍依然存在大量噪聲。經驗模態分解是很好的降噪技術手段，能將信號的成分從高頻到低頻逐層分解，為了有效的檢測到待測的多頻信號，將指數型雙穩隨機共振輸出信號進行EMD分解得到IMF成分如圖12所示，IMF1～IMF3含有幅度太小的噪聲，IMF4～IMF5含有f信號成分，IMF6含有2f信號成分IMF7～IMF8為剩余分量。將IMF4～IMF6含有信號部分疊加得到降低噪聲后的合成信號如圖13所示，可以明顯的觀察到一次諧波與二次諧波信號，這相對于圖11更加準確的檢測到信號。

圖10 經典雙穩隨機共振輸出Fig.10 The output signal spectrum after CBSR

圖11 指數型隨機共振輸出Fig.11 The output signal spectrum after exponential SR

圖12 EMD分解后IMFs頻譜分量Fig.12 The spectrum of IMFs decomposed by EMD

因此在多頻信號大噪聲環境下，可以先通過隨機共振進行部分能量轉移，再將輸出信號進行降噪處理提取待測信號成分進行合成，能得到比隨機共振后更加明顯信號成分。

圖13 合成信號頻譜Fig.13 The synthetic signal spectrum

5 工程應用

滾動軸承故障信號具有非平穩性、調制性、微弱性，由于常常淹沒在強大的背景噪聲中難以發現和提取[25]。軸承故障出現對工程應用有不可估量的損失常常發生在內圈、外圈、滾動體等。軸承故障信號存在往往以倍頻形式，故障越嚴重倍頻成分越多。本文提出基于指數型雙穩隨機共振系統與EMD聯合檢測微弱信號方法并對軸承的內圈故障信號進行驗證。首先，用二次采樣對不滿足絕熱近似理論的故障信號進行預處理；其次，在將信號送入隨機共振系統后，參數尋優后找到最佳參數對信號進行隨機共振；再次，將隨機共振后信號通過經驗模態分解，觀察故障頻率、二倍頻信號和三倍頻信號；最后，含有故障信號及其倍頻信號的模態分量進行合成輸出，實現故障信號檢測。本文故障信號采用Case Western Reserve University電氣工程實驗室[26]公開數據，選用型號為6205-2RS JEM SKF的深溝球軸承為實驗對象，其主要參數如表1所示。軸承轉速r=1 796 r/min(29.93 Hz)，采樣頻率fs=12 kHz，二次采樣頻率5 Hz，采樣點數N=6 000，軸承的故障特征頻率如表2所示。外圈的故障特征頻率為162.2 Hz，其二倍頻為324.4 Hz。三倍頻率為486.6 Hz。

表1 滾動軸承主要計算參數Tab.1 The main computation parameters of the rolling bearing

表2 滾動軸承故障特征頻率Tab.2 Rolling bearing fault feature frequency

故障信號時域圖與頻譜圖如圖14(a)所示，軸承轉動由于故障之處產生周期性沖擊信號，2 000～4 000頻段存在干擾信號，故障信號存在低頻段。將信號二次采樣預處理后送入參數為(a,V,R,k)=(2.530 1,3.524 2,1.758 1,0.5)指數型雙穩隨機共振系統，隨機共振后輸出信號頻譜如圖14(b)所示，特征160 Hz顯現為故障頻率(誤差在允許范圍內)，但其二倍頻、三倍頻特征沒有顯現。將隨機共振后輸出信號進行EMD分解如圖15所示，IMF4～IMF3成分含有特征信號，IMF3～IMF2含有二倍頻信號，IMF2含有三倍頻信號。圖16為IMF2～IMF4分量進行疊后合成信號，由此可見大噪聲存在時，隨機共振都無法識別多頻成分，結合經驗模態分解后可使沒被檢測出的信號得到檢測。

(a)輸入內圈故障信號

(b)指數型隨機共振輸出信號

圖15 EMD分解后IMFs頻譜分量Fig.15 The spectrum of IMFs decomposed by EMD

圖16 合成信號Fig.16 The synthetic signal

6 結 論

本文將簡諧勢阱與GP勢阱系統組合成一種新型的指數雙穩系統，將此勢阱模型應用在Duffing振子方程，進行隨機共振研究。

(1)在小參數條件下阻尼系數k對隨機共振輸出信噪比影響不大，在信號一定條件下存在某一噪聲強度使得隨機共振發生；

(2)隨著a(或R)增加，輸出信噪比大致趨勢都是先增后減，因此調節參數可以達到最佳共振效果；

(3)為了證明新型指數型雙穩隨機共振系統對弱信號檢測的有用性，構造沖擊衰減信號以及諧波振動信號，系統參數選擇采用粒子群算法，以輸出噪比增益為目標函數，找到最優參數(a,V,R)分別驗證了所提系統的有用性；

(4)為了更好的檢測出多頻諧波信號提出一種隨機共振后經驗模態分解信號檢測方法，并將所提方法應用在軸承故障信號檢測中。實驗結果表明，通過將隨機共振輸出信號進行EMD處理，可得到除了故障頻率，還有故障頻率倍頻信號。這種先通過隨機共振后進行EMD結合兩者的優勢，能更加準確的檢測故障信號。