改進的自適應噪聲總體集合經驗模態分解在光譜信號去噪中的應用

李曉莉，李成偉

(哈爾濱工業大學 電氣工程及自動化學院，黑龍江 哈爾濱 150001)

?

改進的自適應噪聲總體集合經驗模態分解在光譜信號去噪中的應用

李曉莉，李成偉*

(哈爾濱工業大學 電氣工程及自動化學院，黑龍江 哈爾濱 150001)

針對近紅外無創血糖檢測過程中噪聲對血糖濃度模型精度和穩定性的影響，提出用自適應噪聲總體集合經驗模態分解方法實現近紅外光譜信號的去噪；同時，根據原始信號曲率和分解后本征模態函數(IMFs)曲率間的離散弗雷歇距離選擇相關模態。首先，將自適應噪聲的總體集合經驗模態分解方法引入近紅外光譜去噪過程，介紹了經驗模態分解、集合經驗模態分解、互補集合經驗模態分解及自適應噪聲總體集合經驗模態分解的基本原理及具體實現過程。然后，應用基于曲率和離散弗雷歇距離的自適應噪聲總體集合經驗模態分解改進算法對仿真信號和光譜信號進行去噪，并將其標準差和信噪比作為評價指標。實驗結果表明：應用提出的方法得到的血糖濃度近紅外光譜數據其標準差為0.179 4，信噪比為19.117 5 dB，實現了信號與噪聲的分離，改善了重構信號質量，具有良好的自適應性，可以有效識別并提取有用信息。

無創血糖檢測；近紅外光譜；信號去噪；自適應噪聲總體集合經驗模態分解；曲率；離散弗雷歇距離

1　引　言

糖尿病是威脅人類健康的重大疾病之一，而且糖尿病患者的數量正以驚人的速度增長，目前，世界上已經有超過一億的人患有糖尿病，預計到2030年將增加到3億。糖尿病患者雖然可以通過調整飲食或是注射胰島素來調節血糖水平，但糖尿病后期會引起嚴重的并發癥，例如：心臟衰竭和失明。所以，糖尿病的預防非常重要，無創血糖檢測方法由于可以減少頻繁檢測所帶來的疼痛感并降低醫療成本，故得到廣泛應用。無創血糖檢測技術包括吸收光譜法、光聲光譜法、光偏振法、熒光法和介電光譜法[1-3]。Arlen Duncan等在1995年利用脈沖光聲光譜法檢測手指血液內的葡萄糖濃度[2]。J.R. Blanco等在2006年研制出一種便攜式葡萄糖傳感器，該設備通過測量特定物質和生物識別系統之間電子交換產生的電流檢測葡萄糖濃度[3]。2011年，Ashok等利用反式激光束測量糖尿病患者的血糖濃度[4]?；诮t外光譜進行無創血糖檢測已經成為國內外專家學者的研究熱點[1]。

在近紅外光譜測量中，噪聲是一個重大的挑戰。目前存在的去噪方法主要有基于模型方法、變換域方法和自適應濾波方法[5-7]。維納濾波器由于具有易于實現和設計的優點，已經被廣泛應用，但是這種線性方法只能用于穩定信號。為了克服這個限制，有學者提出了基于小波閾值的非線性方法。然而小波方法的基函數是固定的，不能對所有的真實信號都匹配。準確地說，這種局限性源于它的非自適應性，一旦確定好基小波，它將被用于分析所有數據，如果選擇的小波分解不合適將會限制其去噪性能。經驗模態分解(Empirical Mode Decomposition， EMD)常被用于分析非穩態和非線性數據[8-9]，它可以將任意信號分解成不同的振蕩成分，即本征模態函數(Intrinsic Mode Functions， IMFs)。這個強大的自適應工具非常適合解決測量領域中的噪聲和頻率估計問題，但EMD是基于相關模態的局部重建，用其進行信號濾波屬于自適應方法，故周期性信號引起模態混疊的問題。集合經驗模態分解(Ensemble Empirical Mode Decomposition， EEMD)可有效地解決這個問題，但同時也引入了新的問題，即在重構信號中含有殘留噪聲。目前，EMD和EEMD兩種方法已經廣泛應用于各個領域[10-14]?；パa集合經驗模態分解(Complementary Ensemble Empirical Mode Decomposition， CEEMD)可以抑制重構余項并消除IMF里的殘余噪聲，但是當參數選擇不當時，會產生錯誤成分導致最后獲得的IMF不能真正滿足IMF的定義。為了克服這個問題，本文提出了自適應噪聲的總體集合經驗模態分解(Complete Ensemble Empirical Mode Decomposition with Adaptive Noise， CEEMDAN)算法來獲得準確的原始信號重構和純凈的分解模態光譜。

