采用粒子群算法的自適應變步長隨機共振研究

張仲海，王 多，王太勇，林錦州，蔣永翔

(1.天津大學 機構理論與裝備設計教育部重點實驗室，天津 300072;2.航天晨光股份有限公司研究院，南京 211100)

自從Benzi等［1］研究古氣象冰川問題提出隨機共振(Stochastic Resonance，SR)概念以來，SR現象受到了廣泛的關注。隨機共振現象是一種非線性現象，它在一定條件下，將部分高頻噪聲能量轉移到低頻信號上，在降低噪聲的同時能夠使淹沒于噪聲中的弱信號得到共振加強，極大地提高輸出信噪比，從而實現從強烈噪聲干擾中檢測微弱信號的目的。傳統的隨機共振受到絕熱近似理論［2］的限制，只適用于小參數信號(信號幅值、信號頻率、噪聲強度遠小于1)，這極大地制約了隨機共振在工程實際中的應用。

為了突破隨機共振只能適應于小參數的限制，許多學者進行了有益探索，取得不少階段性研究成果［3－6］，為隨機共振應用于工程實際提供了可能。例如，文獻［7］提出了二次采樣隨機共振(TSSR)方法，即按照一定的變換尺度將較高采樣信號變換為較低的二次采樣信號，以滿足小參數條件，隨機共振輸出后，再按變換尺度恢復至較高采樣信號。文獻［8］提出一種移頻變尺度隨機共振(FRSR)方法，即通過移頻和尺度變換等手段壓縮分析頻率，從而使之滿足絕熱近似小參數條件。文獻［9］則通過雙穩系統參數的調節，從而達到共振狀態。以上方法為隨機共振應用于工程實測信號提供了理論基礎，但在實際應用中，如何實現參數的自適應選取是一個難題。針對自適應隨機共振的問題，文獻［10］提出了一種基于近似熵測度的自適應隨機共振方法，在固定雙穩系統結構參數b的條件下，對變步長(或二次采樣)隨機共振［8，10］系統參數a和計算步長h進行優化，當參數b選取不合適時，則無法達到真正意義上最優的共振狀態。文獻［11］提出了一種基于遺傳算法的自適應隨機共振算法，能對雙穩系統的結構參數a、b進行優化，但在處理大參數信號時需要在優化前手動進行移頻和尺度變換等操作，當調制頻率fc和變尺度壓縮率R選取不恰當時會影響隨機共振的輸出效果。對于某一確定的含噪信號，需要對隨機共振系統進行多參數的同步調節，當且僅當雙穩系統參數、信號頻率與噪聲相互之間達到最佳協同時，才能達到最佳的共振狀態。

變步長隨機共振能突破絕熱近似理論對于小參數的限制，而應用于大參數條件下的工程實測信號，但如何對結構參數a，b和計算步長h進行自適應選取仍是一個難題。粒子群算法［12－13］作為一種多變量全局優化方法，沒有遺傳算法的交叉和變異，因此效率高、速度快，具有很強的工程應用價值。本文以變步長雙穩隨機共振系統為研究對象，利用粒子群算法的全局并行搜索優化能力，以雙穩系統輸出的信噪比作為粒子群算法的適應度函數，對變步長雙穩隨機共振系統的結構參數a，b和計算步長h進行同步優化，使a，b和h等三個參數達到最佳的協同，從而最優地實現對大參數條件下微弱信號的自適應檢測。

1 變步長隨機共振的基本原理

1.1 雙穩隨機共振概述

在隨機共振的研究中，非線性雙穩系統通常由勢函數表示:

其中，a和b是雙穩系統的結構參數。

在不考慮噪聲的情況下得到一維非線性雙穩系統的確定性動力學系統方程:

由定態方程dx/dt=0得到雙穩系統的三個解為一個不穩定的定態解x1=0和兩個穩定的定態解x2，3=從物理意義上理解，雙穩系統的響應輸出可以表現為一個粒子在雙穩勢阱內的運動，如圖1所示。

由圖1可以看出，系統參數a和b不僅調節雙穩系統勢壘ΔU=a2/4b，而且還改變著雙穩系統兩個勢阱之間的距離

圖1 非線性雙穩系統勢函數Fig.1 Non-linear bi-stable system potential function

當以單周期正弦信號Asin(2πf0t)和白噪聲n(t)為輸入信號時，對應的Langevin方程為:

這里，

其中D是噪聲強度。

