黃澤徽 李亞安 陳哲 劉戀

(西北工業(yè)大學(xué)航海學(xué)院, 西安 710072)

1 引 言

混沌系統(tǒng)的初值敏感性和噪聲免疫性使它在微弱信號檢測領(lǐng)域得到了廣泛的關(guān)注, 為微弱信號的檢測提供了新的思路. 美國的Birx 博士在他的博士論文中最先提出將混沌理論應(yīng)用于弱信號的檢測. 王冠宇等[1]在國外學(xué)者研究工作的基礎(chǔ)上,利用Duffing 振子對白噪聲背景下的微弱正弦信號進行了檢測研究, 實現(xiàn)了對頻率已知的待測信號幅值的估計, 檢測信噪比可以達到–26 dB, 推動了信號檢測領(lǐng)域的發(fā)展. 李月等[2?5]研究了強噪聲下混沌振子對微弱正弦信號、方波信號的檢測, 色噪聲背景下混沌振子對微弱正弦信號、方波信號的檢測, 指出混沌振子系統(tǒng)對零均值噪聲具有很強的免疫力. 之后大量學(xué)者對Duffing 檢測系統(tǒng)進行改進,賴智慧等[6]對基于Homles-Duffing 振子的混沌檢測進行了改進, 提出變尺度微弱信號檢測方法; 叢超等[7]提出了一種基于適應(yīng)步長型間歇混沌振子的信號檢測方法, 用一個振子系統(tǒng)可以實現(xiàn)對任意頻率任意相位的微弱周期信號的搜索檢測; 牛德智等[8]針對同頻微弱信號檢測時存在的盲區(qū), 提出了一種策動力移相法予以消除; 陳志光等[9]利用間歇混沌現(xiàn)象對頻率未知信號進行檢測并取得較好的效果; 時培明等[10]提出了一種基于雙耦合Duffing混沌振子與變尺度相結(jié)合的微弱信號檢測新方法,將其用于檢測任意多頻微弱信號具有明顯優(yōu)勢.通過大量學(xué)者的不懈努力, Duffing 檢測不斷趨于成熟.

Duffing 混沌系統(tǒng)提供了與傳統(tǒng)方法不同的檢測途徑, 在實際應(yīng)用中, 混沌系統(tǒng)躍變閾值的確定是利用Duffing 系統(tǒng)準確進行微弱信號檢測的關(guān)鍵, 但是傳統(tǒng)的時序圖方法不利于計算機的自動處理. Lyapunov 指數(shù)法算法復(fù)雜, 且容易受到噪聲的影響. 近年來有大量學(xué)者對閾值求解方法進行了研究, 提出了相圖分割法[11]、二分法[12]、0-1 方法[13,14]、基于龐加萊界面的定量判別方法[15]等來計算混沌系統(tǒng)躍變閾值, 但仍具有一定的局限性. 本文將熵引入到閾值求解, 熵可以用來表征信號復(fù)雜度, 已有不少學(xué)者利用熵對信號復(fù)雜度進行分析, 梁滌青等[16]利用小波包能量熵來判斷混沌序列的復(fù)雜度;楊孝敬等[17]利用模糊近似熵對磁共振信號復(fù)雜度進行分析; 陳祥龍等[18]利用多尺度樣本熵對壓縮機工作時的振動信號進行了分析; 王鴻姍等[19]利用小波包樣本熵通過分析聲信號的復(fù)雜度進行異常音提取. 本文發(fā)現(xiàn)Duffing 系統(tǒng)在不同狀態(tài)下多尺度熵的明顯差別, 利用這一現(xiàn)象通過分析系統(tǒng)多尺度熵與策動力幅值的變化關(guān)系, 提出一種基于多尺度熵的Duffing 混沌系統(tǒng)躍變閾值確定方法, 并通過仿真計算說明其可行性.

2 Duffing 混沌系統(tǒng)微弱信號檢測原理

Duffing 振子檢測系統(tǒng)模型為:

其中k是阻尼比;x?x3是非線性恢復(fù)力;rcos(wt)是周期性策動力,r是策動力幅值,h是目標信號的幅值,w是待測信號的頻率,z為白噪聲.

調(diào)節(jié)策動力幅值使r從0 逐漸增大, 系統(tǒng)將會遍歷同宿軌道狀態(tài)、分叉狀態(tài)、混沌狀態(tài)、臨界混沌狀態(tài), 當r大于某一值rd時, 系統(tǒng)由臨界混沌狀態(tài)進入大尺度周期狀態(tài), 相軌跡圖將發(fā)生明顯的變化, 如圖1 所示, 其中rd便是本文所要求解的混沌系統(tǒng)躍變閾值. 使用Duffing 混沌系統(tǒng)檢測微弱信號正是基于這一躍變現(xiàn)象: 將系統(tǒng)的策動力幅值r調(diào)為rd, 再將待測信號添加到系統(tǒng), 觀察待測信號添加前后系統(tǒng)的相軌跡是否發(fā)生躍變, 如果發(fā)生了躍變則說明待測信號中有目標信號的存在. 可知利用這一現(xiàn)象檢測微弱信號時, 躍變閾值的確定非常關(guān)鍵, 躍變閾值選取過小, 將會嚴重影響到系統(tǒng)的檢測信噪比; 躍變閾值選取過大, 會導(dǎo)致系統(tǒng)狀態(tài)發(fā)生躍變, 造成對檢測結(jié)果的誤判.

現(xiàn)有計算混沌系統(tǒng)躍變閾值的方法存在一定的局限性, 直觀的相軌跡分析方法簡單直觀、判別方便并且不需要復(fù)雜的計算, 但是該方法受到諸多因素的影響, 并且在判別過程中會受到實驗人員主觀因素的影響, 精度低誤差大; Melnikov 和Shilnikov方法可以計算Duffing 系統(tǒng)的混沌范圍, 但是都只能說明當策動力幅值處于這個范圍時, 系統(tǒng)可能處于混沌狀態(tài), 而不能確定系統(tǒng)一定處于混沌狀態(tài);Lyapunov 指數(shù)用來度量動力學(xué)性態(tài)的規(guī)則性程度, 描述了系統(tǒng)軌跡收斂或發(fā)散的比率, 但是精確的Lyapunov 指數(shù)難以得到且容易受到噪聲的影響. 為此本文基于多尺度熵提出一種簡單有效的計算系統(tǒng)躍變閾值的方法.

圖1 系統(tǒng)狀態(tài)躍變 (a) 臨界混沌狀態(tài); (b) 周期狀態(tài)Fig. 1. System state transition: (a) System state transition;(b) periodic state.

3 基于多尺度熵的混沌系統(tǒng)躍變閾值確定方法

熵是系統(tǒng)混亂或無序程度的度量, 系統(tǒng)混亂或者無序的程度越高, 熵值越大, 反之越小. 分析熵值隨策動力幅值變化關(guān)系后發(fā)現(xiàn), 系統(tǒng)由混沌狀態(tài)進入周期狀態(tài)后, 系統(tǒng)熵值明顯變小且趨于平穩(wěn),本文利用該現(xiàn)象基于多尺度熵((multi-scale sample entropy, MsEn)來求解系統(tǒng)躍變閾值.

3.1 多尺度熵算法

多尺度樣本熵(MsEn)用于描述在不同尺度上時間序列的無規(guī)則程度, 包含參數(shù)m,s,τ, 其中m是嵌入維數(shù),s是相似系數(shù),τ是尺度因子. 多尺度樣本熵由兩部分構(gòu)成: 一是對輸入的時間序列按照尺度τ進行粗細化操作得到新的時間序列; 一是計算尺度τ下新時間序列的樣本熵(sample entropy,SampEn). 具體步驟為[20]:

步驟1對原始時間序列按照尺度τ進行粗細化操作

假定有一離散時間序列x1,x2,x3,··· ,xM共有M個點, 在尺度為τ時對時間序列進行粗-斷點變換得到新的時間序列為:

經(jīng)過粗-斷點變換后得到的新時間序列長度為N=M/τ.

步驟2計算尺度τ下新時間序列的樣本熵

1)將新的時間序列從Y(τ)(1)到Y(jié)(τ)(N?m+1)組成一組m維矢量, 第i個矢量Y τ(i) 為:

式中i=1～N ?m+1 .

2)定義d[Y τ(i),Y τ(j)]為矢量Y τ(i) 與矢量Y τ(j) 對應(yīng)元素中差值最大的一個, 即:

3)給出相似系數(shù)s, 計算的數(shù)目Di, 及該數(shù)目與總的矢量數(shù)目N ?m的比值, 記作

4)將嵌入維數(shù)變成m+1 , 重復(fù)前面的過程,得到和

5)時間序列為有限長時, 得到在尺度τ下的樣本熵為

對原始時間序列計算每個尺度τ所對應(yīng)的樣本熵值, 得到多尺度樣本熵的計算公式為

3.2 閾值求解

研究中發(fā)現(xiàn)系統(tǒng)混沌態(tài)和周期態(tài)輸出的時間序列多尺度熵值明顯不同, 且當系統(tǒng)進入周期態(tài)后, 系統(tǒng)輸出序列的多尺度熵趨于平穩(wěn), 基于此通過描述系統(tǒng)多尺度熵值隨策動力幅值變化情況可以確定系統(tǒng)躍變閾值.

3.2.1 方波信號檢測系統(tǒng)躍變閾值求解

求解方波信號檢測系統(tǒng)的躍變閾值, 系統(tǒng)所采用檢測方程為:

我想我當時的樣子一定很傻，一定呆得像只木雞，但你知道的，在那個愛做夢，對漂亮女孩想入非非的年齡，一般人應(yīng)該都是這種反應(yīng)的。

取w= 1 rad/s,k= 0.5, 初值x0= (1, 1), 時間序列長度N= 30000, 尺度因子t= 10, 嵌入維數(shù)m= 2, 相似容限s= 0.1SDx(時間序列標準差), 分析多尺度熵值隨策動力幅值r變化情況如圖2 所示. 對圖2 進行分析可知, 系統(tǒng)策動力幅值r?0.605 , 多尺度熵值較大且會隨著策動力幅值r波動,r >0.605 后多尺度熵值趨于平穩(wěn)且較之前明顯變小, 根據(jù)熵的含義及系統(tǒng)變化規(guī)律可以判定r >0.605 后系統(tǒng)處于周期狀態(tài), 所以確定該系統(tǒng)的躍變閾值為rd=0.605. 注意在r=0.603 時, 系統(tǒng)多尺度熵值為極小值, 如果只以熵值的大小作為評判標準, 忽略平穩(wěn)性的評判, 會認為此時系統(tǒng)已處于周期狀態(tài), 但從圖中可以看到此時系統(tǒng)熵值并未趨于平穩(wěn), 在之后系統(tǒng)的熵值仍在波動, 所以僅根據(jù)熵值大小來求解閾值并不準確. 對該系統(tǒng)進行仿真實驗, 通過分析系統(tǒng)相軌跡圖對比該方法計算結(jié)果的準確性, 結(jié)果如圖3 所示, 策動力幅值r =0.605 時系統(tǒng)處于左圖所示的混沌狀態(tài), 幅值增加到0.606 時系統(tǒng)便躍變至右圖所示的周期狀態(tài), 通過仿真實驗確定的系統(tǒng)躍變閾值為0.605,多尺度熵方法所求結(jié)果與其一致. 可見, 多尺度熵方法可以準確確定方波信號檢測系統(tǒng)躍變閾值.

圖2 方波信號檢測系統(tǒng)多尺度熵變化情況Fig. 2. Variation of multi-scale entropy in square wave signal detection system.

圖3 仿真實驗求解方波信號檢測系統(tǒng)閾值Fig. 3. Simulation experiment to solve the threshold of square wave signal detection system.

