脈沖星信號的經驗模態分解模態單元比例萎縮消噪算法*

王文波 張曉東 汪祥莉

1 引言

脈沖星是一種高密度、高速自轉且具有強磁場的中子星,其最重要的特征是自轉具有超高的穩定性和均勻性,其毫秒級的脈沖輻射周期被稱為是自然界中最穩定的時鐘[1].近年來,脈沖星在計時、定位、導航、電波傳播學和天體物理學等多個領域日益展現出巨大的應用潛力,受到越來越廣泛的關注和研究[2-4].然而,由于脈沖星在距離地球上千甚至幾十萬光年的遙遠深空中,其輻射脈沖信號在傳播和接收過程中不可避免地受到空間介質、宇宙輻射、大氣層色散延遲以及地球運動等眾多因素的影響,導致所接收到的脈沖輻射信號極弱且常常淹沒在噪聲之中[5,6].因此,對脈沖輻射信號中的噪聲進行有效濾除將直接影響到脈沖累積輪廓和脈沖到達時間的計算精度,是脈沖星信號后繼研究、應用最根本的前提和基礎.

小波分析由于具有良好的時頻局域化特性、多分辨率、去相關性等特點,被廣泛應用于脈沖星信號的消噪中,取得了較好的消噪效果[7,8].但在應用小波方法消噪時,需要預先選定小波基和分解層數,相同條件下選用不同的小波基和分解層數,對去噪結果影響很大,特別是小波基函數的選擇,對去噪結果有決定性的影響.而如何選擇小波基和最優分解層數是小波分解中一個很難解決的問題,這給利用小波進行信號去噪帶來了很大的不便[9,10].經驗模態分解(empiricalmodedecomposition,EMD)是一種新的完全數據驅動的自適應信號分解算法[11],可以把數據分解成具有物理意義的一組內蘊模態函數(intrinsic mode function,IMF)分量.EMD在分解信號時,不需要預先確定基底和分解層數,而是根據信號自身特性,自適應地確定基底和最優分解層數.對于非線性和非平穩信號,EMD的分解結果比小波分解結果更清晰、準確[12,13],分解出的IMFs能夠充分保留信號本身所固有的非線性和非平穩特征.

對于非線性和非平穩信號,EMD消噪可以獲得比小波消噪更好的效果[14-17],脈沖星信號是典型的非平穩信號[18],本文將EMD應用于脈沖星信號的消噪中.利用EMD對非線性、非平穩信號消噪已得到了較多的研究.Boudraa和Cexus[15]提出了部分重構的EMD消噪算法,但部分重構消噪算法中將前幾項IMF當作噪聲項直接刪除,丟失了較多的細節信息而且噪聲不能被有效去除.Olufemi等[16]提出了基于系數閾值的EMD消噪算法,但基于系數閾值的消噪算法破壞了IMF中模態單元的完整性,影響了去噪的效果.文獻[17,19]提出了基于模態單元閾值的EMD消噪算法,將IMF中兩個過零點間的波動單元-模態單元作為一個整體,以模態單元作為基本單位,利用閾值方法進行去噪,保持了IMF中模態單元的完整性,有效提高了去噪效果.但該算法中對極值小于閾值的模態單元直接刪除,不考慮其中所包含的細節信息;對于極值大于閾值的模態單元直接保留,不考慮其中所包含的噪聲,勢必導致部分信號細節信息被誤刪且噪聲不能被完整去除,影響了去噪效果的進一步提高.本文在EMD模態單元閾值去噪算法的基礎上,提出了一種基于模態單元比例萎縮的EMD去噪方法,并應用于脈沖星信號的消噪處理.所提出的方法中以模態單元為基本單位構造最優比例萎縮因子,對IMF中的每個模態單元都進行比例萎縮去噪,在保持模態單元完整性的基礎上,有效去除每個模態單元中所含的噪聲.分別采用小波硬閾值去噪、比例萎縮小波去噪、EMD模態單元閾值去噪和本文方法對脈沖星信號進行消噪對比實驗,結果表明本文方法在噪聲去除和信號細節特征保持兩方面都有較好的提高.