2　EMD方法及其改進版本

2.1EMD方法

EMD通常將信號分解為小數量的IMFs，IMF應該滿足以下2個條件：(1)極值點和過零點的數量相等或最多相差一個點。(2)任一點的上下包絡線均值需為零[15]。EMD是一種通過連續減去包絡均值去除震蕩的自適應方法。對于信號x(t) 而言，EMD算法包括以下步驟[16]：

(1)利用三次樣條插值法順序連接局部最大值(最小值)點，獲得上(下)包絡線。

(2)定義上下包絡線的均值m(t)。

(3)提取瞬時局部振蕩h(t)=x(t)-m(t) 值。

(4)重復步驟(1)～(3)，直到m(t)接近于零，則h(t)是一個IMF，記為c(t)。

(5)計算余項r(t)=x(t)-c(t)。

(6)余項r(t)作為x(t)，重復步驟(1)～(5)產生下一個IMF和余項，直到余項成為單調函數。

因此，原始信號x(t)可以重構為：

(1)

式中：ci(t)為第i個IMF；rn(t)為第n個余項。

2.2EEMD方法

由于EMD算法是利用局部極值構建包絡線的，信號的中斷會引起模態混疊，因此學者們提出了EEMD算法以解決此問題。EEMD算法[15]描述如下：

(1)令xi(n)=x(n)+wi(n)，式中wi(n)(i=1，…，l)為不同的高斯白噪聲。

2.3CEEMD方法

在CEEMD算法中，在原始信號中加入成對的白噪聲產生兩個IMFs集合，因此可以推出兩個信號和白噪聲組合的混合集合[16]：

(2)

其中：S是原始數據；N是添加的白噪聲；M1是原始數據和正噪聲的和；M2是原始數據和負噪聲的和。最終的IMF是IMFs和正負噪聲的集合。由于成對的噪聲可以有效消除最終白噪聲余項，因此CEEMD可以節省計算時間。

2.4CEEMDAN方法

為了比較EMD、EEMD、CEEMD和CEEMDAN 4種方法的分解效果，將理想心電(ECG)信號加上5 dB噪聲作為仿真信號，如圖1所示。利用上述4種方法進行分解，圖2和表1分別給出了4種方法分解后的重構誤差和標準差值，由仿真實驗可知，EMD和CEEMDAN的重構誤差及標準差值要小于EEMD和CEEMD，但是由于EMD方法具有模態混疊的缺點，因此CEEMDAN方法更適合于非穩定信號的分解。

圖1　ECG信號波形

(a)EMD

(b)EEMD

(c)CEEMD

(d)CEEMDAN

METHODEMDEEMDCEEMDCEEMDANSD6.0657e-150.3661.611e-146.0811e-15

3　改進的CEEMDAN方法

3.1曲線的曲率

曲率是指在一條曲線或不同曲線上不同點的彎曲程度[18]。每條曲線的曲率不同，曲率本身是曲率半徑的倒數。對于一般的函數情況，y=f(x)，曲率k為：

(3)

對于一條給定的曲線f(t)=(x(t)，y(t))，曲率k為：

(4)

根據實驗結果可知，曲率可以準確揭示不同頻率曲線的波峰或波谷的屬性和位置(如圖3所示)。另外，每個曲線都有其獨有的曲率。因此，利用EMD方法及其改進方法分解后的模態可以看作是具有不同頻率的不同曲線，即每個IMF的曲率不同。進而可以利用曲率展示不同模態的內在屬性，并用于相關模態選擇。

(b)信號曲率(b)Curvature of signal

3.2離散弗雷歇距離

