Duffing振子變尺度未知多頻微弱信號(hào)提取方法研究

黃繼堯，陳長(zhǎng)興，凌云飛，藺向陽

(空軍工程大學(xué)基礎(chǔ)部，陜西西安 710043)

0 引言

Duffing振子是一種用于信號(hào)處理的非線性混沌振子模型，它在混沌臨界狀態(tài)下對(duì)噪聲具有免疫性，但對(duì)初始條件極為敏感[1]。可以利用Duffing振子的這些特性來檢測(cè)強(qiáng)噪聲背景下的微弱信號(hào)。自20世紀(jì)后期引入混沌振子系統(tǒng)以檢測(cè)微弱信號(hào)[2-4]以來，諸多學(xué)者在此領(lǐng)域進(jìn)行了深入的研究。謝濤等證明了混沌振子檢測(cè)微弱信號(hào)的可靠性，為混沌檢測(cè)微弱信號(hào)提供了理論論據(jù)[5]。

李月等分析了不同類型噪聲對(duì)混沌系統(tǒng)的影響，證明了混沌系統(tǒng)對(duì)噪聲的免疫性[6]；秦衛(wèi)陽等分析了Duffing方程的派生系統(tǒng)在同步過程中能夠?qū)崿F(xiàn)檢測(cè)2種信號(hào)中的微小差別[7]；劉海波等提出一種利用單一振子實(shí)現(xiàn)以Duffing振子大周期態(tài)為基礎(chǔ)的測(cè)量頻率的方法[8]。常見的混沌狀態(tài)判定方法有Lyapunov指數(shù)法、龐加萊截面法、分維數(shù)計(jì)算法、Kolmogorov熵法、Melnikov法等[9]。與傳統(tǒng)的線性信號(hào)檢測(cè)方法相比，基于Duffing振子的微弱信號(hào)檢測(cè)方法可以實(shí)現(xiàn)非常低的信噪比(SNR)，信噪比門限最低可達(dá)到-100 dB，并可得到傳統(tǒng)檢測(cè)方法無法得到的微弱待測(cè)信號(hào)，這為低信噪比的微弱信號(hào)檢測(cè)提供了一種新的方法[10]。

然而目前大多數(shù)基于Duffing振子的微弱信號(hào)檢測(cè)研究都是針對(duì)單頻諧波信號(hào)的檢測(cè)識(shí)別[11-14]，而很少有研究關(guān)注微弱多頻信號(hào)的檢測(cè)。因此對(duì)基于Duffing振子的微弱多頻信號(hào)檢測(cè)進(jìn)行深入研究是十分有必要的。本文提出了一種結(jié)合信號(hào)頻率截取預(yù)處理和Duffing振子檢測(cè)微弱多頻信號(hào)的方法。

1 Duffing振子微弱信號(hào)檢測(cè)原理

Holmes型Duffing方程[15]一般形式為

(1)

式中：k為阻尼比；F為周期性策動(dòng)力幅值；ω為角頻率；Fcosωt為周期性策動(dòng)力。

當(dāng)k保持不變時(shí)，隨著F的增大，存在一個(gè)臨界策動(dòng)力振幅Fd。當(dāng)F小于Fd時(shí)，系統(tǒng)處于混沌狀態(tài)；當(dāng)F大于Fd時(shí)，系統(tǒng)處于大周期態(tài)。在Duffing振子方程式(1)中，臨界值Fd≈0.825?；煦绾痛笾芷趹B(tài)的變化及策動(dòng)力臨界值如圖1所示。

(a)混沌態(tài)

(b)周期態(tài)

(c)臨界點(diǎn)圖1 系統(tǒng)歷經(jīng)的狀態(tài)變化及分岔臨界圖

當(dāng)Duffing系統(tǒng)處于從混沌狀態(tài)轉(zhuǎn)變?yōu)榇笾芷趹B(tài)的臨界狀態(tài)時(shí)，就統(tǒng)計(jì)特性而言，它不受高斯白噪聲的影響。由于對(duì)參數(shù) 的敏感性和對(duì)噪聲的免疫性，式(1)可以進(jìn)一步表示為

(2)

式中：s(t)為待測(cè)信號(hào)；n(t)為噪聲；s(t)+n(t)為帶有噪聲的待測(cè)信號(hào)。

