基于EMD-R/S分析的太赫茲光譜降噪

魏博熠 羅鑒鵬 張立臣

(廣東工業(yè)大學(xué)計(jì)算機(jī)學(xué)院 廣東 廣州 510006)

0 引 言

太赫茲時(shí)域光譜(Terahertz Time Domain Spectrum)可以展示物質(zhì)分子級別的化學(xué)鍵相互作用而引起的能級躍遷變化，反映了物質(zhì)豐富的物理化學(xué)等信息，通過太赫茲光導(dǎo)天線、太赫茲探測器等器件獲取太赫茲時(shí)域光譜信號，從而獲得被測物質(zhì)的信息。太赫茲時(shí)域光譜系統(tǒng)發(fā)射一個(gè)太赫茲脈沖，這個(gè)太赫茲脈沖一般包含半個(gè)到幾個(gè)振蕩周期。太赫茲探測器直接測量太赫茲波的電場隨時(shí)間的變化的信息。這個(gè)時(shí)域信息不但包含了光強(qiáng)，還測量了相位信息。此外，在太赫茲頻段范圍內(nèi)，對應(yīng)著有機(jī)物分子的偶極子轉(zhuǎn)動(dòng)和振動(dòng)躍遷以及大分子間的骨架振動(dòng)等[1-2]。根據(jù)這個(gè)特性，太赫茲光譜技術(shù)在毒品檢測[3-5]、食品質(zhì)量檢測[6-9]等領(lǐng)域有廣泛應(yīng)用。近年來，在太赫茲時(shí)域光譜系統(tǒng)的基礎(chǔ)上，進(jìn)一步研究太赫茲光譜的數(shù)據(jù)處理方法，對將太赫茲技術(shù)推向更廣闊的應(yīng)用領(lǐng)域有重大意義。

使用太赫茲時(shí)域光譜系統(tǒng)測量物質(zhì)，不可避免會(huì)混雜一些噪聲，對含噪信號進(jìn)行降噪處理很有必要。太赫茲光譜的噪聲可以分為背景噪聲和太赫茲波所攜帶的噪聲。這兩類噪聲表現(xiàn)為白噪聲和模糊噪聲的混合形式。目前對于太赫茲時(shí)域光譜信號的降噪處理主要是采用小波變換的方法[10-12]。小波變換的方法對于周期性的平穩(wěn)信號比較有效，而由于非平穩(wěn)信號包含尖峰或者突變信息，小波變換處理非平穩(wěn)信息有時(shí)效果欠佳。而且，小波變換處理信號需要選擇合適的小波基，對分解層數(shù)的選擇也需要依靠先驗(yàn)經(jīng)驗(yàn)，因此這種數(shù)據(jù)處理方法存在一定的局限性。經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition,EMD)是Huang等[13]在1998年提出的處理非線性和非平穩(wěn)時(shí)間序列信號的一種方法。與已有的方法相比，EMD直觀、直接并且能自適應(yīng)不同的信號，因?yàn)镋MD分解算法基于信號數(shù)據(jù)產(chǎn)生基函數(shù)。EMD基于一個(gè)簡單的假設(shè)：任何數(shù)據(jù)都包含了震蕩的不同簡單背部模態(tài)，所有不同的模態(tài)數(shù)據(jù)疊加最終得到復(fù)雜的數(shù)據(jù)。對時(shí)序信號采用EMD分解,結(jié)合信號物理背景，可以得到信號清晰的物理意義。

對太赫茲光譜信號進(jìn)行經(jīng)驗(yàn)?zāi)B(tài)分解后，可以得到若干個(gè)固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)分量，每個(gè)分量含有不同的頻率成分。本文使用Hurst指數(shù)識別出IMF分量中的白噪聲分量，將其從原始信號中剔除，從而獲得較為準(zhǔn)確的信號，方便后續(xù)處理?；赗/S分析法的Hurst指數(shù)是英國水文學(xué)家H.E Hurst建立的用于判斷時(shí)間序列數(shù)據(jù)遵從隨機(jī)游走還是有偏的隨機(jī)游走過程的指標(biāo)。統(tǒng)計(jì)學(xué)上一般認(rèn)為當(dāng)Hurst指數(shù)等于0.5時(shí)，該信號接近于高斯白噪聲信號。根據(jù)太赫茲噪聲的特點(diǎn)，采用Hurst指數(shù)識別噪聲會(huì)是一個(gè)合適的選擇。

