王夢(mèng)蛟 周澤權(quán) 李志軍 曾以成

1)(湘潭大學(xué)信息工程學(xué)院,湘潭 411105)

2)(湘潭大學(xué)物理與光電工程學(xué)院,湘潭 411105)

1 引 言

混沌是一種類隨機(jī)的無規(guī)則運(yùn)動(dòng),是指在確定性非線性系統(tǒng)中不外加隨機(jī)因素所產(chǎn)生的內(nèi)秉隨機(jī)行為[1].已有研究發(fā)現(xiàn),在電子、氣象、水文以及通信等領(lǐng)域中普遍存在混沌現(xiàn)象[2].混沌理論已被廣泛應(yīng)用于保密通信[3]、微弱信號(hào)檢測(cè)[4]和圖像加密[5]等領(lǐng)域.由于實(shí)測(cè)混沌信號(hào)通常受到噪聲的污染,使得刻畫混沌信號(hào)特征的不變系統(tǒng)參數(shù)的計(jì)算變得十分困難[6].有效抑制混沌信號(hào)中的噪聲是對(duì)混沌信號(hào)進(jìn)行有效分析和處理的前提.傳統(tǒng)的線性濾波和譜分析方法無法實(shí)現(xiàn)對(duì)混沌信號(hào)中的噪聲進(jìn)行有效抑制[7,8],因此,針對(duì)混沌信號(hào)的特征研究相應(yīng)的噪聲抑制方法具有重要意義.

目前,混沌信號(hào)的去噪方法主要有以下五類:局部投影(local projection,LP)方法[9,10],該方法基于相空間重構(gòu)的思想,采用LP策略對(duì)噪聲進(jìn)行估計(jì),這類方法在高信噪比的情況下能取得很好的效果,但鄰域半徑的選取會(huì)隨著信噪比的降低而增大,從而導(dǎo)致去噪性能的降低;小波閾值(wavelet thresholding,WT)方法[11,12],該方法利用含噪信號(hào)經(jīng)小波變換后信號(hào)主要成分的小波系數(shù)顯著大于噪聲成分的小波系數(shù),通過閾值策略對(duì)小波系數(shù)進(jìn)行處理,從而將信號(hào)和噪聲分離,但該方法的去噪性能受小波基和分層數(shù)的影響很大,這降低了該方法的自適應(yīng)性;經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)方法[13,14],該方法利用EMD的數(shù)據(jù)驅(qū)動(dòng)自適應(yīng)分解特性克服了WT方法中必須先選定小波基和分解層數(shù)的問題,但這類方法存在閾值難以確定的問題;局部最小二乘多項(xiàng)式擬合方法[15,16],該方法先將含噪信號(hào)進(jìn)行分段平滑預(yù)處理,然后將預(yù)處理后的每段信號(hào)進(jìn)行最小二乘多項(xiàng)式擬合,從而實(shí)現(xiàn)信號(hào)去噪,這類方法本質(zhì)上是一種自適應(yīng)低通濾波方法,信噪比提升性能有限;協(xié)同濾波方法[17],該方法受塊匹配三維濾波的啟發(fā)[18,19],基于混沌信號(hào)具有無限嵌套形式的自相似結(jié)構(gòu),通過對(duì)觀測(cè)信號(hào)進(jìn)行分組、協(xié)同濾波和聚合處理實(shí)現(xiàn)去噪.相對(duì)于其他方法,協(xié)同濾波方法充分利用了混沌信號(hào)的自相似結(jié)構(gòu)特征,具有更高的信噪比提升性能.協(xié)同濾波方法的主要濾波參數(shù)有塊寬、搜索窗長、搜索窗移動(dòng)步長和分組塊數(shù)量,文獻(xiàn)[17]對(duì)濾波參數(shù)的選取進(jìn)行了分析,但只是依據(jù)實(shí)驗(yàn)分析給出了濾波參數(shù)的取值范圍,而沒有提出濾波參數(shù)自動(dòng)優(yōu)化選取的方法,這降低了該方法的自適應(yīng)性.