3.2.2 正弦信號檢測系統(tǒng)躍變閾值求解

求解正弦信號檢測系統(tǒng)的躍變閾值, 利用Duffing 系統(tǒng)(1)式, 取w = 1 rad/s, k=0.5 , 初值x0= (1, 1), 對于多尺度熵算法取時間序列長度N = 30000, 尺度因子 τ=10, 嵌入維數(shù) m=2 ,相似容限 s=0.1SDx, 分析該系統(tǒng)多尺度熵值隨策動力幅值r 變化情況, 結(jié)果如圖4 所示, 可知該系統(tǒng)的躍變閾值為 rd=0.826 . 通過圖5 所示的仿真實驗可以確定該系統(tǒng)的閾值為0.826, 與多尺度熵方法所求結(jié)果一致.

圖4 正弦信號檢測系統(tǒng)多尺度熵變化情況Fig. 4. Variation of multi-scale entropy in sinusoidal signal detection system.

圖5 仿真實驗求解正弦信號檢測系統(tǒng)閾值Fig. 5. Simulation experiment to solve the threshold of sinusoidal signal detection system.

3.2.3 真實信號檢測系統(tǒng)躍變閾值求解

為了驗證多尺度熵方法對真實信號檢測系統(tǒng)的閾值求解效果, 選取一組真實的艦船信號作為樣本數(shù)據(jù), 真實待檢測信號的波形如圖6 所示, 頻域波形如圖7 所示, 分析可知真實水聲信號包含頻率為50.27 Hz 的正弦信號. 在系統(tǒng)中內(nèi)置頻率為50.27 Hz 的正弦信號, 求得該系統(tǒng)閾值為rd=0.825 ,結(jié)果如圖8 所示, 進行仿真實驗也得到rd=0.825的閾值, 如圖9 所示.

圖6 真實水聲信號Fig. 6. Real underwater acoustic signal.

圖7 真實水聲信號頻譜Fig. 7. Spectrum of real underwater acoustic signals.

圖8 真實信號檢測系統(tǒng)多尺度熵值變化情況Fig. 8. Changes in multi-scale entropy of real signal detection system.

將策動力幅值r調(diào)節(jié)至0.825, 也即使系統(tǒng)處于臨界混沌狀態(tài), 將真實水聲信號加入檢測系統(tǒng),系統(tǒng)相軌跡圖由圖10(a)所示混沌狀態(tài)躍變至圖10(b)所示周期狀態(tài), 成功實現(xiàn)對實測水聲信號中目標信號的檢測. 所以多尺度熵方法可以很準確地計算系統(tǒng)閾值.

圖9 仿真實驗求解真實水聲信號檢測系統(tǒng)閾值Fig. 9. Simulation experiment to solve the threshold of real underwater acoustic signal detection system.

圖10 對真實水聲信號的檢測 (a) 系統(tǒng)未添加真實信號; (b)系統(tǒng)添加真實信號Fig. 10. Detecting real underwater acoustic signals: (a) The system did not add a real signal; (b) system adds real signal.

4 最大多尺度熵算法

圖11 5rad/s 正弦信號檢測系統(tǒng)多尺度熵變化情況Fig. 11. Variation of multi-scale entropy in 5 rad/s sinusoidal signal detection system.

圖12 仿真實驗求解5 rad/s 正弦信號檢測系統(tǒng)躍變閾值Fig. 12. Simulation experiment to solve the threshold of 5 rad/s sinusoidal signal detection system.

利用系統(tǒng)(1)將頻率參數(shù)改為5 rad/s, 其他參數(shù)不變, 得到多尺度熵變化情況如圖11 所示,分析可知多尺度熵方法求得系統(tǒng)閾值為rd=0.825 .進行仿真實驗, 仿真結(jié)果如圖12 所示, 可知仿真實驗求得的系統(tǒng)閾值為rd=0.826 , 與多尺度熵方法求解結(jié)果有偏差, 說明多尺度熵方法存在一定問題.