2 基于模態單元比例萎縮的EMD去噪模型

EMD的主要目的是將待分析信號分解為一系列表征時間尺度的IMF分量,要求IMF分量必須滿足兩個條件：IMF的極值點個數與過零點個數不超過1;由極大值點和極小值點確定的包絡線均值為零.EMD分解結束后,原始信號x(t)可被表示為一組IMF和一個余項之和：

其中imfk表示第k個IMF分量,rK(t)表示余項,K表示分解出的IMF的總層數.每個IMF的極大、極小值點形成的上、下包絡線關于時間軸局部對稱,反映了信號內部的固有波動特性.

在IMF中,一次完整的波動由三個過零點和一個波峰、一個波谷構成,相鄰的兩個過零點之間只有一個極值點.由兩個過零點和一個極值點所構成的信號部分是IMF最基本的組成單元,它是構成信號固有波動的最小完整單元,一般將其稱為IMF的模態單元[16],模態單元中極值點的系數絕對值稱為模態振幅.對于第k層內蘊模態函數imfk,如果把imfk中的第i個模態單元記為z(k)i,則imfk可被表示為其中 k=1,2,···,K,本文中模態單元的模態振幅記為,的極值點記為.

在基于系數閾值的EMD去噪算法中[19],借鑒小波閾值去噪的思想構造imfk的去噪閾值Tk,然后利用該閾值對imfk系數進行處理.若令hk(t)=imfk(t),基于IMF系數的硬閾值去噪算法可表示為

在基于IMF系數的閾值去噪算法中,IMF中所有小于閾值的系數都被置0,導致模態振幅較大的模態單元中的部分點也被置零(如圖1(c)所示),破壞了模態單元局部固有波動的完整性,影響了去噪的效果.

基于模態單元閾值的EMD去噪算法中[16,20],將模態單元作為基本單位,針對模態單元構造閾值Tk,利用Tk對imfk中的模態單元進行閾值處理.基于模態單元的閾值去噪法可表示為

可以看出在該算法中,模態振幅小于閾值的模態單元被置零,而模態振幅大于閾值的模態單元被完整保留(如圖1(d)所示),沒有破壞IMF中模態單元局部固有波動的完整性,較好地提高了EMD的去噪效果.

但文獻[10,17]的研究表明,imfk的每個部分都分布有不同程度的噪聲：模態振幅較大的模態單元主要由信號的波動構成,但仍含有少量噪聲,直接保留會導致噪聲不能完整去除;而模態振幅較小的模態單元雖然主要由噪聲的波動構成,但仍含有部分信號細節信息,直接置零將不可避免地會損失部分信號細節.為了進一步提高EMD的去噪效果,本文借鑒小波比例萎縮去噪的思想,以模態單元為基本單位,構造合適的比例萎縮因子,對IMF中的每個模態單元進行比例萎縮去噪.

假設hk=imfk=中模態單元zik的比例萎縮因子為θik,則基于模態單元的比例萎縮去噪方法可表示為

基于模態單元比例萎縮的去噪算法沒有破壞模態單元局部固有波動的完整性,而且對振幅較大的模態單元和振幅較小的模態單元都進行了相應的萎縮去噪處理(如圖1(e)所示),在一定程度上改善了模態單元閾值去噪算法的不足,可以更好地去除噪聲并保留信號細節信息.信號經EMD模態單元比例萎縮去噪后可表示為

式中,Kd表示需要濾波的IMF數目,對于高斯白噪聲,通常取Kd=min(8,K),即去噪時僅對前8個IMF進行濾波,而其余的IMF直接保留.在基于模態單元的比例萎縮去噪中,構造合適的模態單元比例萎縮因子θki非常關鍵.若想構造合適的比例萎縮因子,最大程度的恢復原始信號,必須了解IMF的統計模型.本文在下一節中對含噪脈沖星信號IMF的統計特性進行討論.

