基于Duffing振子的微弱信號參數(shù)估計

冀常鵬，許素娜，冀雯靖

(1.遼寧工程技術(shù)大學(xué) 電子與信息工程學(xué)院，遼寧 葫蘆島 125105；2.遼寧省露天礦山裝備工程研究技術(shù)中心，遼寧 阜新 123009；3.遼寧工程技術(shù)大學(xué) 研究生院，遼寧 葫蘆島 125105；4.香港中文大學(xué)(深圳) 理工學(xué)院，廣東 深圳 518172)

0 引 言

隨著信息技術(shù)的迅速發(fā)展，檢測被強噪聲覆蓋的微弱信號[1]逐漸成為一個很重要的課題。Brix[2]首先提出將混沌理論應(yīng)用在微弱信號檢測方面，之后許多專家開始對其展開研究并獲得了諸多成果。在洛倫茲振子、Duffing振子、Rossler振子等混沌振子[3]中，Duffing振子為最典型的一種，因具有初始條件敏感性、零均值噪聲免疫性[4-7]等特征被廣泛應(yīng)用在微弱信號檢測領(lǐng)域。

文獻[8]將待測信號經(jīng)過預(yù)處理后加入到Duffing振子系統(tǒng)進行檢測，實驗結(jié)果表明，此方法能夠提高檢測效果，但未對微弱信號的參數(shù)進行估計；文獻[9]所構(gòu)造的檢測系統(tǒng)仿真模型，能夠有效檢測出淹沒在強噪聲中的微弱正弦信號，獲得相對滿意的仿真效果，但臨界閾值是通過直接觀察相軌跡圖獲得的，缺乏量化判斷依據(jù)；文獻[10]提出將Duffing振子檢測系統(tǒng)中的非線性恢復(fù)力用-x5+x7替換，同時提出尺度變換方法，但其所用的單Duffing振子僅能檢測出和內(nèi)置周期驅(qū)動力頻率非常接近的信號；雖然文獻[11]利用Duffing振子陣列能夠?qū)θ我馕⑷跣盘柕念l率進行估計，但所需振子數(shù)量多，增加了復(fù)雜度。

針對上述文獻不足，提出結(jié)合Duffing振子和頻譜分析的微弱信號頻率檢測方法，在得到頻率的基礎(chǔ)上，給出可以估計微弱信號幅度和初相位的方法，通過Matlab仿真平臺對所提方法進行驗證并完成對微弱信號的有效檢測和高精度參數(shù)估計，仿真過程中用Lyapunov指數(shù)方法觀察檢測系統(tǒng)是否發(fā)生相變。

1 Duffing振子檢測原理和相變判據(jù)

1.1 Duffing振子檢測原理

用于微弱信號檢測的某一Duffing振子[12]為

(1)

(1)式中：k是阻尼系數(shù)；-x(t)+x3(t)是非線性恢復(fù)力；fcos(ωt)是內(nèi)置周期驅(qū)動力，f和ω分別為其幅值和頻率。

經(jīng)過理論分析和大量的仿真實驗可知(1)式所構(gòu)成的檢測系統(tǒng)對微弱方波信號較為敏感，對微弱正弦信號不太敏感，由于我們主要對微弱正弦信號進行檢測，因此，改進(1)式的-x(t)+x3(t)由-x3(t)+x5(t)代替，此時動力學(xué)方程為Homes-Duffing方程

(2)

固定k值，當(dāng)f=0時，可知檢測系統(tǒng)的相平面有3個奇點，原點(0,0)是其鞍點，(+1,0)，(-1,0)為其2個中心相點，檢測系統(tǒng)的相點最終會停留在其中一個，這取決于系統(tǒng)的初始參數(shù)；當(dāng)f≠0時，檢測系統(tǒng)會呈現(xiàn)復(fù)雜的動力學(xué)形態(tài)，若f很小，檢測系統(tǒng)的相點會圍繞(+1,0)或(-1,0)做衰減周期運動，當(dāng)f增加時，檢測系統(tǒng)會依次呈現(xiàn)出同宿軌道、倍周期分叉狀態(tài)，直到f大于fc(混沌閾值)時，檢測系統(tǒng)開始出現(xiàn)混沌狀態(tài)，檢測系統(tǒng)在之后的一段時間都保持混沌狀態(tài)，直到f即將增加到fd(臨界閾值)時，檢測系統(tǒng)呈現(xiàn)圖1a的臨界混沌狀態(tài)，f≥fd時檢測系統(tǒng)呈現(xiàn)圖1b的大尺度周期狀態(tài)。