近年來,各種熵算法在時(shí)間序列分析中得到了廣泛的應(yīng)用[20?22],其中排列熵算法因其較好的魯棒性和快速簡捷的優(yōu)點(diǎn)受到了廣泛的關(guān)注.本文針對(duì)文獻(xiàn)[17]中混沌信號(hào)協(xié)同濾波去噪算法的參數(shù)優(yōu)化問題,利用排列熵分析重構(gòu)混沌信號(hào)的復(fù)雜度,通過計(jì)算不同濾波參數(shù)情況下重構(gòu)信號(hào)的排列熵確定最優(yōu)濾波參數(shù),實(shí)現(xiàn)濾波參數(shù)的自動(dòng)最優(yōu)化,達(dá)到最優(yōu)濾波效果;分析濾波參數(shù)對(duì)去噪效果的影響及其原因;通過仿真實(shí)驗(yàn)對(duì)濾波參數(shù)自動(dòng)最優(yōu)化準(zhǔn)則的可行性進(jìn)行分析.

2 去噪算法原理

文獻(xiàn)[17]首次將協(xié)同濾波應(yīng)用于混沌信號(hào)去噪,提出混沌信號(hào)協(xié)同濾波去噪算法.該算法分為相似塊分組、協(xié)同濾波和聚合重構(gòu)三個(gè)步驟,其去噪流程如圖1所示.

圖1 混沌信號(hào)協(xié)同濾波去噪流程圖Fig.1.Flowchart of the denoising algorithm for chaotic signals based on collaborative filtering.

2.1 相似塊分組

記參考?jí)K為R,任意塊S與R的相似度采用兩者的歸一化距離來度量:

其中,符號(hào)‖·‖2表示求2范數(shù);w為塊寬;d為相似塊之間的距離,該距離越小則塊S與R的相似度越高.設(shè)搜索窗長為l(l?w),在搜索窗內(nèi)以參考?jí)KR為中心搜索與其距離最小的m個(gè)塊構(gòu)成分組group(R),將其以m×w的二維數(shù)組形式保存,即group(R)∈Rm×w.以步長δ將參考?jí)K從觀測(cè)信號(hào)的始端向末端移動(dòng),搜索窗也隨之移動(dòng),對(duì)在每個(gè)位置獲得的相似組進(jìn)行記錄,并將分組中每個(gè)塊的位置進(jìn)行標(biāo)記.

2.2 協(xié)同濾波

協(xié)同濾波分為信號(hào)變換、閾值處理和信號(hào)逆變換三步.

首先,將分組進(jìn)行二維變換:

然后,對(duì)變換系數(shù)進(jìn)行閾值處理,將小于或等于閾值的系數(shù)置零而將大于閾值的系數(shù)保留,閾值操作用符號(hào)HT(·)表示定義為其中,λ(R)為分組group(R)的閾值,其值由文獻(xiàn)[23]提出的VisuShrink方法確定.VisuShrink閾值定義為

其中,σ為觀測(cè)信號(hào)所含噪聲的標(biāo)準(zhǔn)差,其值通過求系數(shù)矩陣G(R)的中位數(shù)絕對(duì)偏差進(jìn)行估計(jì):

最后,通過逆變換得到濾波后的分組:

2.3 聚合重構(gòu)

由于相似塊之間存在重疊,一個(gè)信號(hào)點(diǎn)通常會(huì)同時(shí)屬于多個(gè)相似塊,信號(hào)點(diǎn)的重構(gòu)通過聚合所有包含該點(diǎn)的相似塊在該位置的濾波輸出實(shí)現(xiàn),聚合方式采用算術(shù)平均:

其中,xS,R(n)為分組group(R)中相似塊S在信號(hào)點(diǎn)n的濾波輸出,

3 參數(shù)優(yōu)化準(zhǔn)則

混沌信號(hào)協(xié)同濾波去噪算法充分利用了混沌信號(hào)的自相似結(jié)構(gòu)特征,將非局部均值方法與變換域方法成功結(jié)合,具有良好的信噪比提升性能.該算法的主要濾波參數(shù)為塊寬、搜索窗長、搜索窗移動(dòng)步長和分組塊數(shù)量,濾波參數(shù)的選取決定了濾波的效果.經(jīng)實(shí)驗(yàn)分析發(fā)現(xiàn)濾波參數(shù)的選取受到信號(hào)特征、采樣頻率和噪聲水平的影響,因此在實(shí)際應(yīng)用中最優(yōu)濾波參數(shù)是不斷變化的,怎樣實(shí)現(xiàn)濾波參數(shù)的自動(dòng)最優(yōu)化決定了算法的自適應(yīng)性.