我們分析是時間序列段的隨機選取造成了閾值計算的偏差, 第三部分中計算采用的時間序列是在整個Duffing 序列中隨機選取的長度為30000的子序列段, 在一般狀態(tài)下隨機選取的子序列段可以代表整個時間序列的復(fù)雜度, 但是當系統(tǒng)在接近或處于臨界混沌狀態(tài)時會有部分序列段已經(jīng)處于周期狀態(tài), 如圖13 所示, 該狀態(tài)下隨機選取的時間序列段不能代表整個時間序列的狀態(tài).

圖13 系統(tǒng)臨界混沌狀態(tài)Fig. 13. Critical chaotic state of the system.

圖13 表示系統(tǒng)正處于混沌狀態(tài), 但右圖中a段序列為周期序列, 取用a序列進行多尺度熵計算會得到較小熵值, 誤認為系統(tǒng)已經(jīng)進入周期狀態(tài), 導(dǎo)致所求系統(tǒng)躍變閾值偏小, 這便是在前述計算中多尺度熵方法計算結(jié)果與仿真實驗計算結(jié)果存在偏差且偏小的原因. 針對這一現(xiàn)象考慮在整個序列中尋找復(fù)雜度最大的子序列及其對應(yīng)的多尺度熵值, 用復(fù)雜度最大的子序列代表整個時間序列, 當系統(tǒng)復(fù)雜度最大的子序列對應(yīng)的熵值處于較小值且趨于平穩(wěn)時, 則系統(tǒng)必定已經(jīng)進入周期狀態(tài), 可更加精確計算系統(tǒng)閾值. 為此引入遺傳算法,利用遺傳算法尋找Duffing 序列最大多尺度熵.

4.1 遺傳算法

遺傳算法(genetic algorithm, GA)是基于生物的自然選擇和遺傳機理而形成的一種全局尋優(yōu)算法, 其本質(zhì)是一種基于概率的隨機搜索算法. 一般認為GA 的計算流程為[20]:

1) 問題解的遺傳表示;

2) 產(chǎn)生初始染色體;

3) 設(shè)計適應(yīng)度函數(shù), 根據(jù)適應(yīng)值對個體進行優(yōu)劣判斷;

4) 進行選擇操作, 選出適應(yīng)值高的染色體, 使它們成為新一代種群中的染色體;

5) 對新種群進行交叉操作, 產(chǎn)生新的染色體;

6) 進行變異操作, 避免算法陷入局部最優(yōu)解的情況;

7) 對新的種群重復(fù)進行選擇、交叉、和變異操作;

8) 經(jīng)過給定次數(shù)的迭代或滿足給定的條件后,把最好的染色體作為優(yōu)化問題的最優(yōu)解.

4.2 最大多尺度熵算法步驟

基于遺傳算法的思想, 本文提出一種最大多尺度熵算法, 算法具體步驟為:

1) 利用Duffing 系統(tǒng)方程迭代產(chǎn)生一組時間序列X=[x1,x2,x3,··· ,xm] , 將X作為一個種群,時間序列的長度為m;

2) 將X分成l個染色體, 分別為CC=[C1,C2,C3,··· ,Cl] , 每個染色體的長度為n, 染色體的數(shù)目l=m/n.CC中前i個染色體構(gòu)成的種群C=[C1,C2,C3,··· ,Ci] 作為算法的初始種群, 循環(huán)次數(shù)為N=l ?i;

3) 選取多尺度樣本熵作為適應(yīng)度函數(shù), 通過計算種群中每個染色體的多尺度樣本熵值來判斷各個染色體的優(yōu)劣, 作為選擇的依據(jù);

4) 對初始種群C=[C1,C2,C3,··· ,Ci] 進行適應(yīng)度測試, 計算每個染色體的適應(yīng)度值, 得到種群中每個染色體所對應(yīng)的適應(yīng)度值分別為M=[M1,M2,M3,··· ,Mi] ;

