汪洋百慧， 鄭永軍， 羅 哉

(中國計量大學(xué) 計量測試工程學(xué)院，浙江 杭州 310018)

1 引 言

在微弱信號檢測中，有用信號總淹沒于強噪聲中，這阻止了有用信號的提取，因此人們總想濾除噪聲[1,2]。但經(jīng)過大量研究表明噪聲在某些條件下反而可以產(chǎn)生積極有利的作用[3～5]。其中隨機共振在微弱周期信號、噪聲和非線性系統(tǒng)的協(xié)同作用下增強輸出信號的幅值，以此達(dá)到微弱信號檢測的目的[6,7]。目前為止，大部分研究人員仍將重點放在整數(shù)階模型上。若系統(tǒng)更加復(fù)雜，仍使用傳統(tǒng)的整數(shù)階模型將會存在很大的局限性，而分?jǐn)?shù)階模型能很好地解決以上問題[8～10]。

隨機共振通常發(fā)生在布朗粒子反常擴散過程中，具有非馬爾可夫特性。由于分?jǐn)?shù)階微積分特點是在時間上的記憶性和空間上的關(guān)聯(lián)性[11～19]，因此可用來描述這類現(xiàn)象，并且基于分?jǐn)?shù)階系統(tǒng)的信噪比等衡量指標(biāo)更優(yōu)異。最常用的分?jǐn)?shù)階微積分定義是Grünwald-Letnikov(G-L)定義，Riemann-Liouville(R-L)定義和Caputo定義。然而有學(xué)者指出上述分?jǐn)?shù)階微積分的內(nèi)核是非局部的，但具有奇異性[20]。Atangana和Baleanu在2016年給出了一種新形式的分?jǐn)?shù)階微積分，有效解決了這一問題[21]。其中借助廣義Mittag-Leffler(M-L)函數(shù)作為非局部和非奇異核，使其具有均方位移的交叉性質(zhì)并保證了信息在延展的開始和結(jié)束時的一致性[22～24]。

本文根據(jù)Atangana和Baleanu新提出的分?jǐn)?shù)階微積分，構(gòu)建了雙穩(wěn)系統(tǒng)的分?jǐn)?shù)階Langevin方程，該系統(tǒng)解決了奇異性問題，能夠更好地描述隨機共振的動力學(xué)特性。同時基于該方程，在Simulink環(huán)境下對其進(jìn)行仿真實驗，采用控制單一變量法，分別在有無噪聲時調(diào)節(jié)分?jǐn)?shù)階求導(dǎo)階次、固定分?jǐn)?shù)階求導(dǎo)階次不變調(diào)節(jié)噪聲強度大小這3種情況下，分析隨機共振的動態(tài)特性。

2 Atangana-Baleanu分?jǐn)?shù)階微積分

2.1 Atangana-Baleanu(簡稱A-B)分?jǐn)?shù)階微積分的定義和性質(zhì)

當(dāng)f∈H1(c,b)，b>c，α∈[0,1]時，基于Caputo的A-B分?jǐn)?shù)階導(dǎo)數(shù)的定義如式(1)所示：

(1)

M(α)的定義如式(2)所示[25]：

(2)

式(1)中Eα是在整個復(fù)平面上由以下冪級數(shù)定義的一個Mittag-Leffler(M-L)函數(shù)，其定義如式(3)所示：

(3)

當(dāng)α=1時，式(3)可以寫為式(4)所示：

(4)

基于Caputo的A-B分?jǐn)?shù)階導(dǎo)數(shù)更加嚴(yán)謹(jǐn)。但當(dāng)α=0時，原函數(shù)不存在。為了避免這個問題，提出了以下的定義。

當(dāng)f∈H1(c,b)，b>c，α∈[0,1]時，基于R-L的A-B分?jǐn)?shù)階導(dǎo)數(shù)的定義如式(5)所示：

(5)

上述定義將有助于討論現(xiàn)實世界中的問題，并且在使用拉普拉斯變換解決初始條件下的某些物理問題時也將具有很大優(yōu)勢。經(jīng)過簡單計算后得到基于Caputo的A-B分?jǐn)?shù)階導(dǎo)數(shù)的Laplace變換如式(6)所示：

(6)

基于R-L的A-B分?jǐn)?shù)階導(dǎo)數(shù)的Laplace變換如式(7)所示：

(7)

基于上述定義，分?jǐn)?shù)階控制系統(tǒng)可以由式(8)表示：

(8)

式(8)可以寫為式(9)所示的傳遞函數(shù)的形式：

(9)

2.2 分?jǐn)?shù)階導(dǎo)數(shù)的逼近算法