弗雷歇距離(Frechet Distance， FD)利用兩個目標的路徑以及兩條曲線上所有離散點的距離，測量兩條曲線的相似度。弗雷歇距離最著名的一個例子是一個人和一個狗之間由一條狗繩連接，各自沿著兩個路徑行走，假設人和狗可以以不同的速度行走，但是不可以向后行走，弗雷歇距離就是人和狗從起點走向終點的相應路徑中狗繩的最短長度。這種測量方法直觀，而且與豪斯多夫距離等其他相似度測量方法相比，其可以更好地代表折線曲線的相似度。長度分別為m和n的折線曲線P和Q之間的弗雷歇距離可定義為[19]：

(5)

式中：Ψm,n=M([1:m+n],[0:m])×M([1:m+n]，[0:n])，[a:b]={a,a+1,…,b}，對于任何兩個整數a,ba≤b。

3.3模態選擇方法

對于原始信號y(t)=x(t)+n(t)，x(t)為有用信號，n(t)為噪聲信號。測量信號y(t)可以被分解為n個模態(IMF1，IMF2，…，IMFn)和一個余項，所有的模態按照頻率從高到低排列。一般而言，噪聲存在于前幾個IMFs中，純凈信號存在于后幾個IMFs中，因此在噪聲信號和有用信號之間會存在一個臨界模態將信號分解為：

(6)

本文中所提改進方法的主要思想是找出可以區分噪聲信號模態和純凈信號模態的相關模態序號K。由于所有的模態都有其固有的特性，所以每個模態的曲率曲線均不相同。模態的曲率可以記為C_IMFs，信號可以分解為

(7)

式中：C_IMFi是每個模態的曲率，n是模態的數量。有用信號的C_IMFs曲線波形與原始信號的C_IMFs曲線波形相似，與此相反，噪聲信號的C_IMFs 曲線波形與原始信號的C_IMFs曲線波形不同。因此，可以用曲線波形的突變來區分噪聲模態和信號模態，此時K值的求取很關鍵。通過計算各個C_IMFs與原始信號曲率曲線的離散弗雷歇距離可以獲得K值。具體步驟如下：

(1)利用EMD、EEMD、CEEMD和CEEMDAN 4種方法分解原始信號，獲得一系列IMFs。

(2)計算原始信號和所有模態的曲率，獲得C_IMFs。

(3)計算原始信號曲率和每個C_IMFs之間的離散弗雷歇距離，獲得FDs。

FDj={CIMF0，CIMFj}(j=1，…，n)，

(8)

式中CIMF0為原始信號y(t)曲率。

(4)相關模態序號K值由FDj確定：

K=argmin{FDj}.

(9)

(5)原始信號重構為：

(10)

3.4算法應用

將改進算法應用于真實信號y(t)，y(t)為由兩個不同頻率組成的周期信號，

y(t)=cos(2πf1t)+sin(2πf2t)，

(11)

式中：f1=4 Hz，f2=6 Hz，數據長度為1 024，在信號中加入信噪比為5 dB的高斯白噪聲(如圖4所示)。根據信號y(t)可知，第7個和第8個模態為有用信號模態(如圖5所示)。不同模態的曲率曲線如圖6所示，通過求取每個C_IMF與原始信號曲率曲線的離散弗雷歇距離可以獲得K值(如圖7所示)。純凈信號和去噪后的重構信號如圖8所示(彩圖見期刊電子版)。

圖4　原始信號波形

圖5　IMFs波形

圖6　C_IMFs波形

圖7　C_IMFs的離散弗雷歇距離

圖8　去噪后重構信號和純凈信號

4　結果與討論

4.1仿真信號實驗結果

對信號y(t)去噪后進行重構處理，對比EMD-CURVATURE-FD， EEMD-CURVATURE-FD， CEEMD-CURVATURE-FD，CEEMDAN-CURVATURE-FD 4種方法的結果(如圖9所示彩圖見期刊電子版)。輸入信噪比(SNFin)為-15～15 dB，間隔為3 dB。輸出信噪比的定義為：

(12)

(13)

圖9　4種方法輸出信噪比對比結果