從粒子動力學角度來考察，式(3)描述了一個過阻尼的質點布朗運動。當A和n(t)均為0，即沒有調制和噪聲作用時，根據系統不同的初始狀態，質點將處于雙穩勢阱中的某一個。當A＞0且足夠小時，在無噪情況下，由于外加信號的驅動作用，整個系統不再處于平衡狀態，勢阱按信號頻率發生周期性傾斜，但是粒子的運動將被限制在某一勢阱內。然而，在引入噪聲的情況下，即使A＜Ac，甚至A?Ac時，質點將越過勢壘而從當前勢阱躍遷到另一勢阱。如果外加噪聲的強度合適，這種大幅度躍遷與周期驅動力達到很好的協同，那么原來兩個勢阱之間的隨機躍遷運動就會變成與周期調制信號頻率相一致的有序躍遷運動，隨機共振現象就發生了。

1.2 變步長隨機共振及其存在的問題

式(3)是一種非線性隨機微分方程，可通過四階Runge-Kutta法進行數值求解，具體算法如下:

其中，n=1，2，…，N。

式(4)中，Sn和xn分別是雙穩系統輸入S(t)=Asin(2πf0t)+n(t)和輸出X(t)的第n個采樣值，h=1/fs(fs)為數值計算步長。

當待測信號為小參數信號，滿足絕熱近似理論時，無需改變計算步長僅通過調節雙穩系統結構參數就達到很好的共振狀態。而當待測信號為大參數信號時，由于絕熱近似理論的限制，僅通過調節系統參數無法產生共振。為了克服絕熱近似理論對小參數的限制，冷永剛等［7］提出了二次采樣隨機共振思想，即對于大參數的待測信號，通過對其進行二次采樣，將其轉化為小參數，經過雙穩系統處理得到共振輸出后，按相應的變換尺度恢復還原到實際信號頻率。設有一大參數的含噪信號，信號頻率為f，采樣頻率為fs，二次采樣頻率為fsr，通過改變式(4)中的計算步長h，使其等于二次采樣頻率的倒數，即h=1/fsr，把原信號頻率f變換成f0=f·fsr/fs，通過式(4)求解，就可得出頻率f0下的共振輸出。在實際數值求解過程中，可先調節計算步長h使雙穩系統達到共振狀態，得到變換后的信號頻率f0，再按變換尺度R=hfs還原恢復實際信號頻率f=Rf0。以上便是變步長隨機共振(SCSR)思想。

大量分析表明，對大參數信號進行變步長隨機共振處理時，通過適當的改變計算步長h，能拓寬隨機共振系統輸出的頻帶寬度。不過僅改變計算步長h不能獲得很好的輸出，還需要結合非線性系統結構參數a和b的調節，以改善非線性系統特性，進而得到待檢測的微弱信號。但是，如何對結構參數a，b和計算步長h進行聯合調節，進而使a，b，h三個參數達到最佳協同，獲得雙穩系統的最優輸出，目前還沒有一個定量的規律可循，只能憑經驗探索性地選取，給隨機共振在工程實際中的應用帶來很大的不便。

2 基于粒子群算法的自適應變步長隨機共振

如何對變步長隨機共振系統的結構參數a，b和計算步長h進行自適應選取是一個難題，本文利用粒子群算法來實現結構參數 a，b和計算步長 h的同步優化。

2.1 粒子群優化算法

Kennedy等［12－13］研究鳥群覓食行為時，受到啟發而提出了粒子群優化算法(Particle Swarm Optimization，PSO)。該算法實現方便，與遺傳算法相比需要設置的參數少，是一種高效、實用的搜索優化算法。

圖2 粒子群算法優化搜索示意圖Fig.2 Searching schematic diagram of PSO

在粒子群優化算法中，每個粒子表征某一確定的待優化問題函數的一個可能解。首先，初始化一組粒子種群，粒子的初始速度和位置隨機產生，之后，粒子群追隨當前最優粒子，不斷更新自己的速度和位置，在多維解空間中搜索，經過若干次迭代找出最優解。在每次迭代中，粒子通過跟蹤其自身的當前最優解(即個體極值Pbest)和整個粒子群的當前最優解(即全局極值Nbest)來更新自己的速度和位置，直到達到設定的最大迭代次數或找到最優解。粒子的優化運動軌跡如圖2所示。