本研究在使用EMD降噪的研究基礎(chǔ)上，提出了一種基于EMD和DFA相結(jié)合的太赫茲光譜降噪方法。該方法對含噪的太赫茲光譜信號進(jìn)行EMD分解，得到多個(gè)頻率信號分量。通過引入Hurst指數(shù)判斷各個(gè)IMF分量的回復(fù)性，將判斷為噪聲的信號從原始信號中剔除，獲得降噪后的太赫茲光譜信號。

1 EMD的理論基礎(chǔ)

EMD算法的目的是將信號分解成一系列的IMF,再通過希爾伯特變換獲得信號的時(shí)頻關(guān)系。經(jīng)驗(yàn)?zāi)B(tài)分解作為一種新的數(shù)據(jù)分析方法，在處理非平穩(wěn)、非線性數(shù)據(jù)上具有明顯的優(yōu)勢，特別適合處理太赫茲時(shí)域信號。在EMD算法中，IMF具有以下的定義：

(1) 整個(gè)數(shù)據(jù)中，極值點(diǎn)和過零點(diǎn)數(shù)目要么相等，要么最多相差1。

(2) 任意一點(diǎn)，由局部極大值定義的包絡(luò)和由局部極小值定義的包絡(luò)的均值等于0。

與簡單的諧波函數(shù)相比，IMF沒有簡單諧波函數(shù)中的恒定振幅和頻率，具有以時(shí)間為函數(shù)變化的振幅和頻率。

以太赫茲時(shí)域光譜信號為例，EMD算法的基本計(jì)算過程如下：

Step1計(jì)算太赫茲時(shí)域光譜信號x(t)的所有極大值點(diǎn)和極小值點(diǎn)，用三次樣條函數(shù)連接所有的極大值點(diǎn)，所得曲線構(gòu)成上包絡(luò)線u(t)；用三次樣條函數(shù)連接所有的極小值點(diǎn)，所得曲線構(gòu)成下包絡(luò)線v(t)。

Step2計(jì)算上下包絡(luò)線的均值，計(jì)算公式為：

(1)

EMD算法的成分分量定義為：

h(t)=x(t)-m(t)

(2)

Step3判斷成分分量h(t)是否滿足IMF的定義。如果h(t)滿足IMF的定義，則h(t)是第一個(gè)篩選出來的IMF分量，記為c1(t)。如果h(t)不滿足IMF的定義，則將h(t)作為原始數(shù)據(jù)，重復(fù)上述兩步，直至計(jì)算得到第一個(gè)IMF分量，該分量同樣記為c1(t)。

Step4將第一個(gè)IMF分量從原始信號中剔除，得到剩余信號，計(jì)算公式如下：

r1(t)=h(t)-c1(t)

(3)

將r1(t)作為新的待分析信號，重復(fù)上述3個(gè)步驟，這樣就得到第二個(gè)IMF分量c2(t)。將第二個(gè)IMF分量從r1(t)中剔除，得到新的余項(xiàng)r2(t)=r1(t)-c2(t)。不斷重復(fù)上述步驟，直至滿足EMD算法的停止準(zhǔn)則。EMD算法的停止準(zhǔn)則使用了Cauchy收斂性準(zhǔn)則，這種檢測要求相鄰兩個(gè)抽取操作的歸一化平方差足夠小，具體定義為：

(4)

式中：T表示信號長度。當(dāng)SDk小于設(shè)定的閾值時(shí)，分解結(jié)束。

2 R/S分析法

R/S分析法是一種用于處理時(shí)間序列的統(tǒng)計(jì)量的方法，其對于時(shí)間序列具備一定的預(yù)測功能，同時(shí)在分形理論中有著重要的作用。R/S分析法的基本思想是對一段時(shí)間序列數(shù)據(jù)進(jìn)行不同尺度的分解，研究其在不同的時(shí)間范圍內(nèi)的統(tǒng)計(jì)性變化規(guī)律，然后將小范圍內(nèi)的統(tǒng)計(jì)變化規(guī)律進(jìn)一步延伸至大范圍時(shí)間，得到整段時(shí)間序列數(shù)據(jù)的統(tǒng)計(jì)學(xué)規(guī)律。當(dāng)然，也有學(xué)者將大范圍時(shí)間變化規(guī)律用于小范圍時(shí)間預(yù)測方面。