表2表明，所提出的基于CEEMDAN的改進算法的SD值最小，即重構誤差最小。結合SNR的對比結果可知CEEMDAN-CURVATURE-FD方法在信號去噪及重構中具有較強的魯棒性。

表2　4種基于EMD方法的信號重構效果

前文中實驗表明CEEMDAN方法優于EMD、EEMD及CEEMD方法，更適合于非穩定信號的分解。針對信號y(t)，利用參考文獻[20](CEEMDAN-MI)和[21](CEEMDAN-NMHD)中提出的去噪方法與本文的基于CEEMDAN的改進算法進行對比分析，去噪重構后的SD值及SNR如表3所示，由結果可知，CEEMDAN-CURVATURE-FD方法去噪后的SD值為0.096 4，SNR為21.170 0 dB，優于參考文獻[20]和[21]中提出的去噪方法處理后的SD值(0.146 2，0.113 7)及SNR(17.083 6 dB，18.892 4 dB)，即重構誤差小，去噪能力強。

表3　不同方法的去噪效果(SNRin=5 dB)

4.2實測信號實驗結果

在近紅外無創血糖檢測實驗中，用葡萄糖溶液暫代血糖溶液，本實驗中所有葡萄糖溶液均為在同一條件下統一配置的連續均勻分布的液體，濃度范圍為50～1 000 mg/dL。所有近紅外光譜實驗數據由傅里葉光譜儀采集，滿足朗伯比爾定律測量原理。

近紅外光譜信號屬于非穩態信號，可以采用基于CEEMDAN的改進算法去除信號中的干擾噪聲，取500 mg/dL葡萄糖溶液作為研究對象，以其光譜數據去噪結果為例，如圖10所示(彩圖見期刊電子版)。

圖10　原始信號及去噪結果(500 mg/dL)

為了驗證上述算法的有效性，引入SNR和SD作為評價指標。由于用含噪信號代替公式(12)和公式(13)中的純凈信號y(n)，所以評價的衡量標準與仿真信號的評價標準相反，即SNR值越小，SD值越大，算法去噪效果越好。圖11(彩圖見期刊電子版)所示為通過不同方法進行光譜信號去噪處理后的SNR值和SD值對比圖，由圖可知，CEEMDAN-CURVATURE-FD方法優于其它3種方法。

圖11　光譜信號的SNR和SD(500 mg/dL)

為了驗證所提方法對光譜數據處理的普遍性，任意抽取10組樣本的光譜數據，利用4種方法分別對其進行去噪處理，SNR和SD值結果如表4和表5所示(表中Method1代表EMD-CURVATURE-FD方法，Method2代表EEMD-CURVATURE-FD方法，Method3代表CEEMD-CURVATURE-FD方法，Method4代表CEEMDAN-CURVATURE-FD方法)。由實驗結果可知，CEEMDAN-CURVATURE-FD方法適用于近紅外光譜數據的去噪處理。

表4不同濃度樣本光譜數據經4種方法處理后的標準差值

Tab.4　　　　SD of spectra data for sample with different concentrations by four methods　(dB)

表5不同濃度樣本光譜數據經4種方法處理后的信噪比

Tab.5　　　　SNR of spectra data for sample with different concentrations by four methods　(dB)

5　結　論

本文提出了一種應用于近紅外光譜的自適應去噪算法。首先利用EMD及其3種改進方法對仿真信號和由傅里葉光譜儀采集的葡萄糖溶液的光譜數據進行分解，然后求得原始信號及分解后IMFs的曲率曲線，并根據原始信號曲率和分解后IMFs曲率間的離散弗雷歇距離選擇相關模態，最后進行信號重構。將SNR和SD值作為評價指標。實驗結果顯示，對于近紅外光譜信號，CEEMDAN-CURVATURE-FD方法處理后的SD值為0.179 4，SNR值為19.117 5 dB，結果優于其它3種方法，該方法適用于紅外光譜數據的去噪。

[1]do AMARAL C E F， WOLF B. Current development in non-invasive glucose monitoring [J].Medical.Engineering&Physics， 2008，30(5)：541-549.

[2]DU NCAN A， HANNIGAN J， FREEBORN S S，etal..A portable non-invasive blood glucose monitor [C].The8thInternationalConferenceonSolid-StateSensorsandActuators， 1995， 2：455-458.