傳統(tǒng)Duffing振子系統(tǒng)檢測原理為：首先令內(nèi)置周期驅(qū)動力幅值f=f0(f0略小于fd)，此時檢測系統(tǒng)呈現(xiàn)臨界混沌狀態(tài)，若f再增加很小的Δf，檢測系統(tǒng)會轉(zhuǎn)為大尺度周期狀態(tài)；然后將待測信號作為內(nèi)置周期驅(qū)動力的補充加入到檢測系統(tǒng)，若待測信號中含有與內(nèi)置周期驅(qū)動力同頻的微弱正弦信號且疊加后幅值大于等于fd，則檢測系統(tǒng)會發(fā)生相變，即待測信號中存在微弱正弦信號，否則不存在；檢測過程中常用直接觀察法判斷系統(tǒng)是否發(fā)生相變。

傳統(tǒng)Duffing振子檢測系統(tǒng)雖然可以有效檢測出某些微弱正弦信號，但存在以下缺點。

1)若使用單一Duffing振子，則僅能檢測和內(nèi)置周期驅(qū)動力頻率非常接近的微弱正弦信號，即微弱正弦信號頻率必須事先知道。當(dāng)微弱正弦信號頻率未知時，雖然可以利用Duffing振子陣列完成檢測，但增加了復(fù)雜度。

2)目前判斷檢測系統(tǒng)是否產(chǎn)生相變的常用方法為直接觀察法，該方法由于是通過肉眼觀察，故存在一定的主觀性、效率低且容易發(fā)生誤判，很難滿足實際工程的需求。

基于以上相關(guān)問題，提出將頻譜分析和Duffing振子相結(jié)合的方法來估計微弱正弦信號的頻率，并利用Lyapunov指數(shù)方法[13]作為檢測系統(tǒng)產(chǎn)生相變的依據(jù)。在該基礎(chǔ)上完成了微弱正弦信號幅度和初相位的估計。

1.2 相變判據(jù)

混沌檢測是通過Duffing振子系統(tǒng)產(chǎn)生相變來檢測微弱信號的，因此，準(zhǔn)確判斷檢測系統(tǒng)的狀態(tài)至關(guān)重要。雖然Melnikov方法[14]可以作為產(chǎn)生相變的判據(jù)，但精確度不是很高，同時為避免直接觀察法的缺點，采用Lyapunov指數(shù)方法來判斷檢測系統(tǒng)是否產(chǎn)生相變。

Lyapunov指數(shù)方法判斷Duffing振子系統(tǒng)狀態(tài)的原理為若檢測系統(tǒng)的最大Lyapunov指數(shù)值Lyapunovmax>0，代表系統(tǒng)呈現(xiàn)混沌狀態(tài)。若檢測系統(tǒng)的Lyapunovmax<0，則代表其呈現(xiàn)大尺度周期狀態(tài)，因此，當(dāng)Lyapunovmax>0跳變?yōu)長yapunovmax<0時所對應(yīng)的內(nèi)置周期驅(qū)動力幅值為檢測系統(tǒng)臨界閾值fd。