在實際應(yīng)用或仿真實驗中需要對分?jǐn)?shù)階傳遞函數(shù)進(jìn)行求解，由于其形式比較復(fù)雜，直接求解較為困難。本文對傳統(tǒng)的Oustaloup算法進(jìn)行改進(jìn)，用其無限逼近于分?jǐn)?shù)階傳遞函數(shù)，這對于后續(xù)求解分?jǐn)?shù)階控制系統(tǒng)有著較好的提升作用。

改進(jìn)后的Oustaloup算法基于如式(10)所示的傳遞函數(shù)：

H(s)=sα

(10)

在給定的頻段[ωb,ωh]內(nèi)，式(10)可以寫為式(11)所示[26]：

(11)

式中： 0<α<1，s=jω，b>0，d>0。

在頻率ωb<ω<ωh范圍內(nèi)，利用泰勒公式對式(11)展開得到如式(12)所示：

(12)

式中：p(s)的表達(dá)式如式(13)所示：

(13)

忽略式(13)中二階及其無窮小量，sα可以寫為式(14)所示：

(14)

式中：無理小數(shù)部分可以用連續(xù)時間有理模型近似表示為式(15)所示：

(15)

根據(jù)實零點和極點的遞歸分布，第k個零點和極點可以寫為式(16)所示：

(16)

綜上，sα的整數(shù)階近似公式為式(17)所示：

(17)

式中：K為系統(tǒng)增益，如式(18)所示：

(18)

經(jīng)過大量實驗驗證，一般取b=10，d=9。

3 分?jǐn)?shù)階雙穩(wěn)態(tài)隨機共振系統(tǒng)

3.1 隨機共振的基本原理

圖1是一個未加外力的雙穩(wěn)系統(tǒng)在過阻尼情況下的勢函數(shù)，表達(dá)式如式(19)所示：

圖1 雙穩(wěn)系統(tǒng)的勢函數(shù)Fig.1 Potential function of bistable system

(19)

式中：a>0；b>0；x為系統(tǒng)輸入；V(x)為系統(tǒng)輸出；勢阱點為±xm=±(a/b)1/2；勢壘點x=0；勢壘高度為ΔV=a2/4b。

由圖2(a)所示，在未受外力的作用下，雙穩(wěn)系統(tǒng)的兩個勢阱相同且關(guān)于縱坐標(biāo)對稱。

將噪聲加至系統(tǒng)后，布朗粒子會在兩勢阱之間來回躍遷，其躍遷速度可用克萊莫斯(Kramers)速率表示，如式(20)所示[27]：

(20)

式中：D為噪聲強度；ΔV為勢壘高度。

如圖2(a)到2(b)的變化可以看出，若系統(tǒng)存在微弱周期信號作用時，其勢阱不再對稱。此時勢壘高度改變了，從而Kramers速率也發(fā)生了變化。當(dāng)Kramers速率為外加力頻率2倍時，系統(tǒng)能夠產(chǎn)生隨機共振現(xiàn)象[28]。再將微弱信號加至系統(tǒng)，勢函數(shù)將會如圖2所示的變化規(guī)律不斷轉(zhuǎn)變。

圖2 隨機共振現(xiàn)象原理圖Fig.2 Schematic diagram of stochastic resonance phenomenon

上述雙穩(wěn)態(tài)系統(tǒng)的隨機共振現(xiàn)象，其動力學(xué)特性可以用Langevin函數(shù)表示，如式(21)所示[29]:

(21)

式中：A0cos(ωt+φ)為周期性調(diào)制信號；ξ(t)為噪聲信號，通常為白噪聲信號。

3.2 分?jǐn)?shù)階Langevin方程

大量研究表明，在相似條件下，基于分?jǐn)?shù)階微分方程的隨機共振系統(tǒng)具有更加豐富的動力學(xué)現(xiàn)象[30～34]，所以本小節(jié)將整數(shù)階Langevin方程推廣至分?jǐn)?shù)階下。

(22)

(23)

在粘性介質(zhì)中，布朗粒子的運動是一個非馬爾可夫過程，在時間和空間上都具有相關(guān)性。為了解決這個問題，提出了廣義的Langevin方程來描述布朗粒子的反常擴散，如式(24)所示：

(24)

由于阻尼力和噪聲都源于布朗粒子與周圍介質(zhì)分子之間的隨機碰撞，因此隨機力ξ(t)和阻尼核函數(shù)γ(t)滿足以下漲落耗散定理:

〈ξ(t)〉=0,〈ξ(t)ξ(t′)〉=kBTγ(t-t′)

(25)

式中: kB為玻爾茲曼常數(shù)；阻尼核函數(shù)可以寫為：

(26)

將式(26)代入式(24)中：

(27)