5) 按照得到的適應(yīng)度值, 從大到小對染色體進行排序, 重新排序后得到的種群為對應(yīng)的適應(yīng)度值為. 對重新排序后的種群C′按照隨機數(shù)w(0

6) 得到由i個新的染色體組成的種群, 返回到步驟(4)進行循環(huán)迭代計算, 直到滿足循環(huán)條件;

7) 得到Duffing 時間序列中熵值最大的子序列及其對應(yīng)的熵值.

4.2.1 最大多尺度熵算法有效性

第三部分5 rad/s 正弦信號檢測系統(tǒng)中, 策動力幅值r=0.826 , 系統(tǒng)處于混沌狀態(tài), 多尺度熵方法求得MsEn = 0.0778, 利用最大多尺度熵算法計算r=0.826 時Duffing 序列最大多尺度熵值,取Duffing 序列長度為400000, 每個子序列長度為30000, 初始染色體個數(shù)為6, 交叉點為隨機數(shù), 每次循環(huán)淘汰的染色體數(shù)目為1, 循環(huán)結(jié)束后得到6 個復(fù)雜度較大的子序列段, 其對應(yīng)的熵值結(jié)果如表1 所列.

表1 Duffing 子序列段熵值( r=0.826 )Table 1. Entropy value of the Duffing subsequence segment ( r=0.826 ).

分析表1 可知, 程序經(jīng)過遺傳進化找到了4 組熵值大于0.0778 的混沌子序列, 最大多尺度熵值為0.1188.r=0.827 時, 系統(tǒng)處于周期狀態(tài), 多尺度熵方法求得MsEn=0.0754, 利用最大多尺度熵算法計算r=0.827 時Duffing 子序列熵值, 結(jié)果如表2 所列.

表2 Duffing 子序列段MsEn 值( r=0.827 )Table 2. Entropy value of the Duffing subsequence segment( r=0.827 ).

從表2 中看到, 經(jīng)過遺傳進化找到了5 組大于0.0754 的子序列, 最大MsEn 為0.0822, 但因r= 0.827 時系統(tǒng)已經(jīng)處于周期狀態(tài), 所求結(jié)果與多尺度熵方法所求結(jié)果相差不大. 結(jié)合表1 和表2可以知道, 最大多尺度熵算法可以很有效地尋找到復(fù)雜度較大的子序列, 得到最大多尺度熵, 可更加精確地求解系統(tǒng)的躍變閾值.

針對5 rad/s 正弦信號檢測系統(tǒng)的偏差問題,現(xiàn)用最大多尺度熵方法求該系統(tǒng)閾值, 分析偏差是否仍舊存在. 取Duffing 序列長度為400000, 每個子序列長度為30000, 初始染色體數(shù)目為6, 交叉點為隨機數(shù), 結(jié)果如圖14 所示, 分析可以確定系統(tǒng)躍變閾值為0.826, 與圖12 仿真實驗所求閾值0.826 相同, 所以最大多尺度熵方法可以很好地解決序列段的選取問題, 準確求解系統(tǒng)躍變閾值.

圖145 rad/s 正弦信號檢測系統(tǒng)最大多尺度熵變化情況Fig. 14. Variation of maximum multi-scale entropy of 5 rad/s sinusoidal signal detection system.

5 結(jié) 論

本文針對Duffing 系統(tǒng)躍變閾值難以確定這一問題進行了研究, 首次提出利用Duffing 系統(tǒng)周期態(tài)和混沌態(tài)多尺度熵值差異明顯這一現(xiàn)象, 通過分析多尺度熵與策動力幅值的關(guān)系對閾值進行求解, 對于該方法存在的問題, 結(jié)合遺傳算法進行了改進, 通過對正弦信號與方波信號檢測系統(tǒng)的閾值求解, 證明該方法可以快速準確得到系統(tǒng)躍變閾值, 解決了閾值難以快速準確得到的問題, 為混沌振子檢測低信噪比信號的實際應(yīng)用奠定了很好的基礎(chǔ).