文獻(xiàn)[17]中分析了混沌信號(hào)協(xié)同濾波去噪算法主要濾波參數(shù)的選取,給出了各濾波參數(shù)的取值范圍,然而經(jīng)實(shí)驗(yàn)分析發(fā)現(xiàn)在不同的信號(hào)特征、采樣頻率和噪聲水平情況下該取值范圍并不具有普適性,某些情況下甚至?xí)顾惴ǖ男旁氡忍嵘阅艿陀谝延兴惴?這降低了該算法的自適應(yīng)性.相對(duì)于其他熵算法,排列熵算法能以更快的計(jì)算速度和更高的準(zhǔn)確度刻畫混沌時(shí)間序列的復(fù)雜度,且具有較好的魯棒性[20,24].為了進(jìn)一步完善混沌信號(hào)協(xié)同濾波去噪算法,提高其自適應(yīng)性,本文基于排列熵算法提出一種混沌信號(hào)協(xié)同濾波參數(shù)自動(dòng)優(yōu)化準(zhǔn)則.

3.1 排列熵算法

設(shè)一長度為N的時(shí)間序列{x(t),t= 1,2,···,N},對(duì)其進(jìn)行相空間重構(gòu)得到

其中,d為嵌入維數(shù),L為延遲時(shí)間,將重構(gòu)向量Xt中的d個(gè)分量按升序重新排列:

當(dāng)出現(xiàn)x(t+(ji1?1)L)=x(t+(ji2?1)L)的情況時(shí),按j值的大小排列,即對(duì)ji1＜ji2有x(t+(ji1?1)L)≤x(t+(ji2?1)L).因此,每個(gè)重構(gòu)向量Xt可得到一組符號(hào)序列

其中h=1,2,···,k,由于d個(gè)不同的符號(hào)共有d!種排列可能,所以k≤d!.統(tǒng)計(jì)每種序列出現(xiàn)的次數(shù)并計(jì)算其概率,記為P1,P2,···,Pk,則依據(jù)Shannon信息熵定義時(shí)間序列{x(t),t=1,2,···,N}的k種不同符號(hào)序列的排列熵為[25]

當(dāng)Pz=1/d!時(shí),Hp(d)達(dá)到最大值ln(d!).為便于分析,利用ln(d!)對(duì)排列熵進(jìn)行歸一化處理,即

嵌入維數(shù)d和延遲時(shí)間L的選取影響排列熵算法的有效性,參考文獻(xiàn)[20,24]中對(duì)d和L選取的討論,并綜合實(shí)驗(yàn)分析得出當(dāng)d=5,L=1時(shí)排列熵能有效刻畫重構(gòu)混沌信號(hào)的復(fù)雜度,本文后續(xù)部分計(jì)算排列熵均取d=5,L=1.

3.2 參數(shù)自適應(yīng)選取

混沌信號(hào)的復(fù)雜度能通過排列熵進(jìn)行有效的刻畫,記排列熵值為PE,PE越大表示混沌信號(hào)具有越高的復(fù)雜度,反之則表示混沌信號(hào)的復(fù)雜度越低.當(dāng)混沌信號(hào)受到噪聲污染時(shí)其復(fù)雜度會(huì)升高,且噪聲水平越高復(fù)雜度越高.下面對(duì)典型混沌信號(hào)及其在不同噪聲水平情況下的排列熵進(jìn)行分析.

Lorenz系統(tǒng)方程為[26]

當(dāng)參數(shù)取典型值α=10,γ=28,b=8/3時(shí)系統(tǒng)處于混沌狀態(tài).采用四階龍格-庫塔方法對(duì)方程進(jìn)行求解,初值取[0.5,0.5,0.5],步長取0.01,取系統(tǒng)狀態(tài)變量x產(chǎn)生Lorenz混沌信號(hào),并對(duì)信號(hào)加上零均值高斯白噪聲,產(chǎn)生信噪比SNR為15,20,25,30,35 dB的含噪Lorenz混沌信號(hào).純凈信號(hào)和不同噪聲水平信號(hào)的PE如表1所列,隨著噪聲水平的增加,含噪混沌信號(hào)的復(fù)雜度也相應(yīng)地增加.

