劉運江， 王輔忠， 劉 露

(天津工業(yè)大學(xué)理學(xué)院 天津，300387)

引 言

Benzi等[1]提出了隨機共振的概念，現(xiàn)已成為強噪聲背景下弱信號提取的重要工具之一。20世紀(jì)90年代，Hu等[2]提出的絕熱近似理論、文獻[3]提出的線性響應(yīng)理論和文獻[4]提出的駐留時間分布理論為經(jīng)典隨機共振理論奠定了理論基礎(chǔ)，標(biāo)志著經(jīng)典隨機共振理論由誕生走向成熟。隨機共振已從經(jīng)典隨機共振理論發(fā)展到了非線性理論[5]、大參數(shù)調(diào)節(jié)理論[6]和多穩(wěn)態(tài)理論[7]等非經(jīng)典隨機共振理論。

近幾年，以Levy噪聲作為背景噪聲在隨機共振提取微弱信號領(lǐng)域的研究備受關(guān)注。Levy分布噪聲[8]服從中心極限定理和大數(shù)定理，保留了自然隨機噪聲的特征特點，是更為普遍的噪聲，但由于Levy噪聲的函數(shù)表達式為超越方程，無法對其進行直接求解，因此對Levy噪聲的研究相對于高斯白噪聲比較緩慢。Zolotarey等[9]對Levy噪聲的分布函數(shù)進行求解和證明，使Levy噪聲得到了完整的解析和表達。賀利芳等[10]對Levy噪聲背景下冪函數(shù)型雙穩(wěn)態(tài)系統(tǒng)微弱信號的提取進行了研究。文獻[11]研究了Levy噪聲下一類周期勢系統(tǒng)的振動共振。張剛等[12-13]研究了Levy噪聲驅(qū)動下的冪函數(shù)型單穩(wěn)態(tài)和指數(shù)型單穩(wěn)態(tài)隨機共振的特性。

沖擊信號是電機、軸承和齒輪等機械設(shè)備中常見的故障信號，包含了設(shè)備重要的狀態(tài)信息，對這些沖擊信息的提取在一定程度上可以預(yù)防機械設(shè)備故障的發(fā)生和確保產(chǎn)品的合格生產(chǎn)，對工業(yè)生產(chǎn)來講具有重要意義。因此，強噪聲背景下弱沖擊信號的檢測是當(dāng)前研究的熱點之一。Li等[14]研究了自適應(yīng)級聯(lián)隨機共振系統(tǒng)對齒輪故障信號的特征提取。Lei等[15]利用欠阻尼隨機共振系統(tǒng)對早期軸承故障信號進行了提取。Liu等[16]對新型非線性系統(tǒng)軸承故障信號的提取進行了研究。這些文獻中利用的信號模型均為沖擊信號。

以上研究均以高斯白噪聲作為背景噪聲，忽略了實際噪聲的隨機脈沖性和復(fù)雜性。筆者提出了利用Levy噪聲作為背景噪聲，研究沖擊信號的非對稱雙穩(wěn)態(tài)隨機共振機理，利用KI作為信號處理的衡量標(biāo)準(zhǔn)，實現(xiàn)強噪聲背景下弱沖擊信息的提取。實驗結(jié)果表明，該方法不僅可以將機械故障沖擊信號的頻域特征信息成功提取出來，還可以獲取其時域特征信息。

1 Levy噪聲模型

Levy噪聲服從alpha穩(wěn)定分布理論，alpha穩(wěn)定分布又稱非高斯穩(wěn)態(tài)分布和重尾分布，是唯一滿足廣義中心極限定理的分布，其拖尾以平方律衰減。Levy噪聲標(biāo)準(zhǔn)參數(shù)系下的特征函數(shù)表達式[17]為

(1)