PSO算法通常的數學描述為:設在某D維空間中，種群 X=(x1，…，xi，…，xm)由 m 個粒子組成。xi(t)=(xi1，xi2，…，xiD)T和 vi(t)=(vi1，vi2，…，viD)T分別表示t時刻種群中第i個粒子的位置和速度;Pbesti(t)=(Pbesti1，Pbesti2，…，PbestiD)T表示 t時刻第 i個粒子的個體極值，t時刻整個種群的全局極值則由Nbest(t)=(Nbest1，Nbest2，…，NbestD)T表示。那么，粒子 xi將按式(5)來更新其速度和位置。

式中，j=1，2，…，D;t為當前進化代數;r1，r2為均勻分布于［0，1］之間的隨機數;c1，c2為學習因子，通常取c1=c2=2;w為慣性權重，即保持原來速度的系數。

慣性權重的計算公式如下:

式中;wmax為慣性權重上限;wmin為慣性權重下限;t為當前進化代數;Tmax為最大進化代數。

粒子速度和位置更新方式如圖3所示。

圖3 粒子速度和位置更新方式示意圖Fig.3 Concept of modification of a searching point by PSO

圖3中，vPbesti表示粒子在解空間中朝著個體極值進行搜索的速度，對應于c1r1(t)［Pbesti(t)－xi(t)］速度部分;vNbesti表示粒子在解空間中朝著全局極值進行搜索的速度，對應于c2r2(t)［Nbest(t)－xi(t)］速度部分。

2.2 目標優化函數

粒子群算法是根據各個粒子的適應度大小來調整進化搜索能力的。而適應度函數和目標優化函數是相關的。此處采用的適應度函數即為目標優化函數——隨機共振輸出信噪比。本文的適應度函數為:

式中:sr(a，b，h)為變步長隨機共振的輸出結果;SNRout(sr(a，b，h))表示隨機共振輸出的信噪比。

已知隨機共振輸出信噪比定義［14］如下:

式中:F0為信號頻率;S(F0)為信號功率;P為系統總功率，包括信號功率和噪聲功率;P－S(F0)即為噪聲功率。

設輸入信號為Asin(2πF0t)+n(t)，該含噪信號經采樣頻率為Fs的采樣得到長度為L的離散序列Zl。Zl經過二次采樣頻率為Fsr的變步長隨機共振，輸出信號sr(a，b，h)中頻率分量F'0=F0Fsr/Fs對應于輸入頻率F0。設F'0分量的單邊譜幅值為 X(k0)，且有 k0=LF'0/Fsr=LF0/Fs，由式(8)可得:

2.3 基于粒子群算法的自適應變步長隨機共振流程

粒子群優化算法容易實現且具有全局并行搜索優化能力，本文利用粒子群優化算法以2.2節所述的隨機共振輸出信噪比為適應度函數，對變步長隨機共振系統的結構參數a、b和計算步長h等三個參數進行同步優化。算法的流程如圖4所示。

下面給出其實現的具體步驟:

步驟1:種群初始化。設置種群數量、參數a、b和h的搜索范圍以及最大進化代數Tmax，最大搜索速度取最大調整步長的10%～20%，這里的最大調整步長指的是所設定的粒子位置范圍的上限值減去粒子位置范圍的下限值所得的差值。隨機初始化搜索點的位置xi(0)和速度vi(0)，設置各個粒子的Pbesti坐標為其當前位置xi(0)，并計算出其相應的個體極值，記錄整個粒子群中個體極值最大的粒子序號，設置Nbest為該最大粒子的當前位置。

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

步驟3:粒子的更新。根據式(5)更新所有粒子的速度和位置。

圖4 粒子群優化算法流程圖Fig.4 Flowchart of PSO

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

步驟5:檢測結果。根據對a、b和h優化輸出的最優解，對原始信號進行變步長隨機共振，尺度恢復后得到最終的微弱信號檢測結果。

3 仿真數據分析

基于粒子群算法的自適應變步長隨機共振通過自動聯合調節結構參數a，b和步長h，得到了自適應條件下雙穩系統的最優輸出。下面利用仿真數據證明該自適應方法的有效性。

設輸入信號為u(t)=A0sin(2πf0t)+n(t)，其中，A0=0.2，f0=20 Hz，添加均值為0、方差 D=3.1 的高斯白噪聲 n(t)，采樣頻率 fs=2 048 Hz，采樣點數 n=2 048。該輸入信號的原始時域波形及其頻譜分別如圖5(a)和圖5(b)所示，由于強噪聲的加入，從圖5(a)上很難發現周期成分，在圖5(b)上也很難辨認出60 Hz的頻率分量。