表1 不同噪聲水平Lorenz混沌信號(hào)的排列熵值Table 1.PE of the Lorenz’s chaotic signal at different noise levels.

Chen系統(tǒng)方程為[27]

當(dāng)參數(shù)取典型值a=35,b=3,c=28時(shí)系統(tǒng)處于混沌狀態(tài).采用四階龍格-庫塔方法對(duì)方程進(jìn)行求解,初值取[0.5,0.5,0.5],步長取0.01,取系統(tǒng)狀態(tài)變量x產(chǎn)生Chen混沌信號(hào),并對(duì)信號(hào)加上零均值高斯白噪聲,產(chǎn)生SNR為15,20,25,30,35 dB的含噪Chen混沌信號(hào).純凈信號(hào)和不同噪聲水平信號(hào)的PE如表2所列,隨著噪聲水平的增加,含噪混沌信號(hào)的復(fù)雜度也相應(yīng)地增加.

表2 不同噪聲水平Chen混沌信號(hào)的排列熵值Table 2.PE of the Chen’s chaotic signal at different noise levels.

文獻(xiàn)[17]給出的濾波參數(shù)選取范圍為:搜索窗長l取4000左右,搜索窗移動(dòng)步長δ取塊寬w的1/2到1/4左右,分組塊數(shù)量m取30,塊寬w的取值不應(yīng)小于100.然而,經(jīng)實(shí)驗(yàn)分析發(fā)現(xiàn),在不同的信號(hào)特征、采樣頻率和噪聲水平情況下,對(duì)濾波參數(shù)搜索窗長、搜索窗移動(dòng)步長和分組塊數(shù)量取文獻(xiàn)[17]給出的范圍能取得比較好的濾波效果,但關(guān)鍵濾波參數(shù)塊寬的選取卻不在所給的范圍之內(nèi),重構(gòu)混沌信號(hào)中所含噪聲越少則其復(fù)雜度越低,本文基于排列熵算法的混沌信號(hào)協(xié)同濾波參數(shù)自動(dòng)優(yōu)化準(zhǔn)則為:取搜索窗長l=4000,搜索窗移動(dòng)步長δ=w/4,分組塊數(shù)量m=30,采用不同塊寬w對(duì)含噪混沌信號(hào)進(jìn)行去噪,計(jì)算重構(gòu)混沌信號(hào)的PE,最小PE對(duì)應(yīng)的塊寬即為最優(yōu)濾波塊寬.設(shè)重構(gòu)混沌信號(hào)為(t),則最佳濾波塊寬wopt為

其中,符號(hào)HP[·]表示求排列熵.

4 仿真結(jié)果與分析

信號(hào)特征、采樣頻率和噪聲水平影響濾波參數(shù)的選取,為了驗(yàn)證本文所提方法的有效性,本節(jié)分別將該方法應(yīng)用于不同信號(hào)特征、采樣頻率和噪聲水平情況下濾波參數(shù)的選取,為了方便分析,記濾波前后信號(hào)的信噪比分別為SNRin,SNRout.

4.1 不同信號(hào)特征

經(jīng)實(shí)驗(yàn)分析發(fā)現(xiàn),混沌信號(hào)協(xié)同濾波去噪算法具有低通濾波性質(zhì),濾波參數(shù)塊寬和低通截止頻率成反比關(guān)系.考慮到不同混沌信號(hào)的頻譜寬度不同,因此針對(duì)不同的混沌信號(hào)要選取不同的塊寬,以Lorenz混沌信號(hào)和Chen混沌信號(hào)為例,其頻譜如圖2所示.