利用Wolf法求時間序列的Lyapunovmax的基本思想是對于混沌時間序列x(t)，對其進行m維的相空間重構(gòu)，即吸引子上的某點可以由{x(t),x(t+τ),…,x(t+[m-1]τ)}獲得，假如選擇點{x(t0),x(t0+τ),…,x(t0+[m-1]τ)}作為基準(zhǔn)軌道上的初始點，然后再選擇距離該點最近一點作為鄰近軌道的初始點，用L(t0)表示它們之間的距離。當(dāng)時間為t1時，距離變?yōu)長′(t1)，為避免其溢出，再選擇另一數(shù)據(jù)點，但應(yīng)滿足2個條件：①與基準(zhǔn)軌道上的對應(yīng)演化點距離要非常??；②基準(zhǔn)軌道上對應(yīng)點到新數(shù)據(jù)點和到舊數(shù)據(jù)點所構(gòu)成的2向量夾角也必須非常小。若不存在滿足以上2個條件的點，就繼續(xù)使用原數(shù)據(jù)點，不進行更新，直到基準(zhǔn)軌道取遍全部的數(shù)據(jù)，因此，最大Lyapunov指數(shù)λmax為

(3)

(3)式中，M是替代的總步數(shù)。

2 改進方法

2.1 頻譜法Duffing振子檢測頻率

頻譜法Duffing振子檢測頻率的思想是在傳統(tǒng)Duffing振子檢測系統(tǒng)基礎(chǔ)上，摒棄內(nèi)置周期驅(qū)動力，將待測信號直接輸入到檢測系統(tǒng)，然后對檢測系統(tǒng)的輸出信號作頻譜分析獲得微弱正弦信號的頻率。

結(jié)合傳統(tǒng)Duffing振子所構(gòu)造的檢測模型為

(4)

(4)式中，s(t)=acos(ωt)+n(t)是待測信號，其中，n(t)是高斯白噪聲。因為沒有內(nèi)置周期驅(qū)動力，此時檢測系統(tǒng)只由待測信號驅(qū)動，所以不用考慮微弱正弦信號與內(nèi)置周期驅(qū)動力是否同頻，故可以檢測任意頻率的微弱正弦信號，避免了利用多振子陣列檢測增加復(fù)雜度的問題。為滿足當(dāng)僅有待測信號存在時，檢測系統(tǒng)能夠發(fā)生相變，待測信號幅度a需大于臨界閾值fd。

在進行頻譜分析時，有2種方法可以選擇，分別為傅里葉分析方法和小波分析方法[15]。由傅里葉分析方法得到的頻譜圖，待測信號的幅度和其他諧波分量幅度之間存在明顯的差異，因此,能求出待測信號頻率(即微弱正弦信號的頻率)。

將待測信號輸入到Duffing振子檢測系統(tǒng)，對輸出量作傅里葉變換，得到以頻率為橫坐標(biāo)，各諧波幅度值為縱坐標(biāo)的頻譜圖，由頻譜圖能求出待測信號的頻率。

對所提方法進行實驗驗證，(4)式中k值取0.5，ω值取10 rad/s，參數(shù)設(shè)置完成后開始仿真，傅里葉分析通過編程實現(xiàn)，并以坐標(biāo)形式給出實驗數(shù)據(jù)，檢測系統(tǒng)的輸出信號頻譜圖如圖2。

2.2 待測信號幅度和初相位估計方法

由2.1節(jié)可完成待測信號頻率的估計，在此基礎(chǔ)上，繼續(xù)對待測信號幅度和初相位進行估計，此時Duffing振子檢測系統(tǒng)內(nèi)置周期驅(qū)動力不置零，其頻率設(shè)置為待測信號頻率，然后結(jié)合Lyapunov指數(shù)方法求出臨界閾值fd，該方法檢測模型為

(5)

(5)式中：f0cos(ωt+θ)為檢測系統(tǒng)內(nèi)置周期驅(qū)動力；acos(ωt+φ)+n(t)是待測信號。此時,檢測系統(tǒng)呈現(xiàn)臨界混沌狀態(tài)，將含有與內(nèi)置周期驅(qū)動力同頻的微弱信號加入到檢測系統(tǒng)，檢測系統(tǒng)會呈現(xiàn)大尺度周期狀態(tài)。通過理論研究(5)式的Duffing振子進入大尺度周期狀態(tài)的條件，若能發(fā)生相變，則(5)式中總的周期驅(qū)動力幅值應(yīng)滿足

(6)

若用fm代表檢測系統(tǒng)剛好發(fā)生相變時所對應(yīng)的內(nèi)置周期驅(qū)動力幅值，則

(7)