圖1 EMD去噪算法模型 (a)測試信號Doppler第4層IMF(imf4);(b)imf4中所選部分放大圖;(c)系數閾值去噪;(d)模態單元閾值去噪;(e)模態單元比例萎縮去噪

3 含噪脈沖星信號IMF的統計特性分析

3.1 含噪脈沖星信號IMF中噪聲分布形式

含噪脈沖星輻射信號y經消色散處理后,可以近似看作是有用信號與高斯白噪聲混合組成,即

式中,x為有用脈沖星信號,d為服從正態分布的高斯白噪聲.含加性高斯白噪聲的信號經EMD分解后,IMF中噪聲的分布形態尚沒有完善的理論體系,但通過大量的實驗分析,一般認為IMF中的噪聲仍服從加性分布[13,21].因此,脈沖星輻射信號經EMD分解后,第k層內蘊模態函數hk=imfk可設為

其中sk表示imfk中所含的信號信息,nk表示所含噪聲,且sk和nk相互獨立.

3.2 含噪脈沖星信號IMF中信號項和噪聲項的統計模型

為了對imfk中的模態單元構造最優比例萎縮因子,需要知道信號項sk和噪聲項nk的分布形式.純零均值高斯白噪聲經EMD分解后,其IMF系數仍服從零均值正態分布[17],即

式中σ2k,n為nk的方差.但imfk中信號項部分sk的分布形式尚沒有具體的研究結果,本文通過實驗的方式分析sk的統計特性.由于EMD分解后IMF中信號和噪聲混雜在一起,難以直接對信號項的統計特性進行分析,本文從含噪脈沖星信號的IMF系數hk出發,對信號項sk的統計特性進行間接分析.

在脈沖星信號B1953+29的標準輪廓中添加白噪聲后進行50次EMD分解,取50次分解結果的平均值作為最終結果,統計各層IMF的系數分布,并與高斯分布進行比較.圖2給出了含噪脈沖星信號第2—5層IMF系數的直方圖和相應的高斯分布模擬逼近曲線,直觀上可以看出高斯分布可以較好地逼近IMF系數的分布.為了進一步驗證零均值高斯分布對含噪脈沖星信號IMF系數的逼近程度,對IMF的系數分布進行偏度、峰度檢驗[22].偏度、峰度檢驗法是正態性檢驗的經典方法,隨機變量S的偏度v1和峰度v2是指S的標準變化量[S-E(S)]/的三階和四階矩,即

其中,E(S)和D(S)分別表示S的期望和方差.當隨機變量S服從正態分布時,v1=0且v2=3.設S1,S2,···,Sn是來自總體S的樣本,則樣本偏度v1和樣本峰度v2可由下式進行簡單計算：v1=B3/B,v2=B4/B22,式中Bk(k=2,3,4)表示樣本的k階中心矩.

圖2 含噪脈沖星信號IMF系數的高斯逼近 (a)imf2系數直方圖及高斯逼近;(b)imf3系數直方圖及高斯逼近;(c)imf4系數直方圖及高斯逼近;(d)imf5系數直方圖及高斯逼近

如果含噪脈沖星信號的IMF系數符合零均值高斯分布,則其樣本偏度接近于0,而樣本峰度接近于3,計算其第2—5層IMF的樣本偏度和樣本峰度,結果如表1所示.

表1 含噪脈沖星信號IMF的樣本偏度和峰度

從表1可以看出,含噪信號各層IMF的樣本偏度和樣本峰度都比較接近于0和3,樣本偏度絕對值的均值約為0.0167,樣本峰度的平均值約為3.0013.因此,高斯模型可以較好地模擬含噪脈沖星信號IMF系數的分布特性.本文對B2351+61和B1953+29等脈沖星信號分別添加不同強度的白噪聲進行多次實驗,得到了類似的實驗結果,所以可假設含噪脈沖星信號的IMF系數近似服從零均值高斯分布,即hk～N(0,σ2h,k).由于sk=hk-nk,且hk,nk都服從零均值高斯分布,由概率知識可知imfk中信號項部分sk也服從零均值高斯分布,即

