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

魏博熠 羅鑒鵬 張立臣

(廣東工業大學計算機學院 廣東 廣州 510006)

0 引 言

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

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

對太赫茲光譜信號進行經驗模態分解后，可以得到若干個固有模態函數(Intrinsic Mode Function，IMF)分量，每個分量含有不同的頻率成分。本文使用Hurst指數識別出IMF分量中的白噪聲分量，將其從原始信號中剔除，從而獲得較為準確的信號，方便后續處理。基于R/S分析法的Hurst指數是英國水文學家H.E Hurst建立的用于判斷時間序列數據遵從隨機游走還是有偏的隨機游走過程的指標。統計學上一般認為當Hurst指數等于0.5時，該信號接近于高斯白噪聲信號。根據太赫茲噪聲的特點，采用Hurst指數識別噪聲會是一個合適的選擇。

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

1 EMD的理論基礎

EMD算法的目的是將信號分解成一系列的IMF,再通過希爾伯特變換獲得信號的時頻關系。經驗模態分解作為一種新的數據分析方法，在處理非平穩、非線性數據上具有明顯的優勢，特別適合處理太赫茲時域信號。在EMD算法中，IMF具有以下的定義：

(1) 整個數據中，極值點和過零點數目要么相等，要么最多相差1。

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

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

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

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

Step2計算上下包絡線的均值，計算公式為：

(1)

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

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

(2)

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

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

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

(3)

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

(4)

式中：T表示信號長度。當SDk小于設定的閾值時，分解結束。

2 R/S分析法

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

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

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

(5)

Step2計算任一子區間均值的累積離差：

(6)

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

Step3計算極差：

RIk=max(Xk)-min(Xk)

(7)

Step4計算子區間的標準差：

(8)

Step5對每一個RIk由對應的SIk進行標準化，得到RIk/SIk。

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

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

(9)

式中：H表示Hurst指數；c為常數。

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

(10)

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

3 基于EMD-R/S分析的光譜降噪實驗

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

圖1 乳糖壓片太赫茲時域光譜

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

圖2 EMD分解結果

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

圖3 原始時域信號與重建時域信號對比

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

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

4 算法性能評價

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

(11)

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

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

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

5 結 語

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