(7)式中，fm，fd均可得到，若還可以建立一個有關(guān)a和φ的類似(7)式的方程，則通過解二元方程組可完成對待測信號的幅度和初相位的估計。以此為出發(fā)點，通過改變檢測系統(tǒng)的內(nèi)置周期驅(qū)動力的初相位得到2個方程組。具體思想是首先用頻譜法Duffing振子在完成待測信號的頻率估計后，然后設(shè)置(5)式中的檢測系統(tǒng)參數(shù)，將待測信號分別輸入到內(nèi)置周期驅(qū)動力初相位為0和π的Duffing振子檢測系統(tǒng)中，求出這2個系統(tǒng)發(fā)生相變時所對應(yīng)的內(nèi)置周期驅(qū)動力幅值，分別用f1和f2表示。結(jié)合(7)式得到的2個方程組為

(8)

解方程組(8)可得

(9)

由(9)式能求出a和φ，至此完成了對待測信號的頻率、幅度和初相位的估計。

3 仿真實驗和結(jié)果分析

若對系統(tǒng)的輸出信號用離散傅里葉變換(discrete Fourier transformation, DFT)算法進行頻譜分析，則復(fù)雜度為序列長度的平方，即o(n2)。對DFT進行改進，形成一套高效計算方法，即快速傅里葉變換算法(fast Fourier transformation, FFT)。其思想是將長序列分解為若干短序列進行DFT計算，然后通過若干旋轉(zhuǎn)因子的復(fù)數(shù)乘法和復(fù)數(shù)加法合成最終的結(jié)果。FFT算法中，若序列長度是2的冪次，可將序列長為N的DFT分割為2個長為N/2的子序列的DFT，稱為基2-FFT?？焖俑道锶~變換算法僅需要o(nlogn)的計算復(fù)雜度，在此，利用FFT對系統(tǒng)的輸出信號進行頻譜分析。

將待測信號s1(t)，s2(t)和s3(t)分別加入到內(nèi)置周期驅(qū)動力為0的Duffing振子檢測系統(tǒng)，對其輸出信號作頻譜分析。給出這3組實驗利用快速傅里葉變換算法得到頻譜圖的運行時間，如表1，輸出信號的頻譜圖如圖3。

由圖3并結(jié)合2.1節(jié)相關(guān)知識對待測信號頻率進行估計，同時利用Duffing振子陣列方法對待測信號s1(t)，s2(t)和s3(t)進行檢測，根據(jù)這3個待測信號的頻率范圍，將第1個振子頻率取1，第2個振子頻率取1.03，即后面振子頻率為其前一個振子頻率的1.03倍，依次類推，因為檢測頻率在1～5 rad，由1.0354=4.934，1.0355=5.082，故需設(shè)置55個振子。

表1 本文算法檢測時間Tab.1 Algorithm detection time in this paper

通過觀察時序圖,發(fā)現(xiàn)其中3組間歇混沌周期最長，分別為第1個和第2個振子、第23和24振子和第46和47個振子，如圖4。

由于第1個和第2個振子分別所對應(yīng)的間歇混沌周期大致相等，因此，取它們對應(yīng)頻率的平均值作為待測信號的近似值，同理，計算第23和第24個振子以及第46和第47個振子所對應(yīng)頻率的平均值。最后將頻譜法Duffing振子和Duffing振子陣列方法對待測信號頻率估計結(jié)果進行對比，如表2。

由表2可知，Duffing振子在信噪比為-43.01 dB的條件下仍能有效檢測出微弱正弦信號，在對信號頻率進行估計時，頻譜法Duffing振子和Duffing振子陣列方法相比具有較高的精度，相對誤差僅為10-3，復(fù)雜度低，易操作等優(yōu)點。

表2 2種方法頻率估計結(jié)果Tab.2 Two methods of frequency estimation results

在得到待測信號的頻率估計值后，繼續(xù)對待測信號的幅度和初相位進行估計，以待測信號s1(t)為例，待測信號s2(t)和s3(t)檢測原理與s1(t)相同，僅需將檢測系統(tǒng)的內(nèi)置周期驅(qū)動力角頻率設(shè)置為2.000 566 rad/s和3.998 619 rad/s。