4 模態單元比例萎縮因子的確定

本文采用線性最小均方誤差(linear minimum mean square error estimation,LMMSE)準則計算imfk去噪時的最優比例萎縮因子[22].令hk=imfk,由(2)式可知hk=sk+nk,設去噪時的比例萎縮因子為θk,則去噪后的信號

為了使?hk與原始真實信號sk間的均方誤差E[(sk-?hk)2]最小,根據LMMSE理論,應有E[(sk-?hk)hTk]=0,從而可求出

考慮到信號sk與噪聲nk都服從零均值高斯分布且相互獨立,所以

因此

但如果對imfk中的每個系數都按照自身的比例萎縮因子進行去噪,會造成一個模態單元內出現多個極值的情況,將破壞模態單元固有波動的完整性.在模態單元zik內,極值點e(ik)反映了模態單元最本質的波動特征和信息,因此可將極值點e(ik)的萎縮因子θ(e(ik))作為整個模態單元zik的萎縮因子,按此方法對imfk進行比例萎縮去噪,即

此時,由于整個模態單元zki內采用相同的萎縮因子θ(e(k)i),因此不會破壞模態單元局部固有波動的完整性.本文中噪聲標準差σn,k采用(5)式進行估計,

式中HH表示hk的高頻子帶小波系數,而信號方差σ2s,k采用文獻[24]中提出的方法進行估計,即

這里M表示所選窗口中信號的長度,本文中取M=5.

5 實驗分析

本文中選用B1953+29和B2351+61脈沖星輻射信號進行實驗分析,數據來源于歐洲脈沖星網絡數據庫EPN(the European Pulsar Network Data Archive).脈沖星信號B1953+29和B2351+61的標準輪廓如圖3(a)和圖4(a)所示,在標準輪廓中分別添加一定量的高斯白噪聲模擬含噪脈沖星信號,B1953+29加噪后的信號如圖3(b)所示,其信噪比為10;B2351+61加噪后的信號如圖4(b)所示,其信噪比為20.分別采用硬閾值小波去噪法[24](wavelet hard threshold denoising,WHD)、比例萎縮小波去噪法[23](wavelet proportion shrinking denoising,WPSD)、EMD模態單元閾值去噪法(mode-cell threshold denoising,MTD)和本文提出的EMD模態單元比例萎縮去噪法(mode-cell proportion shrinking denoising,MSD)分別對含噪脈沖星信號進行消噪處理,消噪后的信號如圖3(c)—(f)和圖4(c)—(f)所示.對消噪結果采用以下四個參數進行聯合評價[7]：信噪比、均方根誤差、峰值相對誤差、峰位誤差.其中峰值相對誤差和峰位誤差的定義公式如下.

1)峰值相對誤差：REPV=(Vo-Vd)/Vo·100%,式中Vo,Vd分別表示脈沖星輻射信號標準輪廓的脈沖峰值和消噪后信號的脈沖峰值.

2)峰位誤差：EPP=|Po-Pd|,式中Po,Pd分別表示脈沖星輻射信號標準輪廓波的脈沖峰位值和消噪后輻射信號的脈沖峰位值.