其中：參數(shù)α，β，σ和μ的取值分別為[0,2]，[-1,1]，[0,1]和[-∞.TIF,+∞]；α為特征指數(shù)，決定分布尾的衰減速率，當(dāng)α=1時，函數(shù)服從柯西分布，當(dāng)α=2時函數(shù)服從高斯分布；β為偏斜參數(shù)，當(dāng)β=0時圖形左右對稱，當(dāng)β<0與β>0時圖形分別向左和向右偏斜；σ為尺度參數(shù)，決定著分布關(guān)于μ的離散程度；μ為位置參數(shù)，當(dāng)μ增大時，圖形整體向左平移，反之向右平移。

文獻[18]利用均勻分布和指數(shù)分布進行變量替代，并證明了Levy分布隨機變量生成的原理以及函數(shù)表達式

(2)

其中：V服從區(qū)間為(-π/2,π/2)的均勻分布；W服從均值為1的指數(shù)分布。

Sα,β和Bα,β的定義表達式為

設(shè)定β=0，σ=1，μ=0，不同特征指數(shù)α下Levy噪聲的時域分布如圖1所示。可以看出，Levy噪聲的沖擊性隨著特征指數(shù)α的減小而增強。當(dāng)α=2時，已無尖峰沖擊特征，其概率密度分布曲線與高斯白噪聲相同。圖2為不同特征指數(shù)α所對應(yīng)的Levy噪聲概率密度分布曲線圖。圖像中的曲線中部越尖，其分布的拖尾性越嚴重，時域分布的沖擊性越顯著。

圖1 不同特征參數(shù)α所對應(yīng)Levy噪聲的時域分布圖(β=0，σ =1，μ=0)Fig.1 Time domain distribution map of Levy noise corresponding to different characteristic parameters α (β=0，σ =1，μ=0)

圖2 不同特征參數(shù)α所對應(yīng)Levy噪聲的概率密度分布曲線(β=0，σ =1，μ=0)Fig.2 Probability density distribution curve of Levy noise corresponding to different characteristic parameters α (β=0，σ =1，μ=0)

2 系統(tǒng)理論與檢測指標(biāo)

2.1 雙穩(wěn)態(tài)隨機共振理論模型

雙穩(wěn)系統(tǒng)是隨機共振系統(tǒng)中最為經(jīng)典的一種非線性系統(tǒng)，雙勢阱結(jié)構(gòu)的朗之萬方程(LE)可以作為描述雙穩(wěn)態(tài)隨機共振系統(tǒng)的典型模型，其方程[19]為

(3)

其中：ε(t)為沖擊信號；η(t)為Levy噪聲，當(dāng)噪聲特征指數(shù)α=2時，η(t)為高斯白噪聲。

勢函數(shù)U(x)的函數(shù)表達式為

(4)

其中：a和b為雙勢阱函數(shù)的系統(tǒng)參數(shù)。

沖擊信號的函數(shù)表達式為

(5)

其中：A為幅值；τ為半峰寬度，τ越小其尖峰性越強；t0為沖擊時刻。

非對稱雙穩(wěn)態(tài)隨機共振系統(tǒng)所對應(yīng)的郎之萬方程式為

(6)

其中：C為非對稱因子。

通過調(diào)節(jié)C的值可以實現(xiàn)左右諧振腔不對稱的調(diào)控。設(shè)定C=-0.1，其余參數(shù)與圖3(a)保持相同，非對稱雙勢阱函數(shù)曲線如圖3(b)所示。以圖3(b)為系統(tǒng)模型，輸入振幅為0.1，衰減系數(shù)為10的周期沖擊信號，系統(tǒng)的輸入輸出如圖4所示?？梢钥闯觯到y(tǒng)輸出信號的沖擊時刻與系統(tǒng)輸入的相同，對應(yīng)的沖擊信號幅值無明顯差距。

式(6)對應(yīng)的Fokker-Plank方程為

圖3 雙穩(wěn)態(tài)系統(tǒng)勢函數(shù)圖形Fig.3 Potential function graph of bistable system