具體實現(xiàn)步驟如下。

步驟1 將Duffing振子檢測系統(tǒng)的內(nèi)置周期驅(qū)動力角頻率設(shè)置成1.000 911 rad/s，及其他系統(tǒng)參數(shù)設(shè)置完成后利用Lyapunov指數(shù)方法得到臨界閾值fd。

步驟2 將待測信號s1(t)作為內(nèi)置周期驅(qū)動力的補充，輸入初相位分別為0和π的Duffing振子檢測系統(tǒng)。

步驟3 通過編程達(dá)到以下功能：將步驟2中所得到的輸出時間序列，利用FFT算法來計算其平均周期；通過C-C方法確定其延遲時間和嵌入維數(shù)并進行相空間重構(gòu)；最后利用重構(gòu)的相空間，再次計算檢測系統(tǒng)的最大Lyapunov指數(shù)，若其值大于0，則調(diào)整內(nèi)置周期驅(qū)動力的幅值，并重新計算調(diào)整后的最大Lyapunov指數(shù)，直到最大Lyapunov指數(shù)剛好由大于0轉(zhuǎn)變?yōu)樾∮?，記錄此時所對應(yīng)的內(nèi)置周期驅(qū)動力幅值f。

步驟4 初相位為0和π的檢測系統(tǒng)發(fā)生相變時所對應(yīng)的內(nèi)置周期驅(qū)動力幅值分別記為f1和f2，代入(9)式后，得到待測信號s1(t)的幅度和初相位的估計值。

先執(zhí)行步驟1，并給出檢測系統(tǒng)未加入待測信號時，最大Lyapunov指數(shù)和內(nèi)置周期驅(qū)動力幅值f的關(guān)系，如圖5，其中，f為歸一化值。

由圖5可知，內(nèi)置周期驅(qū)動力角頻率為1.000 911 rad/s時所對應(yīng)的fd=0.773，角頻率為2.000 566 rad/s時所對應(yīng)的fd=0.747 1，角頻率為3.998 619 rad/s時所對應(yīng)的fd=0.701。

然后執(zhí)行步驟2—4，最后求出當(dāng)待測信號為s1(t)時所對應(yīng)的f1=0.765 896 5，f2=0.780 038 7；為s2(t)時所對應(yīng)的f1=0.729 712，f2=0.764 353；為s3(t)時所對應(yīng)的f1=0.674 66，f2=0.724 66。再結(jié)合步驟1的結(jié)果，代入(9)式經(jīng)計算得到待測信號s1(t)，s2(t)和s3(t)的幅度，與采用直接觀察法得到的幅度估計值進行對比，結(jié)果如表3。對初相位的計算結(jié)果如表4。

表3 2種方法幅度估計結(jié)果Tab.3 Two methods of amplitude estimation results

由表3可看出，本文方法和直接觀察法相比，對待測信號幅度的估計誤差更小，相對誤差僅為10-4。

由表4可看出，本文方法對待測信號初相位的估計誤差很小，相對誤差僅為10-3。

表4 相位估計結(jié)果Tab.4 Phase estimation results

4 結(jié) 論

單一Duffing振子僅能檢測與內(nèi)置周期驅(qū)動力頻率十分接近的微弱信號，本文從這個問題出發(fā)，提出頻譜式Duffing振子方法，該方法能夠只用單個Duffing振子實現(xiàn)未知頻率信號的有效檢測和頻率估計，且精度可以達(dá)到10-3；在完成頻率估計的基礎(chǔ)上，還給出微弱信號幅度和初相位的估計方法，同時為避免直接觀察法的主觀性、效率低和Melnikov方法的精度低等缺點，采用Lyapunov指數(shù)方法作為檢測系統(tǒng)相變判據(jù)；最后進行3組仿真實驗，通過對仿真數(shù)據(jù)的分析，得出所提方法在信噪比為-43.01 dB的條件下，仍能完成對未知頻率且具有任意初相位的正弦信號的檢測，并對3項參數(shù)具有較高的估計精度。