由圖2可知,相對(duì)于Lorenz信號(hào),Chen信號(hào)的頻譜更寬.由上述分析可知,將混沌信號(hào)協(xié)同濾波去噪算法應(yīng)用于這兩類典型混沌信號(hào)時(shí),對(duì)Chen信號(hào)應(yīng)選取比Lorenz信號(hào)更窄的塊寬.分別給Lorenz信號(hào)和Chen信號(hào)加上零均值高斯白噪聲,信噪比取10 dB,采樣時(shí)間取0.01,取不同寬度的塊寬對(duì)這兩類含噪混沌信號(hào)進(jìn)行去噪,并計(jì)算不同塊寬對(duì)應(yīng)重構(gòu)信號(hào)的排列熵.如圖3和圖4所示,重構(gòu)信號(hào)的排列熵體現(xiàn)了對(duì)應(yīng)濾波塊寬的濾波效果,在當(dāng)前采樣頻率和噪聲水平情況下,Lorenz信號(hào)和Chen信號(hào)的最優(yōu)濾波塊寬分別為100和60,該結(jié)果與上述分析相符.

圖2 Lorenz混沌信號(hào)和Chen混沌信號(hào)頻譜Fig.2.Frequency spectra of the Lorenz’s chaotic signal and Chen’s chaotic signal.

圖3 含噪Lorenz信號(hào)去噪 (a)不同塊寬對(duì)應(yīng)的PE;(b)不同塊寬的去噪效果Fig.3.Denoising of the noisy Lorenz’s signal:(a)The PE corresponding to different block widths;(b)the SNRoutcorresponding to different block widths.

圖4 含噪Chen信號(hào)去噪 (a)不同塊寬對(duì)應(yīng)的PE;(b)不同塊寬的去噪效果Fig.4.Denoising of the noisy Chen’s signal:(a)The PE corresponding to different block widths;(b)the SNRoutcorresponding to different block widths.

4.2 不同采樣頻率

混沌信號(hào)協(xié)同濾波去噪算法在高采樣頻率下能取得更好的去噪效果,對(duì)于不同的采樣頻率最優(yōu)濾波塊寬的選取也不同.采樣時(shí)間ts分別取0.02,0.01,0.0075,0.005產(chǎn)生Lorenz信號(hào),給其加上零均值高斯白噪聲,信噪比取10 dB,取不同寬度的塊寬分別對(duì)不同采樣頻率的含噪信號(hào)進(jìn)行去噪.如圖5和表3所示,采樣頻率越高去噪效果越好,最佳濾波塊寬也越寬.

表3 不同采樣頻率下的最優(yōu)濾波塊寬Table 3.Optimal block width at different sampling frequency.

圖5 采樣頻率對(duì)塊寬的影響 (a)不同采樣頻率對(duì)應(yīng)的PE;(b)不同采樣頻率的去噪效果Fig.5.Effect of sampling frequency on block width:(a)The PE corresponding to different sampling frequency;(b)the SNRoutcorresponding to different sampling frequency.

4.3 不同噪聲水平

混沌信號(hào)在頻域通常主要集中在低頻段,而白噪聲為均勻分布,由于混沌信號(hào)協(xié)同濾波去噪算法具有低通濾波性質(zhì),且濾波參數(shù)塊寬和低通截止頻率成反比關(guān)系,因此高信噪比情況下的最優(yōu)濾波塊寬比低信噪比情況下的最優(yōu)濾波塊寬稍窄采樣時(shí)間取0.01產(chǎn)生Lorenz信號(hào),給其加上零均值高斯白噪聲,信噪比分別取5,10,15,20 dB,取不同寬度的塊寬分別對(duì)不同噪聲水平的含噪信號(hào)進(jìn)行去噪.如圖6和表4所示,噪聲水平越高對(duì)應(yīng)的最優(yōu)濾波塊寬越寬.

圖6 噪聲水平對(duì)塊寬的影響 (a)不同噪聲水平對(duì)應(yīng)的PE;(b)不噪聲水平的去噪效果Fig.6.Effect of noise level on block width:(a)The PE corresponding to different noise levels;(b)the SNRout corresponding to different noise levels.

表4 不同噪聲水平下的最優(yōu)濾波塊寬Table 4.Optimal block width at different noise levels.

5 結(jié) 論