由式(2)可知，當(dāng)式(1)中的F略小于臨界幅值Fd時(shí)，s(t)+n(t)輸入Duffing振子系統(tǒng)，n(t)作為擾動(dòng)項(xiàng)，則可以使用式(2)檢測(cè)待測(cè)信號(hào)。sn(t)=s(t)+n(t)是含噪待測(cè)信號(hào)。當(dāng)s(t)包含角頻率為ω的信號(hào)分量時(shí)，系統(tǒng)將從混沌狀態(tài)轉(zhuǎn)換為大周期態(tài)，這可以從系統(tǒng)響應(yīng)的相平面軌跡獲得，然后便可以獲得待測(cè)信號(hào)的頻率。

傅里葉變換是一種分析信號(hào)的方法，可以把具有周期特征的非正弦信號(hào)分解為多個(gè)疊加的正弦信號(hào)，然后通過對(duì)時(shí)間連續(xù)信號(hào)的等間隔采樣，得到采樣頻率fs，并把采樣值轉(zhuǎn)換為數(shù)字序列。連續(xù)時(shí)間序列x(t)與其對(duì)應(yīng)的連續(xù)非周期頻譜函數(shù)X(f)之間的轉(zhuǎn)換關(guān)系為：

(3)

當(dāng)設(shè)定一個(gè)待測(cè)信號(hào)為hcos(ωt+φ)，其采樣頻率為fs。ω為確定的角頻率，φ是初始相位。假設(shè)采樣信號(hào)向坐標(biāo)軸左側(cè)移動(dòng)dt=1/fs。根據(jù)公式T/dt=2π/dθ，就可以得到dθ=dt·2π/T=ω/fs，T為信號(hào)的周期，dθ為對(duì)應(yīng)相位的增量。因此，將采樣信號(hào)向左挪動(dòng)時(shí)間j等同于將待測(cè)信號(hào)由hcos(ωt+φ)變成hcos(ωt+φ+jdθ)。

在不考慮噪聲，即n(t)=0的情況下，設(shè)定待測(cè)信號(hào)尺度變換后(尺度變換系數(shù)R=ω)與策動(dòng)力有相同的角頻率，則檢測(cè)待測(cè)信號(hào)的Duffing方程為

(4)

改變式(6)右側(cè)項(xiàng)的形式

(5)

從式(4)中，可以推出只有當(dāng)(φ+jdθ)=2mπ(m是整數(shù))時(shí)，與待測(cè)信號(hào)和策動(dòng)力相加的待測(cè)信號(hào)的振幅是最大的。

文獻(xiàn)[16]提出了一種基于Duffing振子的尺度變換微弱信號(hào)檢測(cè)模型，有效地解決了這些問題。檢測(cè)模型如式(6)所示：

(6)

這4個(gè)方程的輸入信號(hào)初始相位的檢測(cè)范圍能夠覆蓋整段頻率(-π，π)，避免了單一方程只能檢測(cè)某一段初始相位的限制。

式(6)解決了基于Duffing振子微弱信號(hào)檢測(cè)的一些關(guān)鍵問題，然而目前僅對(duì)含噪聲的單頻諧波信號(hào)進(jìn)行分析，而不涉及實(shí)際工程應(yīng)用中常見的微弱多頻信號(hào)的檢測(cè)。因此本文將研究基于Duffing振子的微弱多頻信號(hào)檢測(cè)。

2 基于頻率截獲Duffing振子的微弱多頻信號(hào)檢測(cè)

2.1 信號(hào)頻率截取方法

為了將式(6)中的檢測(cè)模型應(yīng)用于微弱多頻信號(hào)，提出了一種信號(hào)頻率截取的預(yù)處理方法。其原理如下:根據(jù)Shannon采樣定理，試驗(yàn)信號(hào)sn(t)的采樣頻率為fs，則其截止頻率為fs/2，為了便于后續(xù)分析，首先在頻域中截取待測(cè)信號(hào)。截取頻率的下限為fmin，上限為fmax，其滿足0

首先，利用快速傅里葉變換(FFT)將試驗(yàn)信號(hào)sn(t)變換為頻域信號(hào)sn(f);其次，信號(hào)sn′(f)是由sn(f)通過頻率截取得到的，即離散譜數(shù)據(jù)的頻率范圍為(fmin,fmax);最后，將信號(hào)sn′(f)通過快速傅里葉逆變換(IFFT)轉(zhuǎn)變?yōu)闀r(shí)域信號(hào)sn′(t)。

上述頻率截獲過程可以表示為:

(7)

(8)

(9)

頻率截取過程后，sn′(t)只包含滿足fmin

2.2 待測(cè)信號(hào)提取模型

將信號(hào)頻率截獲方法與式(6)檢測(cè)模型相結(jié)合，針對(duì)微弱多頻信號(hào)，提出了一種新的待測(cè)信號(hào)提取模型。其過程具體如下:

(1)對(duì)原始待測(cè)信號(hào)進(jìn)行加噪處理，得到被噪聲污染的待測(cè)信號(hào)sn(t)，采樣頻率為fs， 由強(qiáng)度為D的背景噪聲和微弱待測(cè)信號(hào)s1(t),s2(t),…,sn(t)構(gòu)成，微弱待測(cè)信號(hào)的角頻率和振幅分別為ω1,ω2,…,ωn和h1,h2，…,hn。

(2)確定原始待測(cè)信號(hào)的頻率測(cè)試范圍[ωmin,ωmax]。然后設(shè)定ω=ωmin，接著設(shè)置頻率增量dω，令ωi+1=ωi+dω。

(3)然后對(duì)待測(cè)信號(hào)sn(t)進(jìn)行傅里葉變換后得到信號(hào)sn(f)再進(jìn)行頻段為[ωmin,ωmax]的頻率截取預(yù)處理。新的信號(hào)在預(yù)處理后為sn′(f) ;

(4)設(shè)定尺度系數(shù)Ri=ωi(即ω的值每增長(zhǎng)一次，R的值都會(huì)相應(yīng)的增長(zhǎng)并與之相等)，然后對(duì)信號(hào)進(jìn)行頻率/時(shí)間尺度變換，之后得到在時(shí)間尺度t′下新的待測(cè)信號(hào)sn′(t)。變換后的采樣頻率為fsr=fs/Ri；

(5)頻率截取預(yù)處理和尺度變換后，將新獲得的信號(hào)sn′(t)輸入到式(6)中。然后利用四階Runge-Kutta算法求解各方程，繪制各輸出響應(yīng)的相平面軌跡圖。只要任何一個(gè)相平面軌跡處于大周期態(tài)，則意味著待測(cè)信號(hào)sn(t)中的頻率分量ωk被檢測(cè)到，然后重復(fù)迭代過程，直到頻率滿足ω≥ωmax為止，過程如圖2所示。

圖2 信號(hào)檢測(cè)框圖

從而檢測(cè)到該頻率范圍內(nèi)的所有信號(hào)頻率分量，實(shí)現(xiàn)對(duì)未知多頻率信號(hào)的檢測(cè)。如果dω太小，檢測(cè)未知信號(hào)的計(jì)算量將大幅增加。選擇角頻率為ω=20 rad/s的待測(cè)信號(hào)進(jìn)行測(cè)試，設(shè)定Δω為0.2 rad/s。選擇不同的Δω的倍數(shù)與計(jì)算時(shí)間的關(guān)系如圖3所示，選擇dω為2Δω有較快的運(yùn)算速度，因此選擇dω為2Δω。

圖3 不同的Δω倍數(shù)與計(jì)算時(shí)間的關(guān)系圖

所以當(dāng)檢測(cè)過程完成后，可以覆蓋完整的頻段[ωmin,ωmax]。迭代過程如圖4所示：

圖4 未知多頻信號(hào)檢測(cè)的迭代過程

3 仿真與分析

設(shè)定待測(cè)信號(hào)的方程為

(10)

式中：ω1,ω2,ω3和h1,h2,h3分別為sn(t)中包含的3個(gè)頻率分量的角頻率和振幅；g(t)為均值為0，方差為1的高斯白噪聲。

設(shè)定待測(cè)信號(hào)的初始相位為0，待測(cè)信號(hào)的采樣頻率fs=100 Hz，h1=h2=h3=0.1，ω1、ω2、ω3分別為20，30，50 rad/s，噪聲強(qiáng)度D=1，Δω為0.2 rad/s。設(shè)定[ωmin,ωmax]為(0,60)。繪制待測(cè)信號(hào)的時(shí)域圖和頻域圖如圖5所示，波形顯示3 000個(gè)采樣點(diǎn)。

由于在實(shí)際工程問題中噪聲一般是不可忽略的，因此在時(shí)域圖和頻域圖中，可以輕易發(fā)現(xiàn)由于噪聲的污染，很難找到待測(cè)信號(hào)的頻率的存在。而相圖的檢測(cè)也同樣容易受到噪聲的干擾，如圖6所示。

將ω=20 rad/s的待測(cè)信號(hào)送入式(6)中的第一個(gè)方程，可以發(fā)現(xiàn)圖6(a)所顯示的相圖在噪聲的影響下，較難分辨是出于混沌狀態(tài)還是大周期態(tài)。而圖6(b)所顯示的則是在經(jīng)過尺度變換后消除了噪聲的影響，相圖是處于大周期態(tài)的。在下文中，都將以經(jīng)過尺度變換后的相圖闡述。