圖4 非對稱雙穩(wěn)態(tài)系統(tǒng)的沖擊響應(yīng)Fig.4 Impact response of asymmetric bistable system

(7)

其中：A(x)=ax-bx3+C+ε(t)；B(x)=D；D=σα，為噪聲強度。

當(dāng)α為2時，η(t)退化成高斯白噪聲，此時

〈η(t)η(s)〉=2Dδ(t-s)

(8)

(9)

(10)

其中：Nst為歸一化常數(shù)。

當(dāng)α≠2時，按照文獻[20] 的有限差分方法進行數(shù)值計算，得到不同噪聲強度下密度函數(shù)ρ(s,t)隨x的變化趨勢。圖5為噪聲特征指數(shù)α分別為1.5與1.8時的穩(wěn)態(tài)概率密度函數(shù)曲線。其中：非對稱因子C為 -0.2；噪聲強度D為0.5。

圖5 非對稱雙穩(wěn)態(tài)系統(tǒng)穩(wěn)態(tài)概率密度函數(shù)曲線Fig.5 Steady state probability density function curve of asymmetric bistable system

2.2 檢測指標(biāo)的構(gòu)建

2.2.1 峭度指標(biāo)

峭度指標(biāo)[21]為無量綱，是一個可以衡量一組數(shù)據(jù)離散程度的指標(biāo)，常作為沖擊信號的檢測指標(biāo)，其數(shù)學(xué)表達式為

(11)

其中：x為信號序列；n為信號的長度。

2.2.2 互關(guān)聯(lián)系數(shù)

互關(guān)聯(lián)系數(shù)可以作為兩組序列同步性變化的指標(biāo)，該指標(biāo)不受均值和相位的影響，是衡量兩組序列協(xié)調(diào)性的重要標(biāo)準(zhǔn)之一。其表達式[22]為

(12)

根據(jù)互關(guān)聯(lián)性可知，|Inc|≤1，當(dāng)兩組序列成線性關(guān)系時，|Inc|=1。

2.2.3KI的構(gòu)造

對于多沖擊信號的檢測，單以Kr做為衡量指標(biāo)容易造成漏檢且也不利于源信號的特征檢測。針對峭度指標(biāo)在沖擊信號檢測的不足，筆者構(gòu)造了一種新的無量綱檢測標(biāo)準(zhǔn)(KI)，其表達式為

(13)

利用KI最大化作為沖擊信號的檢測指標(biāo)，對沖擊信號進行檢測。以KI作為信號檢測指標(biāo)，既能夠考慮到?jīng)_擊信號的尖峰特性，又能顧及到?jīng)_擊信號的原始特征，彌補了峭度指標(biāo)作為沖擊信號檢測依據(jù)的不足，同時增強了系統(tǒng)的抗干擾性能。

2.2.4 不同檢測指標(biāo)下單沖擊信號的檢測對比

設(shè)定Levy噪聲參數(shù)α，β，σ和μ分別為1.9，0，1和0，噪聲強度D為0.1，單沖擊信號幅值與衰減系數(shù)分別為0.8和3，采樣頻率和采樣點數(shù)分別為5Hz和5 000，雙穩(wěn)態(tài)隨機共振系統(tǒng)參數(shù)a和b的掃描步長分別為0.01和0.1，掃描區(qū)間均為[0,4]。由于Levy含有較強的隨機性和沖擊性，對系統(tǒng)進行連續(xù)50次重復(fù)實驗并求其均值，得到系統(tǒng)輸出如圖6所示。圖6(a)為以峭度指標(biāo)作為信號質(zhì)量檢測依據(jù)時，對稱雙穩(wěn)態(tài)系統(tǒng)輸出的時域圖形。此時系統(tǒng)最佳參數(shù)a和b分別為1.2和2.9，系統(tǒng)輸出特征信號幅值為0.302。圖6(b)為基于KI檢測指標(biāo)對稱雙穩(wěn)態(tài)系統(tǒng)輸出的最佳時域圖。此時a和b分別為0.74和1.9，其中系統(tǒng)輸出特征信號幅值為0.686，是圖6(a)特征信號幅值的2.27倍。由于Levy噪聲含有較強隨機性和沖擊性，僅以峭度指標(biāo)作為信號質(zhì)量檢測依據(jù)，容易受到噪聲干擾而出現(xiàn)檢測失真現(xiàn)象。如圖6(c)所示，原始沖擊信息已完全淹沒在噪聲之中，以KI作為信號質(zhì)量檢測指標(biāo)未出現(xiàn)此類現(xiàn)象。