[3]BLANCO J R， FERRERO F J， CAMPO J C，etal.. Design of a low-cost portable potentiostat for amperometric biosensors [C]. 2006IEEEInstrumentationandMeasurementTechnologyConferenceProceedings， 2006： 690-694.

[4]UNNIKRISHNA K A， HEMACHANDRAN D， ABHISHEK T K. A survey on non-invasive blood glucose monitoring using NIR [C].2013InternationalConferenceonCommunicationsandSignalProcessing(ICCSP)， 2013：1069-1072.

[5]周玲芳， 陳菲. 基于斜率差值的自適應圖像椒鹽噪聲濾波算法[J]. 液晶與顯示， 2015， 30(4)：695-700.

ZHOU L F， CHEN F. Adaptive slope difference algorithm for filtering salt and pepper noise in image [J].ChineseJournalofLiquidCrystalsandDisplays， 2015，30(4)： 695-700.(in Chinese)

[6]董雪，林志賢， 郭太良. 基于LoG算子改進的自適應閾值小波去噪算法[J]. 液晶與顯示， 2014， 29(2)：275-280.DONG X， LIN ZH X， GUO T L. Improved self-adaptive threshold wavelet denoising analysis based on LoG operator [J].ChineseJournalofLiquidCrystalsandDisplays， 2014，29(2)： 275-280.(in Chinese)

[7]李權， 趙勛杰， 彭青艷，等. 基于主成分分析法的窗口自適應粒子濾波算法[J]. 紅外與激光工程， 2014， 43(10)：3474-3479.

LI Q， ZHAO X J， PENG Q Y，etal.. Windows adaptive particle filter algorithm based on principal component analysis [J].InfraredandLaserEngineering， 2014， 43(10)：3474-3479.(in Chinese)

[8]顧有林， 葉應流， 曹光華，等. EMD和小波變換在低可探測目標檢測中的應用 [J]. 紅外與激光工程， 2015， 44(11)：3494-3499.

GU Y L， YE Y L， CAO G H，etal.. Application of EMD and wavelet transform in low detectable targets detection [J].InfraredandLaserEngineering， 2015，44(11)：3494-3499.(in Chinese)

[9]HUANG N E. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis [J].Proc.R.Soc.Lond.A，1998， 454：903-995.

[10]李欣，梅德慶，陳子辰. 基于經驗模態分解和希爾伯特-黃變換的精密孔鏜削顫振特征提取[J].光學 精密工程， 2011，19(6)： 1291-1297.

LI X，MEI D Q，CHEN Z CH.Feature extraction of chatter for precision hole boring processing based on EMD and HHT [J].Opt.PrecisionEng.， 2011，19(6)： 1291-1297.(in Chinese)

[11]REN Y， SUGANTHAN P N， SRIKANTH N. A comparative study of empirical mode decomposition-based short-term wind speed forecasting methods [J].IEEETransactionsonSustainableEnergy2015， 6(1)： 236-244.

[12]魯麗，顏國正，趙凱，等. 基于集合經驗模態分解的人體結腸動力分析[J].光學 精密工程， 2015，23(6)： 1850-1856.LU L， YAN G ZH， ZHAO K，etal.. Analysis of human colonic motility using EEMD [J].Opt.PrecisionEng.， 2015，23(6)： 1850-1856.(in Chinese)

[13]GAN Y， SUI L， WU J，etal.. An EMD threshold denoising method for inertial sensors [J].Measurement，2014， 49(1)， 34-41.

[14]羅玉昆，羅詩途，羅飛路，等. 激光超聲信號去噪的經驗模態分解實現及改進[J].光學 精密工程， 2013，21(2)： 479-487.LUO Y K， LUO SH T， LUO F L，etal.. Realization and improvement of laser ultrasonic signal denoising based on empirical mode decomposition[J].Opt.PrecisionEng.， 2013，21(2)： 479-487.(in Chinese)

[15]MARIA E T， MARCELO A C， GASTON S，etal.. A complete ensemble empirical mode decomposition with adaptive noise [C]. 2011IEEEInternationalConferenceonAcoustics，SpeechandSignalProcessing(ICASSP)， 2011： 4144-4147.