采用本文提出的優(yōu)化模型進(jìn)行檢測(cè)，在ω=10 rad/s,20 rad/s,30 rad/s,40 rad/s,50 rad/s,60 rad/s，6個(gè)點(diǎn)處進(jìn)行取值，如圖7所示。

(a)時(shí)域圖

(a)未經(jīng)過尺度變換的相平面軌跡

(b)經(jīng)過尺度變換的相平面軌跡 圖6 ω=20 rad/s的待測(cè)信號(hào)的相平面軌跡

通過實(shí)驗(yàn)可以得知，圖7(b)、(c)、(e)所顯示的是分別為ω=20 rad/s、ω=30 rad/s和ω=50 rad/s的頻率分量所對(duì)應(yīng)的相圖，相圖是明顯處于大周期態(tài)的。而圖7(a)、(d)、(f)所顯示的分別是在ω=10 rad/s、ω=40 rad/s和ω=60 rad/s的頻點(diǎn)所對(duì)應(yīng)的相圖，相圖是處于混沌狀態(tài)的，符合實(shí)驗(yàn)設(shè)定的預(yù)期條件。

實(shí)驗(yàn)發(fā)現(xiàn)：

(1)未使用尺度變換，就很難將噪聲污染的相圖分辨出來，存在對(duì)相圖的誤判可能，尺度變換可以有效地消除噪聲對(duì)相圖的干擾，尺度系數(shù)R的存在有利于待測(cè)信號(hào)的提取。

(2)使用本方法成功檢測(cè)到了ω=20 rad/s,ω=30 rad/s，ω=50 rad/s的頻率分量。而在ω=60 rad/s這兩個(gè)頻點(diǎn)則相圖未發(fā)生狀態(tài)轉(zhuǎn)移，證明不存在該頻率的分量，具有普遍性和實(shí)用性。

(a)ω=10 rad/s的相平面軌跡

(b)ω=20 rad/s的相平面軌跡

(c)ω=30 rad/s的相平面軌跡

(d)ω=40 rad/s的相平面軌跡

(e)ω=50 rad/s的相平面軌跡

(f)ω=60 rad/s的相平面軌跡

(3)本方法實(shí)現(xiàn)的微弱信號(hào)檢測(cè)最低信噪比門限可以達(dá)到-24.2 dB，低于文獻(xiàn)[16]中提出的原檢測(cè)模型最低信噪比門限的-13.01 dB。

(4)Δω的選取如果過大，則會(huì)出現(xiàn)漏檢測(cè)現(xiàn)象；Δω的選取如果過小，則檢測(cè)效率就較低，因此需要找一個(gè)平衡點(diǎn)。一般采用比待測(cè)信號(hào)頻率低3個(gè)數(shù)量級(jí)的Δω可以實(shí)現(xiàn)限定頻率范圍內(nèi)全部頻率的檢測(cè)。本文所提出的方法在實(shí)際實(shí)驗(yàn)中采用的Δω是0.2 rad/s，檢測(cè)到了3個(gè)待測(cè)信號(hào)的頻率。

(5)本文所提出的方法與Duffing振子并行處理方式的區(qū)別在于，本文所提出的方法是對(duì)頻域加帶通濾波器的基礎(chǔ)上進(jìn)行而且僅用一個(gè)Duffing振子便可檢測(cè)到限定頻域內(nèi)全部的待測(cè)信號(hào)，而Duffing振子并行處理方式一般采用陣列的方式進(jìn)行檢測(cè)，即需要較多振子。本文提出的方法節(jié)省了振子的數(shù)量，對(duì)運(yùn)算條件需求較低，在實(shí)際實(shí)驗(yàn)處理中也可以達(dá)到較為理想的效果。

4 結(jié)束語

本文研究了一種未知多頻微弱信號(hào)檢測(cè)方法，該方法針對(duì)原檢測(cè)模型只能同時(shí)檢測(cè)某一頻率的待測(cè)信號(hào)，克服了噪聲誘導(dǎo)對(duì)檢測(cè)結(jié)果的不利影響，提出了一種可同時(shí)檢測(cè)多個(gè)未知頻率微弱信號(hào)的Duffing振子優(yōu)化模型，可以實(shí)現(xiàn)提取未知多頻微弱信號(hào)各信號(hào)分量的頻率參數(shù)，具有實(shí)用性和普遍性。仿真結(jié)果表明：該方法為微弱多頻信號(hào)處理提供了一種解決方案，拓寬了Duffing振子在微弱信號(hào)檢測(cè)領(lǐng)域的應(yīng)用。