圖6 不同檢測指標(biāo)下對稱雙穩(wěn)態(tài)系統(tǒng)的輸出Fig.6 Output of symmetric bistable system under different detection indexes

3 Levy噪聲背景下的沖擊信號檢測

3.1 系統(tǒng)輸出KI值隨系統(tǒng)參數(shù)a的變化趨勢

設(shè)定沖擊信號幅值與衰減系數(shù)分別為0.8和3，噪聲參數(shù)α，β，σ和μ分別為1.5，0，1和0，噪聲強度D為0.1，系統(tǒng)參數(shù)b為0.1，采樣頻率和采樣點數(shù)與圖6相同，非對稱因子C為0，此時系統(tǒng)為對稱雙穩(wěn)系統(tǒng)，系統(tǒng)參數(shù)a的掃描區(qū)間為(0.5,8]，步長為0.05。系統(tǒng)輸出KI值隨系統(tǒng)參數(shù)a的變化趨勢如圖7所示?？梢钥闯觯到y(tǒng)輸出KI值隨系統(tǒng)參數(shù)a呈現(xiàn)先增大后減小的趨勢。當(dāng)a為1.1時，KI達到最大值，KI值為0.967。設(shè)定系統(tǒng)參數(shù)a為1.05，其余參數(shù)保持不變，系統(tǒng)輸入與輸出如圖8所示。圖8(a)和(c)中沖擊信號已完全被Levy噪聲淹沒，無法發(fā)現(xiàn)沖擊信號的任何特征信息。從圖8(b)中可以看到所含沖擊信號的信息。圖8(d)中包含了與沖擊信號頻譜相似的特征信息，說明弱沖擊信號在強Levy噪聲背景下已被有效地識別出來。

圖7 系統(tǒng)輸出KI值隨系統(tǒng)參數(shù)a的變化趨勢Fig.7 The change trend of the output KI value of the system with the system parameter a

圖8 系統(tǒng)參數(shù)a為1.05時系統(tǒng)的輸入與輸出Fig.8 System input and output when the system parameter a is 1.05

3.2 系統(tǒng)輸出KI值隨非對稱因子C的變化趨勢

設(shè)定系統(tǒng)參數(shù)a為1.1，非對稱因子C的掃描區(qū)間為[-1,1]，步長為0.01，其余參數(shù)與圖7保持一致，圖9為系統(tǒng)輸出KI值隨非對稱因子C的變化趨勢。為預(yù)防Levy噪聲給系統(tǒng)帶來的微小波動，此數(shù)據(jù)為重復(fù)20次仿真實驗求得的數(shù)據(jù)均值。從圖9 可以看出，KI值隨非對稱因子C的變化趨勢為先增大后減小。當(dāng)C為負數(shù)時，如圖3(b)所示，右勢阱低于左勢阱，整體向右傾斜，此時含噪聲信號與系統(tǒng)沒有達到最佳協(xié)同狀態(tài)，當(dāng)非對稱因子C逐步增大時，右勢阱升高，左勢阱相對降低，此時KI值也逐步增大。當(dāng)C大于0時，右勢阱高于左勢阱。當(dāng)C為0.54時，KI值為1.04并達到最大值，含噪聲的信號與系統(tǒng)達到最佳隨機共振匹配狀態(tài)。隨著C繼續(xù)增加，右勢阱相對左勢阱繼續(xù)增高，打破了最佳協(xié)同狀態(tài)，因此KI值也隨之迅速減小。當(dāng)C增大到0.84左右時，由于此時系統(tǒng)模型已不再適用而產(chǎn)生失真現(xiàn)象(此時輸出為復(fù)數(shù))。