太赫茲時(shí)域光譜信號是一段輻射強(qiáng)度隨時(shí)間變化的數(shù)據(jù)，這個(gè)信號數(shù)據(jù)可以表示成x(t1),x(t2),…,x(tn)的形式，也即為“時(shí)間序列”。臧元章等[14]指出，太赫茲時(shí)域光譜儀的幅值噪聲來源為延時(shí)線的重合抖動(dòng)、采樣抖動(dòng)以及飛秒激光抖動(dòng)，是白噪聲與1/f噪聲的混合。這種混合分形噪聲是非平穩(wěn)過程，理論上并不適合采用傅里葉變換、小波變換等譜分析方法，而EMD算法針對非平穩(wěn)信號，可以將特征信號與噪聲分離。此時(shí)，采用R/S分析法即可判斷哪些是信號，哪些是噪聲，然后重構(gòu)成高信噪比的信號。R/S分析法的實(shí)現(xiàn)步驟如下：

Step1對使用EMD算法分解后長度為N的IMF分量x(t)，以長度n均分這個(gè)IMF分量，得到m個(gè)相鄰的子區(qū)間。對這m個(gè)相鄰的子區(qū)間，可以表示為Ik，k=1,2,…,m。其中N=m×n。任一個(gè)子區(qū)間中的元素可以表示為N(k,n)，即表示任一個(gè)元素在第k個(gè)子區(qū)間，在該區(qū)間的第n位?？梢缘玫饺我蛔訁^(qū)間的均值為：

(5)

Step2計(jì)算任一子區(qū)間均值的累積離差：

(6)

其中，u=1,2,…,n。

Step3計(jì)算極差：

RIk=max(Xk)-min(Xk)

(7)

Step4計(jì)算子區(qū)間的標(biāo)準(zhǔn)差：

(8)

Step5對每一個(gè)RIk由對應(yīng)的SIk進(jìn)行標(biāo)準(zhǔn)化，得到RIk/SIk。

Step6n從4開始計(jì)算，重復(fù)上述5個(gè)步驟，直至n=N，即可得到序列[R/S]n,n=4,5,…,N。Hurst指數(shù)可以表示為：

(R/S)n=(c×n)H

(9)

式中：H表示Hurst指數(shù)；c為常數(shù)。

對式(11)兩邊取自然對數(shù)，可得：

(10)

Hurst指數(shù)是R/S分析的有效統(tǒng)計(jì)量。在分形理論中，H=0.5表明時(shí)間序列可以用隨機(jī)游走來描述；0≤H<0.5表明該時(shí)間序列是1/f噪聲，存在均值回復(fù)的特性；H>0.5表示時(shí)間序列存在長期記憶性。

3 基于EMD-R/S分析的光譜降噪實(shí)驗(yàn)

黃瑞瑞等[15]進(jìn)行了乳糖的太赫茲光譜和紅外光譜的特性研究。本文采用乳糖作為實(shí)驗(yàn)材料進(jìn)行研究，將本文的實(shí)驗(yàn)數(shù)據(jù)與文獻(xiàn)中的實(shí)驗(yàn)結(jié)果進(jìn)行對比，本文測量的乳糖在頻域相似位置出現(xiàn)了吸收峰，證明測量過程是準(zhǔn)確的。以厚度為1.3 mm的乳糖壓片的太赫茲透射譜為例，測量時(shí)溫度為23 ℃，濕度為3.4%。測量得到的太赫茲時(shí)域信號如圖1所示。

圖1 乳糖壓片太赫茲時(shí)域光譜

對乳糖壓片的太赫茲時(shí)域光譜信號采用EMD算法分解得到如圖2所示的結(jié)果。從分解結(jié)果中看到，IMF2分量與原始信號形狀相似，保留了原始信號的大部分信息。其余分量代表了噪聲等不同的物理信息。

圖2 EMD分解結(jié)果

對EMD算法分解得到的IMF分量逐個(gè)使用R/S分析法計(jì)算Hurst指數(shù)。各個(gè)IMF分量的Hurst指數(shù)分別為0.42、0.69、0.76、0.66、0.84、0.79、0.91、0.96、0.99、0.99。從Hurst指數(shù)計(jì)算結(jié)果，可以得到IMF1分量是信號中的噪聲。將噪聲信號從原始信號中剔除，可以得到重建信號。為了方便對比，本文將重建信號向下平移了0.1個(gè)單位，效果如圖3所示。重建后的信號跟原始信號相比更加平滑，毛刺更少。