圖5 大參數數據的仿真結果Fig.5 Simulation results of large parameters data

現利用本文提出的自適應變步長隨機共振方法對該信號進行處理。初始化種群:設置種群數量為40，a、b、h的搜索范圍分別為［0.1，10］、［0.1，1 000］和［0.02，0.2］，最大搜索速度取最大調整步長的20%，即a、b、h 的最大搜索速度分別為 1.98、199.98 和 0.036，最大進化代數為150，最小錯誤標準為1×10－4。從圖5c的收斂曲線可以看出，經過62次迭代，算法收斂，輸出的最優參數分別為 a=6.33、b=217.39、h=0.041 3。將最優參數代入變步長隨機共振系統，對原始信號進行隨機共振處理后分別得到圖5(d)和圖5(e)所示的時域波形和頻譜。從圖5d可以看出，噪聲已經被極大地削弱了，在圖5(e)上可以發現頻率為0.236 7 Hz的分量非常突出。按變換尺度 R=hfs=0.041 3×2 048=84.582 4還原恢復后，可以得到f0=Rf=84.582 4×0.236 7=20 Hz，其頻率正好對應于原始信號中的20 Hz頻率成分。仿真結果表明，將粒子群優化算法用于變步長隨機共振系統能實現對參數a、b和h的自適應最優選取，從而快速有效的檢測出大參數條件下的微弱信號。

4 工程應用

實驗選用型號為6205－2RS JEM SKF的深溝球軸承，該軸承的尺寸和故障頻率如表1和表2如示。

表1 滾動軸承6205－2RS的尺寸參數(厘米)Tab.1 Size parameters of 6205 －2RS

表2 滾動軸承6205－2RS的故障頻率(轉頻的倍數)Tab.2 Defect frequencies of 6205 －2RS

使用電火花加工技術在該軸承內圈上布置了單點故障，故障直徑為0.017 78 cm，根據表2可以算出該滾動軸承內圈故障的特征頻率為156.14 Hz。試驗中，該軸承用于支承電機軸，電機轉速為1 730 r/min，使用加速度傳感器采集振動信號，采樣頻率為Fs=12 kHz，采樣點數n=8 192。圖6(a)和圖6(b)所示的分別為原始采樣信號的時域波形和幅值譜，在圖6(b)的頻譜圖上根本無法辨認156.14 Hz的滾動軸承內圈故障頻率。由于原始信號中沖擊成分較為明顯，對其進行包絡解調處理，所得的時域波形和幅值譜分別如圖6(c)和圖6(d)所示，從圖6(d)的包絡譜中可以看到156 Hz的故障頻率，但譜線很不明顯。

現對該滾動軸承故障信號的包絡進行基于粒子群算法的自適應變步長隨機共振處理。首先初始化粒子種群:設置種群數量為80，a、b和h的搜索范圍分別為［0.01，30］、［0.01，15 000］和［0.002，0.8］，最大搜索速度取最大調整步長的20%，即a、b、h的最大搜索速度分別為 5.998、2 999.998 和 0.159 6，最大進化代數為200，最小錯誤標準為1×10－4。

從圖7的收斂曲線可以看出，經過103次迭代，算法收斂，輸出的最優參數分別為 a=8.15、b=2 220.6、h=0.048 7。將最優參數代入變步長隨機共振系統，對包絡信號進行隨機共振處理后分別得到圖6(e)和圖6(f)所示的時域波形和頻譜。對比圖6(c)和圖6(e)后可以發現，隨機共振處理后的時域波形中的噪聲成分被極大的削弱了。在圖6(f)中可以非常清楚的看到0.266 8 Hz的頻率分量及其二倍頻，按變換尺度R=hFs=0.048 5×12 000=584.4 還原恢復后，可以得到F0=RF=584.4 ×0.266 8≈155.9 Hz，考慮到存在計算舍入誤差，即該頻率正好是軸承內圈的故障特征頻率，這與滾動軸承存在內圈故障的事實相吻合。

圖6 滾動軸承內圈故障的診斷結果Fig.6 Diagnosis results of the fault of inner ring

圖7 PSO算法的最優收斂曲線Fig.7 The optimal convergence curve of PSO algorithm

5 結論

