Atangana-Baleanu分?jǐn)?shù)階雙穩(wěn)系統(tǒng)的隨機(jī)共振現(xiàn)象

汪洋百慧， 鄭永軍， 羅 哉

(中國(guó)計(jì)量大學(xué) 計(jì)量測(cè)試工程學(xué)院，浙江 杭州 310018)

1 引 言

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

隨機(jī)共振通常發(fā)生在布朗粒子反常擴(kuò)散過(guò)程中，具有非馬爾可夫特性。由于分?jǐn)?shù)階微積分特點(diǎn)是在時(shí)間上的記憶性和空間上的關(guān)聯(lián)性[11～19]，因此可用來(lái)描述這類(lèi)現(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ù)階微積分，有效解決了這一問(wèn)題[21]。其中借助廣義Mittag-Leffler(M-L)函數(shù)作為非局部和非奇異核，使其具有均方位移的交叉性質(zhì)并保證了信息在延展的開(kāi)始和結(jié)束時(shí)的一致性[22～24]。

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

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

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

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

(1)

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

(2)

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

(3)

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

(4)

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

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

(5)

上述定義將有助于討論現(xiàn)實(shí)世界中的問(wèn)題，并且在使用拉普拉斯變換解決初始條件下的某些物理問(wèn)題時(shí)也將具有很大優(yōu)勢(shì)。經(jīng)過(guò)簡(jiǎn)單計(jì)算后得到基于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)可以寫(xiě)為式(9)所示的傳遞函數(shù)的形式：

(9)

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

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

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

H(s)=sα

(10)

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

(11)

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

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

(12)

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

(13)

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

(14)

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

(15)

根據(jù)實(shí)零點(diǎn)和極點(diǎn)的遞歸分布，第k個(gè)零點(diǎn)和極點(diǎn)可以寫(xiě)為式(16)所示：

(16)

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

(17)

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

(18)

經(jīng)過(guò)大量實(shí)驗(yàn)驗(yàn)證，一般取b=10，d=9。

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

3.1 隨機(jī)共振的基本原理

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

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

(19)

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

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

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

(20)

式中：D為噪聲強(qiáng)度；ΔV為勢(shì)壘高度。

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

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

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

(21)

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

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

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

(22)

(23)

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

(24)

由于阻尼力和噪聲都源于布朗粒子與周?chē)橘|(zhì)分子之間的隨機(jī)碰撞，因此隨機(jī)力ξ(t)和阻尼核函數(shù)γ(t)滿(mǎn)足以下漲落耗散定理:

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

(25)

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

(26)

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

(27)

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

(28)

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

(29)

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

(30)

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

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

(31)

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

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

如圖3(a)所示為基于A-B分?jǐn)?shù)階導(dǎo)數(shù)時(shí)，分?jǐn)?shù)階求導(dǎo)階次以0.1的步長(zhǎng)從0.1遞增到0.9的輸出響應(yīng)的時(shí)域圖。可以看出α的變化改變了勢(shì)壘高度的大小并使勢(shì)阱點(diǎn)發(fā)生漂移。當(dāng)α大于0.6時(shí)，布朗粒子只能在x=1或x=-1的單一勢(shì)阱處運(yùn)動(dòng)，說(shuō)明其沒(méi)能越過(guò)勢(shì)壘。減小α，當(dāng)α小于0.6時(shí)，可以看出布朗粒子的運(yùn)動(dòng)軌跡發(fā)生了改變，在x=1和x=-1之間運(yùn)動(dòng)。說(shuō)明α值越小，系統(tǒng)的內(nèi)核函數(shù)在時(shí)間上的記憶性越強(qiáng)，布朗粒子就有足夠大的能量可以越過(guò)勢(shì)壘。

圖3 A=0.3時(shí)α變化下系統(tǒng)輸出時(shí)域圖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)的臨界值時(shí)，需要一些額外的能量來(lái)驅(qū)動(dòng)布朗粒子克服勢(shì)壘高度；當(dāng)α小于該臨界值時(shí)無(wú)需外力驅(qū)使即可產(chǎn)生隨機(jī)共振。

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

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

如圖5所示，可以看出系統(tǒng)的臨界階次受信號(hào)幅值的影響，信號(hào)幅值越大臨界階次越大。這是因?yàn)樾盘?hào)幅值的大小會(huì)影響勢(shì)壘高度的大小，幅值越大勢(shì)壘高度越大，則系統(tǒng)的臨界階次也就越大。

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

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

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

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

圖6 外加噪聲時(shí)D不變，α變化下系統(tǒng)輸出的時(shí)域圖和頻域圖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)的噪聲強(qiáng)度固定，當(dāng)α小于0.7時(shí)，功率譜值均小于0.02，甚至還未達(dá)到原始信號(hào)的功率譜；增加α的大小，功率譜值也隨之增大并達(dá)到一個(gè)極大值，在此極大值處，雙穩(wěn)態(tài)系統(tǒng)產(chǎn)生了隨機(jī)共振；繼續(xù)增加α的大小，系統(tǒng)的勢(shì)壘值也在不斷增大，布朗粒子很難再翻越勢(shì)阱，此時(shí)功率譜值反而減小，抑制了隨機(jī)共振的產(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)階次時(shí)改變?cè)肼晱?qiáng)度

由4.2節(jié)中可知噪聲強(qiáng)度固定，α過(guò)小或過(guò)大時(shí)，雙穩(wěn)態(tài)系統(tǒng)均不會(huì)產(chǎn)生隨機(jī)共振?；诖饲闆r，本節(jié)將討論當(dāng)α=0.75時(shí)，加入噪聲并且改變?cè)肼晱?qiáng)度的隨機(jī)共振現(xiàn)象。

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

圖8 外加噪聲時(shí)α不變，D變化下系統(tǒng)輸出的時(shí)域圖和頻域圖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)階次固定，增大噪聲強(qiáng)度，雙穩(wěn)態(tài)系統(tǒng)可以產(chǎn)生隨機(jī)共振，但持續(xù)增大噪聲強(qiáng)度反而會(huì)抑制隨機(jī)共振的產(chǎn)生。當(dāng)噪聲提供的能量遠(yuǎn)大于布朗粒子所需翻越勢(shì)阱的能量時(shí)，輸出信號(hào)淹沒(méi)于噪聲信號(hào)中，輸出信號(hào)的頻譜值減小，反而不會(huì)產(chǎn)生隨機(jī)共振。

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

4 結(jié) 論

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

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

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

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

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