圖3 原始時(shí)域信號與重建時(shí)域信號對比

對物質(zhì)的太赫茲信號分析常常不局限于時(shí)域信號，不少學(xué)者會(huì)將時(shí)域信號轉(zhuǎn)換到頻域進(jìn)行分析。實(shí)驗(yàn)分別對原始信號與重建信號做傅里葉變換。為了方便觀察，重建頻域信號進(jìn)行了向下平移，所得結(jié)果如圖4所示。與原始頻域信號相比，重建頻域信號更平滑。

圖4 原始頻域信號與重建頻域信號對比

4 算法性能評價(jià)

實(shí)驗(yàn)采用小波軟閾值降噪、小波硬閾值降噪和本文提出的方法分別對20個(gè)不同厚度的乳糖壓片進(jìn)行了太赫茲光譜測量，分別計(jì)算原始時(shí)域信號和重建時(shí)域信號前后的相關(guān)系數(shù)、信噪比和峰值信噪比三個(gè)指標(biāo)。本文測量20個(gè)實(shí)驗(yàn)室制造的乳糖壓片，厚度范圍為0.8～1.6 mm，每個(gè)乳糖壓片表面光滑。 乳糖壓片的厚度使用電子游標(biāo)卡尺測量，對同一個(gè)乳糖壓片測量5次取其平均值作為最終測量結(jié)果，測量誤差為0.02 mm。小波去噪常用的方法是小波閾值去噪法，這是本文實(shí)驗(yàn)采用的對比方法。小波閾值去噪法的原理是將信號通過小波分解后，有用的信號小波系數(shù)幅值大，能量大，小波系數(shù)數(shù)量偏少；噪聲信號小波系數(shù)幅值小，能量小，小波系數(shù)多且分布均勻。小波閾值去噪方法可細(xì)分為軟閾值和硬閾值兩種。李鵬鵬[16]使用小波分析對太赫茲信號分析得出小波分解的小波基函數(shù)、分解層數(shù)和閾值選擇對去噪效果有影響，其中的最優(yōu)選擇需要依靠經(jīng)驗(yàn)。本文對乳糖的太赫茲光譜小波分解后，發(fā)現(xiàn)使用db2小波基、分解層數(shù)為兩層時(shí)降噪效果比較好。因此，小波降噪法采用的小波基是db2函數(shù)，分解層數(shù)為兩層。小波降噪法采用的軟閾值通過無偏風(fēng)險(xiǎn)估計(jì)閾值計(jì)算得到，硬閾值采用固定閾值公式結(jié)合經(jīng)驗(yàn)分析取得。固定閾值公式為：

(11)

式中：N表示輸入信號的長度。實(shí)驗(yàn)過程中，本文使用單一控制法控制各種變量，即盡量保證實(shí)驗(yàn)環(huán)境不變，測量不同厚度的乳糖壓片。不同降噪算法的平均效果對比如表1所示。

表1 不同降噪算法的平均效果

由表1可以看出，本文提出的方法在信噪比和峰值信噪比的指標(biāo)上優(yōu)于小波降噪方法。主要的原因是小波降噪依靠人工經(jīng)驗(yàn)選取閾值，閾值選取的效果直接影響降噪算法的效果。而本文提出的方法物理意義清晰，可操作性強(qiáng)，能夠很好地保留信號原有的特征信息。

5 結(jié) 語

本文提出了采用EMD-R/S分析的方法去除太赫茲時(shí)域光譜中的噪聲。該方法利用EMD算法自適應(yīng)分解太赫茲時(shí)域光譜信號，然后利用Hurst指數(shù)對分解后的IMF分量進(jìn)行噪聲判斷，最后將噪聲從原始信號中剔除。該方法屬于數(shù)據(jù)驅(qū)動(dòng)的方法，克服了小波降噪方法依靠人工經(jīng)驗(yàn)決定降噪效果的弊端，物理意義清晰，可操作性強(qiáng)。這對后期的太赫茲時(shí)域光譜應(yīng)用推廣有一定的參考價(jià)值。