本文提出了一種采用粒子群優化算法的自適應變步長隨機共振方法，選用共振輸出的信噪比作為粒子群算法的適應度函數，利用粒子群算法的全局并行搜索優化能力，對變步長隨機共振系統的多個參數進行同步優化，實現了結構參數a、b和計算步長h等三個參數的自適應選取，從而最優檢測出強噪聲背景下的高頻微弱信號。該方法克服了單參數優化及步長選取依賴經驗的缺點，充分體現了聯合調參的思想。將該方法用于仿真微弱信號的檢測及滾動軸承的故障診斷，結果表明所提方法能快速有效地檢測出大參數條件下的微弱信號。

［1］Benzi R，Parisi G，VulpianiA. Theoryofstochastic resonance in climatic chaner［J］.SIAM Journal on Applied Mathematics，1983，43(3):565 －578.

［2］Mcnamara B， Wiesenfeld K， RoyR. Observation of Stochastic Resonance in a Ring Laser［J］.Physical Review Letters.1988，60(25):2626 －2629.

［3］夏均忠，劉遠宏，馬宗坡，等.基于調制隨機共振的微弱信號檢測研究［J］.振動與沖擊，2012，31(3):132 －135，140.XIA Jun-zhong，LIU Yuan-hong，MA Zong-po，et al.Weak signal detction based on the modulated stochastic resonance［J］.Journal of Vibration and Shock.2012，31(3):132 －135，140.

［4］李曉龍，冷永剛，范勝波，等.基于非均勻周期采樣的隨機共振研究［J］.振動與沖擊，2011，30(12):78 －84.LI Xiao-long， LENG Yong-gong，FAN Sheng-bo，et al.Stochastic resonance based on periodic non-uniform sampling［J］.Journal of Vibration and Shock.2011，30(12):78 －84.

［5］郝 研，王太勇，萬 劍，等.基于級聯雙穩隨機共振和多重分形的機械故障診斷方法研究［J］.振動與沖擊，2012，31(8):181－185.HAO Yan，WANG Tai-yong，WAN Jian，et al.Mechanical faultdiagnosis based on cascaded bistable stochastic resonance and multi-fractal［J］.Journal of Vibration and Shock，2012，31(8):181 －185.

［6］趙艷菊，王太勇，徐 躍，等.雙穩隨機共振降噪下的經驗模式分解研究［J］.振動與沖擊，2009，28(3):149－151.ZHAO Yan-ju，WANG Tai-yong，XU Yue，et al.Empirical mode decomposition based on bistable stochastic resonance denoising［J］.Journal of Vibration and Shock，2009 ，28(3):149－151.

［7］冷永剛，王太勇.二次采樣用于隨機共振從強噪聲中提取弱信號的數值研究［J］.物理學報，2003，52(10):2432－2437.LENG Yong-gang，WANG Tai-yong.Numerical research of twice sampling stochastic resonance for the detection of a weak signal submerged in a heavy Noise［J］.Acta Physica Sinica，2003，52(10):2432 －2437.

［8］Tan J，Chen X，Wang J，et al.Study of frequency-shifted and re-scaling stochastic resonance and its application to fault diagnosis［J］.Mechanical Systems and Signal Processing，2009，23(3):811－822.

［9］陳 敏，胡蔦慶，秦國軍，等.參數調節隨機共振在機械系統早期故障檢測中的應用［J］.機械工程學報，2009，45(4):131－135.CHEN Min，HU Niao-qing，QIN Guo-jun，et al.Application of parameter-tuning stochastic resonance for detecting early mechanical faults［J］.Journal of Mechanical Engineering，2009，45(4):131－135.

［10］Qiang L，Taiyong W，Yonggang L，et al.Engineering signal processing based on adaptive step－changed stochastic resonance［J］.Mechanical Systems and Signal Processing，2007，21(5):2267－2279.

［11］Wang J，Zhang Q，Xu G H.Genetic Stochastic Resonance:A new fault diagnosis method to detect weak signals in mechanical systems［J］.Advanced Science Letters，2011，4(6－7):2508－2512.

［12］Kennedy J，Eberhart R.Particle swarm optimization［C］.Perth，Aust:IEEE，1995.

［13］Eberhart R，Kennedy J.A new optimizer using particle swarm theory［C］.New York，NY，USA:IEEE，1995.

［14］萬 頻，詹宜巨，李學聰，等.一種單穩隨機共振系統信噪比增益的數值研究［J］.物理學報，2011，60(4):60－66.WAN Pin，ZHAN Yi-ju，LI Xue-cong，et al.Numerical research ofsignal-to-noise ratio gain on a monostable stochastic resonance［J］.Acta Physica Sinica，2011，60(4):60－66.