又因為式(1)成立的條件是建立在t→τ的時候，即當(dāng)(t-τ)α→0時M-L函數(shù)有：

(28)

所以，式(1)可以寫為：

(29)

由此能夠得到過阻尼情況下的分?jǐn)?shù)階Langevin方程為：

(30)

4 仿真驗證及分析

本章將基于分?jǐn)?shù)階Langevin方程，采用單一變量的方法，研究各個參數(shù)變化對基于新分?jǐn)?shù)階導(dǎo)數(shù)的雙穩(wěn)態(tài)系統(tǒng)隨機共振現(xiàn)象的影響。這里微弱周期信號取F(t)=Acos(2 π ft)，設(shè)系統(tǒng)參數(shù)a=b=1，代入(30)中可得：

(31)

4.1 未加噪聲時改變分?jǐn)?shù)階求導(dǎo)階次

當(dāng)微弱周期信號的幅值A(chǔ)=0.3，頻率為f=0.01 Hz時，無高斯白噪聲輸入的情況下，只改變分?jǐn)?shù)階求導(dǎo)階次。

如圖3(a)所示為基于A-B分?jǐn)?shù)階導(dǎo)數(shù)時，分?jǐn)?shù)階求導(dǎo)階次以0.1的步長從0.1遞增到0.9的輸出響應(yīng)的時域圖?？梢钥闯靓恋淖兓淖兞藙輭靖叨鹊拇笮〔⑹箘葳妩c發(fā)生漂移。當(dāng)α大于0.6時，布朗粒子只能在x=1或x=-1的單一勢阱處運動，說明其沒能越過勢壘。減小α，當(dāng)α小于0.6時，可以看出布朗粒子的運動軌跡發(fā)生了改變，在x=1和x=-1之間運動。說明α值越小，系統(tǒng)的內(nèi)核函數(shù)在時間上的記憶性越強，布朗粒子就有足夠大的能量可以越過勢壘。

圖3 A=0.3時α變化下系統(tǒng)輸出時域圖Fig.3 Time-domain diagram of system output under the change of α when A=0.3

如圖3(b)所示，做了進(jìn)一步分析得到了該系統(tǒng)更加精確的臨界值為0.68。根據(jù)以上結(jié)論可以看出，在未加高斯白噪聲的分?jǐn)?shù)階雙穩(wěn)態(tài)系統(tǒng)中，當(dāng)α大于系統(tǒng)的臨界值時，需要一些額外的能量來驅(qū)動布朗粒子克服勢壘高度；當(dāng)α小于該臨界值時無需外力驅(qū)使即可產(chǎn)生隨機共振。

如圖4所示為信號幅值增加至0.4時，α從0.1以0.1的步長遞增到0.9的輸出響應(yīng)的時域圖。從圖4(a)中可以看出勢壘高度隨著幅值的增加而增大，從圖4(b)中可以看出分?jǐn)?shù)階臨界階次增加至0.92。

圖4 A=0.4時α變化下系統(tǒng)輸出時域圖Fig.4 Time-domain diagram of system output under the change of α when A=0.4

如圖5所示，可以看出系統(tǒng)的臨界階次受信號幅值的影響，信號幅值越大臨界階次越大。這是因為信號幅值的大小會影響勢壘高度的大小，幅值越大勢壘高度越大，則系統(tǒng)的臨界階次也就越大。

圖5 幅值與分?jǐn)?shù)階階次臨界值的關(guān)系圖Fig.5 Relationship between amplitude and critical value of fractional order

4.2 加入噪聲時改變分?jǐn)?shù)階求導(dǎo)階次

由第4.1節(jié)可知當(dāng)α大于系統(tǒng)臨界值時，雙穩(wěn)態(tài)系統(tǒng)沒有產(chǎn)生隨機共振?；诖饲闆r，本節(jié)將討論向雙穩(wěn)態(tài)系統(tǒng)中加入高斯白噪聲，固定噪聲強度D=3.1，外加周期信號幅值A(chǔ)=0.3時是否能產(chǎn)生隨機共振。

如圖6(a)為固定噪聲強度，不同α下系統(tǒng)輸出的時域圖。圖6(b)為其對應(yīng)的頻譜圖?？梢钥闯霎?dāng)α=0.7時，布朗粒子在兩勢阱間做無規(guī)律運動，此時對應(yīng)頻譜值為0.019 96，還未達(dá)到輸入信號頻譜值0.025 03，稱此狀態(tài)為欠共振；當(dāng)α增加至0.85時，布朗粒子逐漸趨于規(guī)律運動，此時的頻譜值0.415 2遠(yuǎn)大于輸入信號頻譜值，雙穩(wěn)態(tài)系統(tǒng)產(chǎn)生隨機共振；α繼續(xù)增加至0.95后，布朗粒子又進(jìn)入無規(guī)則運動，頻譜值較α為0.85時反而小，稱此狀態(tài)為過共振。