首先從視覺效果上對四種消噪方法進行比較.從圖3(c)和圖4(c)可以看出,WHD算法較好地去除了脈沖星信號中的噪聲,消噪后信號的平滑部分與標準輪廓基本符合,但由于對小波閾值的系數全部去除,使得在去除噪聲的同時也損失了較多的細節信息,導致信號的突變部分去噪后與標準輪廓之間的符合程度較差.WPSD算法的去噪效果要優于WHD(如圖3(d)和圖4(d)所示),消噪后信號的平滑部分與標準輪廓符合較好,突變部分與標準輪廓的符合度有一定程度的提高.MTD的消噪效果(如圖3(e)和圖4(e)所示)與WPSD比較接近,在平滑部分和幾個大的波峰處,兩種方法去噪都取得了較好的效果,但在信號劇烈振蕩的地方,MTD比WPSD具有更好的符合效果.本文方法消噪后的結果如圖3(f)和圖4(f)所示,通過比較可以看出,本文方法消噪后平滑部分與WPSD和MTD基本相同,但在脈沖峰值和劇烈振蕩的突變點處,與標準輪廓有更好的逼近、誤差更小,減小了脈沖峰高度值和位置值的偏差,突變點處的細節信息損失更小.因此從視覺分析來看,本文方法去噪效果優于WPSD和MTD,更好地保留了脈沖星信號的細節信息,在信號的突變點處與標準輪廓具有更好的吻合度.

圖3 含噪脈沖星信號B1593+29消噪結果比較 (a)脈沖星信號標準輪廓;(b)含噪脈沖星信號;(c)WHD消噪;(d)WPSD消噪;(e)MTD EMD消噪;(f)本文方法消噪

表2 不同方法消噪后的參數比較

圖4 含噪脈沖星信號B2351+61消噪結果比較 (a)脈沖星信號標準輪廓;(b)含噪脈沖星信號;(c)WHD消噪;(d)WPSD消噪;(e)MTD EMD消噪;(f)本文方法消噪

再對四種方法消噪后的參數進行比較,當消噪后信號的信噪比(SNR)越大,而均方根誤差(RMSE),REPV和EPP越小時,表明消噪效果越好.不同方法消噪后的SNR,RMSR,REPV和EPP如表2所示.從表2可以看出,WHD的效果稍差,除了峰位誤差其他三項指標都不如另外三種方法.WPSD的參數與EMD,MTD相差不大,MTD消噪后的均方根誤差和峰值相對誤差略優于WPSD.本文方法消噪后的脈沖星信號信噪比最大、均方根誤差最小,表明本文方法消噪后的信號能更好地逼近脈沖星標準輪廓;本文方法消噪后的峰位誤差與另外三種方法基本相等,但峰值相對誤差小于另三種方法,表明本文方法去除噪聲的同時較好地保留了信號中脈沖尖峰部分的有用細節信息,能更好地減小脈沖峰值處重要信息的損失.整體比較可知,本文方法取得了更好的消噪效果,在噪聲抑制和脈沖尖峰重要信息的保持等方面都有一定程度的提高,去噪后的參數指標優于WPSD和MTD.

6 結論

將一種基于模態單元比例萎縮的EMD濾波方法應用到脈沖星信號的消噪處理中.脈沖星信號經EMD分解后,通過分析其IMF中信號項和噪聲項的統計特性,以模態單元為基本單位在最小均方誤差準則下構造最優比例萎縮因子,對IMF進行比例萎縮降噪.實驗結果表明：與硬閾值小波消噪法、比例萎縮小波消噪法和基于模態閾值的EMD消噪法相比,本文方法能獲得更好的消噪效果,可以在有效濾除噪聲的同時更好地保留原信號中有用的細節信息,脈沖星信號經本文方法消噪后,具有更高的信噪比和更低的均方根誤差、峰值相對誤差.本文在分析脈沖星信號IMF中噪聲的分布時采用的是加性噪聲模型,尚缺少完善的理論證明,通過理論分析研究IMF中噪聲的分布形式是今后研究的重點;本文中對IMF中信號項部分采用高斯分布進行擬合,是否可以采用其他分布如廣義高斯分布、混合高斯分布等對其進行更精確的擬合逼近,也是今后將繼續探討的一個方向.

[1]Taylor J H 1991 Proc.IEEE 79 1054

[2]Xie Q,Xu L P,Zhang H,Luo N 2012 Acta Phys.Sin.61 119701(in Chinese)[謝強,許錄平,張華,羅楠2012物理學報61 119701]