本文基于排列熵提出了一種混沌信號(hào)協(xié)同濾波參數(shù)自動(dòng)優(yōu)化準(zhǔn)則,深入分析了在不同信號(hào)特征、采樣頻率和噪聲水平情況下混沌信號(hào)協(xié)同濾波中各參數(shù)的選取問題.仿真結(jié)果表明,雖然濾波參數(shù)的選取受到不同信號(hào)特征、采樣頻率和噪聲水平的影響,新的參數(shù)優(yōu)化準(zhǔn)則能有效實(shí)現(xiàn)不同情況下濾波參數(shù)的自動(dòng)最優(yōu)化.該準(zhǔn)則提高了混沌信號(hào)協(xié)同濾波去噪算法的自適應(yīng)性,使其更符合實(shí)際應(yīng)用需求,在各類混沌信號(hào)去噪中具有更廣泛的應(yīng)用價(jià)值.

[1]Lü J H,Lu J A,Chen S H 2002The Analysis and Applications of Chaotic Time Series(Wuhan:Wuhan University Press)pp1–8(in Chinese)[呂金虎,陸君安,陳士華 2002混沌時(shí)間序列分析及其應(yīng)用(武漢:武漢大學(xué)出版社)第1—8頁]

[2]Han M,Xu M L 2013Acta Phys.Sin.62 120510(in Chinese)[韓敏,許美玲 2013物理學(xué)報(bào)62 120510]

[3]Sun J W,Shen Y,Yin Q,Xu C J 2013Chaos23 013140

[4]Li G Z,Zhang B 2017IEEE Trans.Ind.Electron.64 2255

[5]Peng G Y,Min F H 2017Nonlinear Dynam.90 1607

[6]Urbanowicz K,Ho?yst J A 2003Phys.Rev.E67 046218

[7]Feng J C 2012Chaotic Signals and Information Processing(Beijing:Tsinghua University Press)pp32–35(in Chinese)[馮久超 2012混沌信號(hào)與信息處理(北京:清華大學(xué)出版社)第32—35頁]

[8]Badii R,Broggi G,Derighetti B,Ravani M 1988Phys.Rev.Lett.60 979

[9]Cawley R,Hsu G H 1992Phys.Rev.A46 3057

[10]Schreiber T,Richter M 1999Int.J.Bifurcat.Chaos9 2039

[11]Donoho D L 1995IEEE Trans.Inf.Theory41 613

[12]Han M,Liu Y H,Xi J H,Guo W 2007IEEE Signal Process.Lett.14 62

[13]Kopsinis Y,McLaughlin S 2009IEEE Trans.Signal Process.57 1351

[14]Wang W B,Zhang X D,Wang X L 2013Acta Phys.Sin.62 050201(in Chinese)[王文波,張曉東,汪祥莉 2013物理學(xué)報(bào)62 050201]

[15]Tung W W,Gao J B,Hu J,Yang L 2011Phys.Rev.E83 046210

[16]Gao J B,Sultan H,Hu J,Tung W W 2010IEEE Signal Process.Lett.17 237

[17]Chen Y,Liu X Y,Wu Z T,Fan Y,Ren Z L,Feng J C 2017Acta Phys.Sin.66 210501(in Chinese)[陳越, 劉雄英,吳中堂,范藝,任子良,馮久超 2017物理學(xué)報(bào) 66 210501]

[18]Dabov K,Foi A,Katkovnik V,Egiazarian K 2007IEEE Trans.Image Process.16 2080

[19]Yadav S K,Sinha R,Bora P K 2015IET Signal Process.9 88

[20]Hou W,Feng G L,Dong W J,Li J P 2006Acta Phys.Sin.55 2663(in Chinese)[侯威,封國林,董文杰,李建平2006物理學(xué)報(bào)55 2663]

[21]Sun K H,He S B,Yin L Z,A D L·Duo L K 2012Acta Phys.Sin.61 130507(in Chinese)[孫克輝,賀少波,尹林子,阿地力·多力坤2012物理學(xué)報(bào)61 130507]

[22]Yu M Y,Sun K H,Liu W H,He S B 2018Chaos Solitons Fractals106 107

[23]Donoho D L,Johnstone I M 1994Biometrika81 425

[24]He S B,Sun K H,Wang H H 2016Physical A461 812

[25]Bandt C,Pompe B 2002Phys.Rev.Lett.88 174102

[26]Lorenz E N 1963J.Atmos.Sci.20 130

[27]Chen G R,Ueta T 1999Int.J.Bifurcat.Chaos9 1465