設(shè)定非對稱因子C為0.54，其余參數(shù)與圖9保持相同，系統(tǒng)輸出如圖10所示。對比于圖8的C為0，此時系統(tǒng)輸出KI值增長了7.03%，從系統(tǒng)輸出效果上來看，圖10(a)的相對沖擊幅值明顯高于圖8(b)。頻域上，相比于圖8(d)來講，圖10(b)中系統(tǒng)輸出的頻譜較高，且高頻分量衰減的更低，說明無用噪聲能量更多地被轉(zhuǎn)化到了信號之中?？梢姡菍ΨQ系統(tǒng)輸出效果優(yōu)于對稱雙穩(wěn)系統(tǒng)。

圖9 系統(tǒng)輸出KI值隨非對稱因子C的變化趨勢Fig.9 The change trend of the KI value of the system with the asymmetric factor C

圖10 系統(tǒng)最佳參數(shù)時系統(tǒng)輸出的時域與頻域圖Fig.10 Time domain and frequency domain diagram of system output with optimal system parameters

3.3 系統(tǒng)輸出KI值與噪聲強度D的變化趨勢

設(shè)定a為1.1，C為0.54，噪聲強度D的掃描區(qū)間為(0,0.5)，步長為0.01，其余參數(shù)與圖9保持一致。由于Levy噪聲的不穩(wěn)定性，使得系統(tǒng)每次輸出都會發(fā)生微小的波動，為了保證實驗數(shù)據(jù)的合理性，重復(fù)30次實驗并求其均值，得到系統(tǒng)輸出KI值隨噪聲強度D的變化趨勢如圖11所示?？梢钥闯?，KI值隨噪聲強度D增大整體的變化為下降趨勢，但也會出現(xiàn)輕微的浮動，如D為0.19處的KI值明顯高于D為0.18，0.17與0.16處，這時因為此時系統(tǒng)的協(xié)同效應(yīng)相對更加匹配。噪聲強度D大于0.4時，KI值基本趨于0.5以下。

圖11 系統(tǒng)輸出KI值隨噪聲強度D的變化趨勢Fig.11 The variation trend of the output KI value of the system with the noise intensity D

4 工程實例

圖12 外圈故障信號原始數(shù)據(jù)Fig.12 Original data of outer ring fault signal

圖13 并聯(lián)非對稱雙穩(wěn)態(tài)系統(tǒng)設(shè)計流程圖Fig.13 Design flow chart of parallel asymmetric bistable system

圖14 并聯(lián)非對稱雙穩(wěn)系統(tǒng)的輸出Fig.14 Output of parallel asymmetric bistable system

(14)

表1 并聯(lián)非對稱雙穩(wěn)態(tài)各級系統(tǒng)參數(shù)

5 結(jié) 論

非對稱雙穩(wěn)系統(tǒng)可作為沖擊信號檢測的有效系統(tǒng)，性能明顯優(yōu)于對稱雙穩(wěn)系統(tǒng)。KI可作為沖擊信號檢測的有效指標(biāo)，以KI最大化作為信號質(zhì)量改善標(biāo)準(zhǔn)，可將強噪聲背景下的沖擊信號信息檢測出來。利用非對稱雙穩(wěn)態(tài)隨機共振系統(tǒng)與KI檢測指標(biāo)可將Levy噪聲背景下的弱沖擊信號檢測出。工程實際應(yīng)用驗證了該方法的實用性和可行性。