王夢蛟 周澤權 李志軍 曾以成
1)(湘潭大學信息工程學院,湘潭 411105)
2)(湘潭大學物理與光電工程學院,湘潭 411105)
混沌是一種類隨機的無規(guī)則運動,是指在確定性非線性系統(tǒng)中不外加隨機因素所產生的內秉隨機行為[1].已有研究發(fā)現,在電子、氣象、水文以及通信等領域中普遍存在混沌現象[2].混沌理論已被廣泛應用于保密通信[3]、微弱信號檢測[4]和圖像加密[5]等領域.由于實測混沌信號通常受到噪聲的污染,使得刻畫混沌信號特征的不變系統(tǒng)參數的計算變得十分困難[6].有效抑制混沌信號中的噪聲是對混沌信號進行有效分析和處理的前提.傳統(tǒng)的線性濾波和譜分析方法無法實現對混沌信號中的噪聲進行有效抑制[7,8],因此,針對混沌信號的特征研究相應的噪聲抑制方法具有重要意義.
目前,混沌信號的去噪方法主要有以下五類:局部投影(local projection,LP)方法[9,10],該方法基于相空間重構的思想,采用LP策略對噪聲進行估計,這類方法在高信噪比的情況下能取得很好的效果,但鄰域半徑的選取會隨著信噪比的降低而增大,從而導致去噪性能的降低;小波閾值(wavelet thresholding,WT)方法[11,12],該方法利用含噪信號經小波變換后信號主要成分的小波系數顯著大于噪聲成分的小波系數,通過閾值策略對小波系數進行處理,從而將信號和噪聲分離,但該方法的去噪性能受小波基和分層數的影響很大,這降低了該方法的自適應性;經驗模態(tài)分解(empirical mode decomposition,EMD)方法[13,14],該方法利用EMD的數據驅動自適應分解特性克服了WT方法中必須先選定小波基和分解層數的問題,但這類方法存在閾值難以確定的問題;局部最小二乘多項式擬合方法[15,16],該方法先將含噪信號進行分段平滑預處理,然后將預處理后的每段信號進行最小二乘多項式擬合,從而實現信號去噪,這類方法本質上是一種自適應低通濾波方法,信噪比提升性能有限;協同濾波方法[17],該方法受塊匹配三維濾波的啟發(fā)[18,19],基于混沌信號具有無限嵌套形式的自相似結構,通過對觀測信號進行分組、協同濾波和聚合處理實現去噪.相對于其他方法,協同濾波方法充分利用了混沌信號的自相似結構特征,具有更高的信噪比提升性能.協同濾波方法的主要濾波參數有塊寬、搜索窗長、搜索窗移動步長和分組塊數量,文獻[17]對濾波參數的選取進行了分析,但只是依據實驗分析給出了濾波參數的取值范圍,而沒有提出濾波參數自動優(yōu)化選取的方法,這降低了該方法的自適應性.
近年來,各種熵算法在時間序列分析中得到了廣泛的應用[20?22],其中排列熵算法因其較好的魯棒性和快速簡捷的優(yōu)點受到了廣泛的關注.本文針對文獻[17]中混沌信號協同濾波去噪算法的參數優(yōu)化問題,利用排列熵分析重構混沌信號的復雜度,通過計算不同濾波參數情況下重構信號的排列熵確定最優(yōu)濾波參數,實現濾波參數的自動最優(yōu)化,達到最優(yōu)濾波效果;分析濾波參數對去噪效果的影響及其原因;通過仿真實驗對濾波參數自動最優(yōu)化準則的可行性進行分析.
文獻[17]首次將協同濾波應用于混沌信號去噪,提出混沌信號協同濾波去噪算法.該算法分為相似塊分組、協同濾波和聚合重構三個步驟,其去噪流程如圖1所示.
圖1 混沌信號協同濾波去噪流程圖Fig.1.Flowchart of the denoising algorithm for chaotic signals based on collaborative filtering.
記參考塊為R,任意塊S與R的相似度采用兩者的歸一化距離來度量:
其中,符號‖·‖2表示求2范數;w為塊寬;d為相似塊之間的距離,該距離越小則塊S與R的相似度越高.設搜索窗長為l(l?w),在搜索窗內以參考塊R為中心搜索與其距離最小的m個塊構成分組group(R),將其以m×w的二維數組形式保存,即group(R)∈Rm×w.以步長δ將參考塊從觀測信號的始端向末端移動,搜索窗也隨之移動,對在每個位置獲得的相似組進行記錄,并將分組中每個塊的位置進行標記.
協同濾波分為信號變換、閾值處理和信號逆變換三步.
首先,將分組進行二維變換:
然后,對變換系數進行閾值處理,將小于或等于閾值的系數置零而將大于閾值的系數保留,閾值操作用符號HT(·)表示定義為其中,λ(R)為分組group(R)的閾值,其值由文獻[23]提出的VisuShrink方法確定.VisuShrink閾值定義為
其中,σ為觀測信號所含噪聲的標準差,其值通過求系數矩陣G(R)的中位數絕對偏差進行估計:
最后,通過逆變換得到濾波后的分組:
由于相似塊之間存在重疊,一個信號點通常會同時屬于多個相似塊,信號點的重構通過聚合所有包含該點的相似塊在該位置的濾波輸出實現,聚合方式采用算術平均:
其中,xS,R(n)為分組group(R)中相似塊S在信號點n的濾波輸出,
混沌信號協同濾波去噪算法充分利用了混沌信號的自相似結構特征,將非局部均值方法與變換域方法成功結合,具有良好的信噪比提升性能.該算法的主要濾波參數為塊寬、搜索窗長、搜索窗移動步長和分組塊數量,濾波參數的選取決定了濾波的效果.經實驗分析發(fā)現濾波參數的選取受到信號特征、采樣頻率和噪聲水平的影響,因此在實際應用中最優(yōu)濾波參數是不斷變化的,怎樣實現濾波參數的自動最優(yōu)化決定了算法的自適應性.
文獻[17]中分析了混沌信號協同濾波去噪算法主要濾波參數的選取,給出了各濾波參數的取值范圍,然而經實驗分析發(fā)現在不同的信號特征、采樣頻率和噪聲水平情況下該取值范圍并不具有普適性,某些情況下甚至會使算法的信噪比提升性能低于已有算法,這降低了該算法的自適應性.相對于其他熵算法,排列熵算法能以更快的計算速度和更高的準確度刻畫混沌時間序列的復雜度,且具有較好的魯棒性[20,24].為了進一步完善混沌信號協同濾波去噪算法,提高其自適應性,本文基于排列熵算法提出一種混沌信號協同濾波參數自動優(yōu)化準則.
設一長度為N的時間序列{x(t),t= 1,2,···,N},對其進行相空間重構得到
其中,d為嵌入維數,L為延遲時間,將重構向量Xt中的d個分量按升序重新排列:
當出現x(t+(ji1?1)L)=x(t+(ji2?1)L)的情況時,按j值的大小排列,即對ji1<ji2有x(t+(ji1?1)L)≤x(t+(ji2?1)L).因此,每個重構向量Xt可得到一組符號序列
其中h=1,2,···,k,由于d個不同的符號共有d!種排列可能,所以k≤d!.統(tǒng)計每種序列出現的次數并計算其概率,記為P1,P2,···,Pk,則依據Shannon信息熵定義時間序列{x(t),t=1,2,···,N}的k種不同符號序列的排列熵為[25]
當Pz=1/d!時,Hp(d)達到最大值ln(d!).為便于分析,利用ln(d!)對排列熵進行歸一化處理,即
嵌入維數d和延遲時間L的選取影響排列熵算法的有效性,參考文獻[20,24]中對d和L選取的討論,并綜合實驗分析得出當d=5,L=1時排列熵能有效刻畫重構混沌信號的復雜度,本文后續(xù)部分計算排列熵均取d=5,L=1.
混沌信號的復雜度能通過排列熵進行有效的刻畫,記排列熵值為PE,PE越大表示混沌信號具有越高的復雜度,反之則表示混沌信號的復雜度越低.當混沌信號受到噪聲污染時其復雜度會升高,且噪聲水平越高復雜度越高.下面對典型混沌信號及其在不同噪聲水平情況下的排列熵進行分析.
Lorenz系統(tǒng)方程為[26]
當參數取典型值α=10,γ=28,b=8/3時系統(tǒng)處于混沌狀態(tài).采用四階龍格-庫塔方法對方程進行求解,初值取[0.5,0.5,0.5],步長取0.01,取系統(tǒng)狀態(tài)變量x產生Lorenz混沌信號,并對信號加上零均值高斯白噪聲,產生信噪比SNR為15,20,25,30,35 dB的含噪Lorenz混沌信號.純凈信號和不同噪聲水平信號的PE如表1所列,隨著噪聲水平的增加,含噪混沌信號的復雜度也相應地增加.
表1 不同噪聲水平Lorenz混沌信號的排列熵值Table 1.PE of the Lorenz’s chaotic signal at different noise levels.
Chen系統(tǒng)方程為[27]
當參數取典型值a=35,b=3,c=28時系統(tǒng)處于混沌狀態(tài).采用四階龍格-庫塔方法對方程進行求解,初值取[0.5,0.5,0.5],步長取0.01,取系統(tǒng)狀態(tài)變量x產生Chen混沌信號,并對信號加上零均值高斯白噪聲,產生SNR為15,20,25,30,35 dB的含噪Chen混沌信號.純凈信號和不同噪聲水平信號的PE如表2所列,隨著噪聲水平的增加,含噪混沌信號的復雜度也相應地增加.
表2 不同噪聲水平Chen混沌信號的排列熵值Table 2.PE of the Chen’s chaotic signal at different noise levels.
文獻[17]給出的濾波參數選取范圍為:搜索窗長l取4000左右,搜索窗移動步長δ取塊寬w的1/2到1/4左右,分組塊數量m取30,塊寬w的取值不應小于100.然而,經實驗分析發(fā)現,在不同的信號特征、采樣頻率和噪聲水平情況下,對濾波參數搜索窗長、搜索窗移動步長和分組塊數量取文獻[17]給出的范圍能取得比較好的濾波效果,但關鍵濾波參數塊寬的選取卻不在所給的范圍之內,重構混沌信號中所含噪聲越少則其復雜度越低,本文基于排列熵算法的混沌信號協同濾波參數自動優(yōu)化準則為:取搜索窗長l=4000,搜索窗移動步長δ=w/4,分組塊數量m=30,采用不同塊寬w對含噪混沌信號進行去噪,計算重構混沌信號的PE,最小PE對應的塊寬即為最優(yōu)濾波塊寬.設重構混沌信號為(t),則最佳濾波塊寬wopt為
其中,符號HP[·]表示求排列熵.
信號特征、采樣頻率和噪聲水平影響濾波參數的選取,為了驗證本文所提方法的有效性,本節(jié)分別將該方法應用于不同信號特征、采樣頻率和噪聲水平情況下濾波參數的選取,為了方便分析,記濾波前后信號的信噪比分別為SNRin,SNRout.
經實驗分析發(fā)現,混沌信號協同濾波去噪算法具有低通濾波性質,濾波參數塊寬和低通截止頻率成反比關系.考慮到不同混沌信號的頻譜寬度不同,因此針對不同的混沌信號要選取不同的塊寬,以Lorenz混沌信號和Chen混沌信號為例,其頻譜如圖2所示.
由圖2可知,相對于Lorenz信號,Chen信號的頻譜更寬.由上述分析可知,將混沌信號協同濾波去噪算法應用于這兩類典型混沌信號時,對Chen信號應選取比Lorenz信號更窄的塊寬.分別給Lorenz信號和Chen信號加上零均值高斯白噪聲,信噪比取10 dB,采樣時間取0.01,取不同寬度的塊寬對這兩類含噪混沌信號進行去噪,并計算不同塊寬對應重構信號的排列熵.如圖3和圖4所示,重構信號的排列熵體現了對應濾波塊寬的濾波效果,在當前采樣頻率和噪聲水平情況下,Lorenz信號和Chen信號的最優(yōu)濾波塊寬分別為100和60,該結果與上述分析相符.
圖2 Lorenz混沌信號和Chen混沌信號頻譜Fig.2.Frequency spectra of the Lorenz’s chaotic signal and Chen’s chaotic signal.
圖3 含噪Lorenz信號去噪 (a)不同塊寬對應的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信號去噪 (a)不同塊寬對應的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.
混沌信號協同濾波去噪算法在高采樣頻率下能取得更好的去噪效果,對于不同的采樣頻率最優(yōu)濾波塊寬的選取也不同.采樣時間ts分別取0.02,0.01,0.0075,0.005產生Lorenz信號,給其加上零均值高斯白噪聲,信噪比取10 dB,取不同寬度的塊寬分別對不同采樣頻率的含噪信號進行去噪.如圖5和表3所示,采樣頻率越高去噪效果越好,最佳濾波塊寬也越寬.
表3 不同采樣頻率下的最優(yōu)濾波塊寬Table 3.Optimal block width at different sampling frequency.
圖5 采樣頻率對塊寬的影響 (a)不同采樣頻率對應的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.
混沌信號在頻域通常主要集中在低頻段,而白噪聲為均勻分布,由于混沌信號協同濾波去噪算法具有低通濾波性質,且濾波參數塊寬和低通截止頻率成反比關系,因此高信噪比情況下的最優(yōu)濾波塊寬比低信噪比情況下的最優(yōu)濾波塊寬稍窄采樣時間取0.01產生Lorenz信號,給其加上零均值高斯白噪聲,信噪比分別取5,10,15,20 dB,取不同寬度的塊寬分別對不同噪聲水平的含噪信號進行去噪.如圖6和表4所示,噪聲水平越高對應的最優(yōu)濾波塊寬越寬.
圖6 噪聲水平對塊寬的影響 (a)不同噪聲水平對應的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.
本文基于排列熵提出了一種混沌信號協同濾波參數自動優(yōu)化準則,深入分析了在不同信號特征、采樣頻率和噪聲水平情況下混沌信號協同濾波中各參數的選取問題.仿真結果表明,雖然濾波參數的選取受到不同信號特征、采樣頻率和噪聲水平的影響,新的參數優(yōu)化準則能有效實現不同情況下濾波參數的自動最優(yōu)化.該準則提高了混沌信號協同濾波去噪算法的自適應性,使其更符合實際應用需求,在各類混沌信號去噪中具有更廣泛的應用價值.
[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混沌時間序列分析及其應用(武漢:武漢大學出版社)第1—8頁]
[2]Han M,Xu M L 2013Acta Phys.Sin.62 120510(in Chinese)[韓敏,許美玲 2013物理學報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混沌信號與信息處理(北京:清華大學出版社)第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物理學報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物理學報 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物理學報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物理學報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