圖6 外加噪聲時D不變，α變化下系統(tǒng)輸出的時域圖和頻域圖Fig.6 Noise is added, time-frequency diagram of system output under the change of α when D is unchanged

從圖7中可以看出，雙穩(wěn)態(tài)系統(tǒng)的噪聲強度固定，當(dāng)α小于0.7時，功率譜值均小于0.02，甚至還未達(dá)到原始信號的功率譜；增加α的大小，功率譜值也隨之增大并達(dá)到一個極大值，在此極大值處，雙穩(wěn)態(tài)系統(tǒng)產(chǎn)生了隨機共振；繼續(xù)增加α的大小，系統(tǒng)的勢壘值也在不斷增大，布朗粒子很難再翻越勢阱，此時功率譜值反而減小，抑制了隨機共振的產(chǎn)生。

圖7 功率譜值與分?jǐn)?shù)階求導(dǎo)階次的關(guān)系圖Fig.7 Relationship between the power spectrum value and the fractional derivative order

4.3 固定分?jǐn)?shù)階求導(dǎo)階次時改變噪聲強度

由4.2節(jié)中可知噪聲強度固定，α過小或過大時，雙穩(wěn)態(tài)系統(tǒng)均不會產(chǎn)生隨機共振。基于此情況，本節(jié)將討論當(dāng)α=0.75時，加入噪聲并且改變噪聲強度的隨機共振現(xiàn)象。

如圖8(a)為固定α，不同噪聲強度下系統(tǒng)輸出的時域圖。圖8(b)為其對應(yīng)的頻譜圖。可以看出當(dāng)D=1.5時，布朗粒子在某一勢阱中做無規(guī)律運動，此時頻譜值為0.022 49，還未達(dá)到輸入信號頻譜值0.025 03，此時為欠共振狀態(tài)；當(dāng)D增加至2.2時，布朗粒子逐漸趨于規(guī)律運動，此時的頻譜值0.642 3遠(yuǎn)大于輸入信號頻譜值，雙穩(wěn)態(tài)系統(tǒng)產(chǎn)生隨機共振；D繼續(xù)增加至5后，布朗粒子又進(jìn)入無規(guī)則運動，頻譜值較D為2.2時反而小，此時為過共振狀態(tài)。

圖8 外加噪聲時α不變，D變化下系統(tǒng)輸出的時域圖和頻域圖Fig.8 Noise is added, time-frequency diagram of system output under the change of D when α is unchanged

從圖9中可以看出，分?jǐn)?shù)階求導(dǎo)階次固定，增大噪聲強度，雙穩(wěn)態(tài)系統(tǒng)可以產(chǎn)生隨機共振，但持續(xù)增大噪聲強度反而會抑制隨機共振的產(chǎn)生。當(dāng)噪聲提供的能量遠(yuǎn)大于布朗粒子所需翻越勢阱的能量時，輸出信號淹沒于噪聲信號中，輸出信號的頻譜值減小，反而不會產(chǎn)生隨機共振。

圖9 功率譜值與噪聲強度的關(guān)系圖Fig.9 Relationship between power spectrum value and noise intensity

4 結(jié) 論

本文提出了一種基于Atangana-Baleanu分?jǐn)?shù)階微積分的雙穩(wěn)隨機共振系統(tǒng)模型，利用改進(jìn)的Oustaloup算法對其近似化求解，從仿真實驗中可以得到以下結(jié)論：

(1) 當(dāng)系統(tǒng)中只有微弱周期信號作用時，無需外力作用下分?jǐn)?shù)階求導(dǎo)階次小于臨界階次時系統(tǒng)即可產(chǎn)生隨機共振且臨界階次的大小隨周期信號幅值的增大而增大；

(2) 再將高斯白噪聲加入系統(tǒng)中，噪聲強度一定時改變分?jǐn)?shù)階求導(dǎo)階次，分?jǐn)?shù)階求導(dǎo)階次和輸出信號的頻譜值呈非線性關(guān)系，存在一個最佳的分?jǐn)?shù)階求導(dǎo)階次使系統(tǒng)產(chǎn)生隨機共振；

(3) 固定分?jǐn)?shù)階求導(dǎo)階次不變改變噪聲強度，噪聲強度和輸出信號的頻譜值呈非線性關(guān)系，存在一個最佳的噪聲強度使系統(tǒng)產(chǎn)生隨機共振。

綜上所示證明了該方法的有效性，后續(xù)可以將其應(yīng)用至軸承故障檢測的實際工程中。