[16]JIA R Y， JIANN S S. Complementary ensemble empirical mode decomposition： A novel noise enhanced data analysis method [J].AdvancesinAdaptiveDataAnalysis， 2010， 2(2)：135-156.

[17]ANNE H H， PIERRE A， GUILLAUME M，Analysis of laser speckle contrast images variability using a novel empirical mode decomposition： comparison of results with laser Doppler flowmetry signal variability [J].IEEETransactionsonMedicalImaging， 2015， 34(2)：618-626.

[18]COOLIDGE J L. The unsatisfactory story of curvature [J].TheAmericanMathematicalMonthly， 1952， 59(6)： 375-379.

[19]R S， KARTHIK K， CHIRANJIB B. Frechet distance based Approach for searching online handwritten documents [C].NinthInternationalConferenceonDocumentAnalysisandRecognition，(ICDAR)，2007.

[20]HAN L， LI C， LIU H. Feature extraction method of rolling bearing fault signal based on EEMD and cloud model characteristic entropy [J].Entropy， 2015， 17(10)：6683-6697.

[21]LI C， ZHAN L. A hybrid filtering method based on a novel empirical mode decomposition for friction signals [J].MeasurementScience&Technology， 2015， 26(12)：125003.

李曉莉(1984-)，女，黑龍江哈爾濱人，博士研究生，2008年、2011年于東北農業大學分別獲得學士、碩士學位，主要從事無創血糖檢測及胰島素注射劑量調節研究。E-mail： xiaoli72460@163.com

導師簡介：

李成偉(1963-)，男，黑龍江哈爾濱人，博士，教授，1985年、1988年、2000年于哈爾濱工業大學分別獲得學士、碩士、博士學位，主要從事生物醫學工程、智能控制技術研究。E-mail： lcw@hit.edu.cn

(版權所有未經許可不得轉載)

Application of improved complete ensemble empirical mode decomposition with adaptive noise in spectral signal denoising

LI Xiao-li， LI Cheng-wei*

(SchoolofElectricalEngineeringandAutomation，HarbinInstituteofTechnology，Harbin150001，China)

*Correspondingauthor，E-mail：lcw@hit.edu.cn

As the accuracy and stability of a blood glucose level model is affected by the noise in near infrared non-invasive blood glucose detection process, an improved complete ensemble empirical mode decomposition method with adaptive noise was proposed for denoising of near infrared spectroscopy signals. Meanwhile, a mode selection method based on Frechet distance combining with the feature of curve curvature was proposed for the selection of Intrinsic Mode Functions(IMFs). Firstly. the complete ensemble empirical mode decomposition method with adaptive noise was introduced in the denoising processing of near infrared spectroscopy, and the basic principles and concrete realization processes of empirical mode decomposition, ensemble empirical mode decomposition, complementary ensemble empirical mode decomposition and the complete ensemble empirical mode decomposition based on adaptive noise were described. Then, an improved complete ensemble empirical mode decomposition method with adaptive noise based on curvature and discrete Frechet distance was applied in denoising for simulation signals and spectral signals， and their standard deviation and the Signal to Noise Ratio(SNR) were taken as the evaluation indexes. The simulation and experimental results show that the standard deviation of the improved method based on curvature and discrete Frechet distance in the near infrared spectral signal is 0.179 4, and the SNR is 19.117 5 dB, which extracts useful information, realizes the separation of signal and noise, and improves the quality of reconstructed signals. The proposed method has a good adaptability to effectively identify and separate the signal and noise components.

non-invasive blood glucose detection; near infrared spectroscopy; signal denoising; Complete Ensemble Empirical Mode Decomposition with Adaptive Noise(CEEMDAN); curvature; discrete Frechet distance

2016-03-10；

2016-04-25.

哈爾濱市科技創新人才專項資金資助項目(No.2014RFXXJ065)；哈爾濱工業大學理工醫交叉學科基礎研究培育計劃資助項目(No.HIT.IBRSEM.201307)

1004-924X(2016)07-1754-09

O657.33；R587.1

Adoi:10.3788/OPE.20162407.1754