[3]Sheikh S I,Pines D J,Ray P S,Wood K S,Michael N L,Wolff M T 2006 J.Guidance,Control and Dynamics 29 49

[4]Zhu J,Li P Y 2008 Chin Phys.B 17 356

[5]Zhao B S,Hu H J,Sheng L Z,Yan Q R 2011 Acta Phys.Sin.60 029701(in chinese)[趙寶升,胡慧君,盛立志,鄢秋榮2011物理學報60 029701]

[6]Hu H J,Zhao B S,Sheng L Z,Yan Q R,Yang H,Chen B M 2011 Sci.Sin.Phys.Mech.&Astron.41 1015(in chinese)[胡慧君,趙寶升,盛立志,鄢秋榮,楊顥,陳寶梅2011中國科學：物理學、力學、天文學41 1015]

[7]Yan D,Xu L P,Xie Z H 2007 J.Xi’an Jiaotong Univ.41 1193(in chinese)[閻迪,徐錄平,謝振華2007西安交通大學學報41 1193]

[8]Gao G R,Liu Y P,Pan Q 2012 Acta Phys.Sin.61 139701(in Chinese)[高國榮,劉艷萍,潘瓊2012物理學報61 139701]

[9]Zhang H,Xu L P,Xie Q,Luo N 2011 Acta Phys.Sin.60 049701(in Chinese)[張華,許錄平,謝強,羅楠2011物理學報60 049701]

[10]Li Y Q,Li P,Yan X P,Chen H M 2008 Transactions of Beijing Institute of Technology 28 723(in chinese)[李月琴,粟蘋,閆曉鵬,陳慧敏2008北京理工大學學報28 723]

[11]Zhang H,Chen X H,Yang H Y 2011 OGP 46 70(in Chinese)[張華,陳小宏,楊海燕2011石油地球物理勘探46 70]

[12]Huang N E,Shen Z,Long S R,Wu M C,Shih H H,Zheng Q,Yen N C,Tung C C,Liu H H 1998 Proc.of the Royal Society of London A 454 903

[13]Flandrin P,Paulo G 2004 Int.J.Wavel.Multiresol.Inform.Proc.2 1

[14]Gao X Q,Dong W J,Gong Z Q,Zou M W 2005 Acta Phys.Sin.54 3948(in Chinese)[高新全,董文杰,龔志強,鄒明瑋2005物理學報54 3948]

[15]Boudraa A,Cexus J 2007 IEEE Trans.Instrum.Measur.56 2196

[16]Olufemi A,Vladimir A,Auroop R 2011 IEEE Sens.J.11 2565

[17]Kopsinis K,Mclaughlin S 2009 IEEE Trans.Signal Proc.57 1351

[18]Wu Y F,Qiu Y,Yang Y F,Ren X M 2010 Acta Phys.Sin.59 3778(in Chinese)[吳亞鋒,裘炎,楊永鋒,任興民2010物理學報59 3778]

[19]Qu C S,Lu T Z,Tan Y 2010 Acta Autom.Sin.36 67(in Chinese)[曲從善,路廷鎮,譚營2010自動化學報36 67]

[20]Jenet F A,Anderson S B 1998 Publica.Astrono.Soc.Pacif i c 100 1467

[21]Wu Z H,Norden E H 2004 Proc.R.Soc.Lond.A 460 1597

[22]Sheng Z,Xie S Q,Pan C Y 2006 Probability Theory and Mathematical Statistics(Beijing：Higher Education Press)p350(in Chinese)[盛驟,謝式千,潘承毅2006概率論與數理統計(北京：高等教育出版社)第350頁]

[23]Chang S,Yu B,Vetterli M 2000 IEEE Ttrans.Image Proc.9 1532

[24]Marian K 2003 IEEE Signal Proc.Lett.10 324

[25]Donoho D L,Johnstone I M 1995 J.Am.Statist.Association 90 1200