翟哲 李偉凱 李長凱 裴玉

摘要：基于每個IMF自關(guān)聯(lián)函數(shù)的特征，提出了一種新的EMD去噪方法。以檢測苗期玉米葉片葉綠素含量為例，首先對原光譜信號采用SNV+Detrending方法進行預(yù)處理，然后利用該方法對預(yù)處理后的信號進行去噪，并與小波去噪方法和EMD融合小波去噪方法進行對比，最后應(yīng)用偏最小二乘回歸方法進行校正模型的建立。結(jié)果表明：將該方法應(yīng)用到實際近紅外光譜信號去噪中，其預(yù)測集決定系數(shù)（r2）達到0.984，殘差均方根RMSE為0.075，證明該方法在近紅外光譜處理過程中具有很好的去噪效果，建立的校正模型具有較高的魯棒性和推廣性。

關(guān)鍵詞：經(jīng)驗?zāi)B(tài)分解;自適應(yīng);近紅外光譜;葉綠素

中圖分類號： S126 文獻標志碼： A

文章編號：1002-1302（2015）04-0392-05

收稿日期：2014-05-25

基金項目：黑龍江八一農(nóng)墾大學(xué)研究生創(chuàng)新科研項目（編號：YJSCX2013-16BYND）。

作者簡介：翟 哲（1988—），男，黑龍江慶安人，碩士研究生，主要從事數(shù)據(jù)處理研究。E-mail：zhai_zhe@163.com。

通信作者：李偉凱，博士，教授，博士生導(dǎo)師，主要從事光電檢測研究。Tel：（0459）6819009;E-mail：bynd@263.net.cn。

隨著計算機科學(xué)技術(shù)的飛速發(fā)展，近紅外光譜技術(shù)在儀器分析領(lǐng)域受到了有關(guān)專家的高度重視[1]。近紅外光譜技術(shù)之所以能迅速發(fā)展并在各個行業(yè)都有廣泛的應(yīng)用，是因為它有很多的優(yōu)越性[2]。但是由于外界環(huán)境的影響，近紅外光譜儀所采集到的光譜信號，除了包含自身信息外，在測量中還不可避免地得到許多無關(guān)的噪聲信號[3]。因此在使用化學(xué)計量方法建立校正模型時，消除光譜數(shù)據(jù)無關(guān)噪聲在光譜數(shù)據(jù)分析中變得十分關(guān)鍵和必要。本研究將改進經(jīng)驗?zāi)B(tài)分解（empirical mode decomposition，EMD）方法引入到近紅外光譜信號去噪中，旨在探索近紅外光譜信號去噪的新方法。

1 基本思想

EMD方法就是把1個非線性非平穩(wěn)的信號分解為有限個本征模函數(shù)（intrinsic mode function，IMF）分量和1個趨勢項[4]，原始信號x（t）可表示為：

x（t）=∑ni=1ci（t）+rn（t）。

（1）

式中：ci（t）為第i個IMF分量;rn（t）為篩選到最后剩下的趨勢信息;t為時間;n為分解的IMF個數(shù)。各個IMF分量都代表信號從高頻到低頻的分量，通常情況下階數(shù)較小的IMF代表高頻分量和噪聲信號，階數(shù)較大的IMF代表低頻分量，受噪聲影響較小。

EMD其實就是把信號的極值特征尺度作為度量而進行篩選的過程，信號從最小特征尺度實行篩選，因此得到周期最短的IMF[5-6]。之后再進行逐層篩選，最終獲得周期尺度漸次增大的多個IMF，此過程亦顯示出了多分辨辨識的濾波全程。這個途徑是依據(jù)信號分解的IMF分量建構(gòu)濾波函數(shù)，所以能極大保持信號固有的非線性、非平穩(wěn)特征。

假設(shè)x（t）為含噪信號，則經(jīng)過EMD算法的分解后其高通濾波表達如下：

xhp（t）=∑ki=1ci（t），（1

（2）

帶通濾波表達如下：

xdp（t）=∑ki=hci（t），（1

（3）

低通濾波表達如下：

xlp（t）=∑ni=kci（t）+rn（t），（1

（4）

但是，噪聲同信號在IMF分量內(nèi)疊混，可以用EMD閾值去噪的方法[7]。其中在閾值選擇方面，根據(jù)Donoho等給出的小波去噪中的閾值[8]，其中小波去噪軟閾值為：

ysoft（t）=sng[x（t）]·[|x（t）|-δ]，|x（t）|>δ

0，|x（t）|≤δ。

（5）

式中：δ為通用閾值。

在第j層選取δ=σj2lnN。

式中：N為信號的長度;σj為噪聲在第j層的標準差，可利用σj=media/0.674 5進行估計，media為第j層上小波系數(shù)的絕對中值。

因此，將這個等式略加改變后運用至IMF去噪上，即獲公式：

imf ′j=sng[imfj（i）]·[|imfj（i）|-δ]，|imfj（i）|>δj

0，|imfj（i）|≤δj。

（6）

式中：δj為第j個分量imfj的閾值，由下式可以計算：

δj=σj2lnN=media[abs（imfj）]0.674 5×2lnN。

（7）

式中：media為第j個imfj分量上的絕對中值。

2 基于能量最小原則的EMD自適應(yīng)去噪方法

信號經(jīng)由EMD算法進行分解后，可獲得可數(shù)的頻率由高至低的IMF，它們中階數(shù)較低的IMF與信號的高頻成分相對應(yīng)，普遍涵蓋了信號尖銳部分或噪聲;階數(shù)較高的IMF與信號的低頻成分相對應(yīng)，低頻成分中噪聲能量較弱。EMD去噪的首要內(nèi)容即為針對大部分含噪信號，它的主要能量累積在低頻段范圍內(nèi)，越往高頻范圍，它所蘊含的能量便越小。假設(shè)信號區(qū)和噪聲區(qū)的分界點為k，因此肯定有1個IMFk分量，讓這個分量之后的IMF里的信號成為主導(dǎo)模態(tài)，而在這之前，IMF里的噪聲則是主要模態(tài)，因此原始的基于EMD進行去噪的目的即為發(fā)現(xiàn)IMFk。

Boudraa等首先提出連續(xù)均方誤差（consecutive mean square error，CMSE）準則的EMD去噪法，這種方法利用了IMF能量的全局極小值位置作為噪聲起主導(dǎo)作用與信號起主導(dǎo)作用的分界點[9]。文獻[9]和文獻[10]都采用了連續(xù)均方誤差準則這種方法，對信號成分發(fā)揮主導(dǎo)功能模態(tài)和噪聲成分發(fā)揮主導(dǎo)功能模態(tài)進行了區(qū)分，利用反映信號主要結(jié)構(gòu)的模態(tài)對信號進行部分重構(gòu)以實現(xiàn)去噪。文獻[9]將全部IMF能量的最小值點對應(yīng)的模態(tài)當(dāng)作信號成分與噪聲成分的分界點，但是判斷準則可能存在偏差，主要體現(xiàn)在2個方面：首先，若CMSE在全局極小值前面有局部極小值，則局部極小值對應(yīng)的位置應(yīng)加1;其次，如果不存在局部極小值，則全局極小值對應(yīng)的位置應(yīng)加1。因此文獻[10]在文獻[9]的基礎(chǔ)上做了一些改進，采用首個局部極小值點對應(yīng)的模態(tài)作為分界點，若在全部極小值前存在著局部極小值時，則有：

k=argfirstlocal1≤j≤n-1min[CMSE（xj，xj+1）]+1;

（8）

否則：

k=argmin[CASE（xj，xj+1）]+11≤j≤n-1。

（9）

式中：firstlocal表示首個局部極小值。

但是在試驗中不難發(fā)現(xiàn)，文獻[10]中提出的改進方法與其采用的連續(xù)均方誤差有關(guān)，信號模態(tài)和噪聲模態(tài)分選準則仍然存在一定的偏差。對此文獻[11]對基于各IMF與所給信號相關(guān)系數(shù)的大小進行了修正：它考慮到文獻[10]的分選準測，信號在低信噪比情況下存在的問題，因此選擇各階模態(tài)與原始信號之間的相關(guān)系數(shù)作為噪聲模態(tài)與信號模態(tài)的分選準則來代替能量法的連續(xù)均方誤差準則，從而減輕因信號模態(tài)能量較小而被舍棄的概率。各階模態(tài)函數(shù)與原始信號的相關(guān)系數(shù)定義如下：

R[x（ti），IMFj（ti）]=cov[x（ti），IMFj（ti）]cov x（ti）cov IMFj（ti）。

（10）

式中：cov（·）代表協(xié)方差，噪聲模態(tài)與信號模態(tài)的分界點k由公式（11）確定：

k=arg firstlocal1≤j≤nmin{R[x（ti），IMFj（ti）]}+1。

（11）

主要步驟：首先，對原信號進行經(jīng)驗?zāi)B(tài)分解;其次，以各階IMF與原信號的相關(guān)系數(shù)曲線中第1個極小值點的位置作為標準，找到噪聲主導(dǎo)模態(tài)和信號主導(dǎo)模態(tài)的分界點;最后，將IMF中信號起主導(dǎo)作用的模態(tài)成分進行重構(gòu)來實現(xiàn)去噪。

這些方法在相當(dāng)程度上都解決了基于EMD算法進行濾波去噪的自適應(yīng)濾波器終結(jié)階數(shù)k不好界定的問題。如果信噪比較高，此判別法較為有效，但是在噪聲污染太重的低信噪比前提下，這個辦法的功能并不是很穩(wěn)固。這是由于低信噪比時，可用信號能量極少，由EMD算法分解而來的IMF模態(tài)中，有的IMF能量差別較小，即便是某個IMF的能量達到局部最小值，卻也并不一定對噪聲起主要作用。因此，選取此IMF后面的IMF來實現(xiàn)重構(gòu)時，將可能導(dǎo)致有用信息遺失。更重要的是，在某些特殊條件下，難以發(fā)現(xiàn)IMF能量全部極小值，此時這種方法完全行不通。同時，噪聲和有用信號在IMF分量里有時有混疊行為，而如果依然用原方法除去一些IMF，就會誤刪某些有用信號的高頻成分，破壞信號的完整性，從而達不到良好的去噪效果。

基于EMD去噪的主要目的就是尋找信號主導(dǎo)模態(tài)與噪聲主導(dǎo)模態(tài)的分界點k，但是鑒于含有高斯白噪聲的信號應(yīng)用EMD算法進行分解時，由于EMD算法的自適應(yīng)性和強制信號的上下對稱性，使分解而來序數(shù)小的蘊含噪聲的IMF里，高斯白噪聲特征被毀壞，獲得非真正意義上的白噪聲分量。只是白噪聲統(tǒng)計特征仍近似留存，也就是其在零點處自關(guān)聯(lián)函數(shù)為最大，別處雖然不為0，但是相對零點處的自關(guān)聯(lián)函數(shù)值的關(guān)聯(lián)性已很小，極速衰減。對于普通信號，它們的自關(guān)聯(lián)函數(shù)在零點處獲得最大值而在其他地方則并不一定為0，它會隨時間差的變化而改變。鑒于這些特征，本研究依據(jù)每個IMF自關(guān)聯(lián)函數(shù)特征，提出基于能量最小原則的EMD自適應(yīng)去噪方法，我們在斷定噪聲發(fā)揮主要功能的IMF模態(tài)中，它的構(gòu)成成分除噪聲外，還包含少數(shù)有用信號的高頻部分。至此，對于噪聲模態(tài)分量篩選的所有閾值，運用去噪法把想除去的IMF分量實行濾波，再將全部經(jīng)過處置和未處置的IMF分量實行重構(gòu)，從而獲得去噪的信號。其算法步驟為：（1）對含噪信號應(yīng)用EMD算法進行分解，得到n個IMF分量;（2）按照公式10所示，分別計算各階IMF自相關(guān)系數(shù);（3）依照能量的積累高低，由公式11斷定信號主導(dǎo)模態(tài)和噪聲主導(dǎo)模態(tài)分界點k;（4）把噪聲主導(dǎo)模態(tài)的IMF對細節(jié)系數(shù)取1/3進行抑制，獲取去噪后的IMF分量;（5）信號主導(dǎo)模態(tài)IMF分量與去噪后的IMF分量相加，重構(gòu)原信號。

3 結(jié)果與分析

3.1 光譜采集與預(yù)處理

使用北京瑞利分析儀器公司生產(chǎn)的WQF-600N傅立葉變換近紅外光譜儀（波數(shù)范圍3 300～10 000 cm-1，分辨率優(yōu)于4 cm-1）來進行近紅外光譜數(shù)據(jù)的采集，儀器自帶MainFOTS軟件，可將采集到的光譜數(shù)據(jù)保存到計算機中。試驗選取樣品為同一生長期各種不同外觀、形狀的苗期玉米葉片，共計42張，42張苗期玉米葉片的原始光譜結(jié)果見圖1-a。

從圖1-a可以看出，所采集的原始近紅外光譜信號中包含了樣品組成的特征信息，但是基線漂移嚴重，并且還包含由溫度、濕度、時間、樣品背景等影響而造成的許多與信號無關(guān)的噪聲，因此分析偏差較大。所以本研究先對原始光譜數(shù)據(jù)實行預(yù)處理的規(guī)范化工作，然后對規(guī)范化后的光譜數(shù)據(jù)進行去噪處理。

比較常用的數(shù)據(jù)規(guī)范化方法有變量標準化（SNV）、趨勢變換法（Detrending）、SNV+Detrending、中心化（MC）、歸一化等方法，通過對比驗證、比較效果，試驗已經(jīng)證明SNV+Detrending 預(yù)處理方法優(yōu)于其他預(yù)處理方法。因此本研究先對原始光譜先進行數(shù)據(jù)規(guī)范化，這樣就可以消除變異程度從而實現(xiàn)數(shù)據(jù)規(guī)范。圖1-b為采用SNV+Detrending方法預(yù)處理后的光譜結(jié)果。

從圖1-b可以看出：實行規(guī)范化后的光譜基本除去了部分線性或接近線性的背景噪聲對目標光譜的影響;但是經(jīng)過處理后光譜的噪聲并沒有消除，因此還要對預(yù)處理后的光譜進行消噪處理。

3.2 光譜信號去噪

為了對基于能量最小原則的EMD自適應(yīng)去噪效果進行比較，本研究對全部光譜進行去噪。下面采用此方法對10號苗期玉米葉片樣品進行EMD自適應(yīng)去噪處理，10號苗期玉米葉片樣品預(yù)處理后的光譜經(jīng)過EMD算法分解后，得到6階imf（imf1～imf6）分量和1個殘余項（imf7），每個IMF都有不同的振幅和頻率，分解按照頻率從高到低自適應(yīng)地進行，然后自適應(yīng)地選出分界點，進行抑制和重構(gòu)。從去噪后的光譜圖（圖2-c）可以看出，光譜非常光滑，噪聲基本上得到了消除。其中10號樣品預(yù)處理后的光譜圖、6階模態(tài)函數(shù)、去噪后的光譜結(jié)果分別如圖2-a、圖2-b、圖2-c所示。

采用同樣的方法對其他41組樣品預(yù)處理后的光譜圖做了去噪處理，為了作對比，給出了10號樣品預(yù)處理后光譜的小波去噪、EMD融合小波去噪、基于能量最小原則的EMD自適應(yīng)去噪方法的去噪效果圖。從圖3可以看出，基于能量最小原則的EMD自適應(yīng)去噪方法去噪效果最好，幾乎完全消除了噪聲的干擾，在去除大部分噪聲的同時對光譜的峰形沒有太大影響，而且反映原始信號的特征尖峰點得到很好的保留。小波去噪、EMD融合小波去噪的去噪效果也較好，對峰形也沒有太大影響，但母小波的選擇對去噪效果影響很大。相比之下，基于能量最小原則的EMD自適應(yīng)去噪方法是最有效的，光譜平滑性增強，且對光譜的峰形沒有太大影響。與其他2種方法相比，基于能量最小原則的EMD自適應(yīng)去噪方法具

有完全自適應(yīng)性，這就避免了選擇母小波的麻煩，增加了分析結(jié)果的可靠性。

通過對規(guī)范化后的光譜信號進行去噪處理，可以濾除其含有的無關(guān)噪聲信號，從而降低噪聲對光譜分析的干擾，充分提取有效信息，有利于建立最佳數(shù)學(xué)模型，提高分析的準確度。

3.3 玉米葉片葉綠素含量檢測

從編號的42個苗期玉米葉片樣品中隨機抽取6號、9號、10號、17號、26號、32號、35號、37號、39號、40號這10個樣品作為預(yù)測集，用于檢驗方程的可靠程度，其余32個樣品作為校正集，用來進行校正模型的建立，下一步建立校正模型的校正樣本和預(yù)測樣本散點圖。

采用OPUS軟件建立葉綠素偏最小二乘回歸模型，圖4、圖5、圖6為用偏最小二乘回歸法[12]所得到的苗期玉米葉片葉綠素含量校正樣本和預(yù)測樣本的散點圖，其中圖4、圖5、圖6分別是使用小波去噪、EMD融合小波去噪、基于能量最小原則的EMD自適應(yīng)去噪的散點圖，可以看出苗期玉米葉片預(yù)測集和實測集的相關(guān)性。表1列出了3種預(yù)處理方法苗期玉米葉片葉綠素含量預(yù)測值和實測值的比較結(jié)果，可以看出，經(jīng)過基于能量最小原則的EMD自適應(yīng)去噪后，光譜分析的準確性明顯提高。結(jié)果表明，采用基于能量最小原則的EMD自適

應(yīng)去噪方法建立模型對預(yù)測樣品進行預(yù)測的準確性明顯得到了改善。

應(yīng)用決定系數(shù)r2和殘差均方根RMSE作為評價偏最小二乘回歸模型優(yōu)劣的標準。經(jīng)基于能量最小原則的EMD自適應(yīng)去噪后，預(yù)測值與實測值之間的決定系數(shù)r2為0.984，殘差均方根RMSE為0.075。試驗數(shù)據(jù)說明，該模型對苗期玉

米葉片葉綠素含量有較好的預(yù)測效果，提高了葉綠素近紅外檢測的魯棒性。

4 結(jié)論

本研究以苗期玉米葉片近紅外光譜為研究對象，將EMD去噪新方法應(yīng)用到葉綠素近紅外光譜分析中。首先對原始光譜圖用SNV+Detrending方法進行規(guī)范化預(yù)處理，實現(xiàn)數(shù)據(jù)規(guī)范;其次將新方法應(yīng)用到規(guī)范化后的光譜去噪中，并與小波去噪方法和EMD融合小波去噪方法進行比較;最后采用偏最小二乘回歸方法進行了校正模型的建立，測量葉綠素含量。結(jié)果表明，采用基于能量最小原則的EMD自適應(yīng)去噪方法進行近紅外光譜數(shù)據(jù)去噪是可行的，此方法比采用傳統(tǒng)的去噪方法有更高的預(yù)測精度。EMD方法在紅外光譜處理中可以有效地消除噪聲的影響，提取光譜中的有效信息，在光譜數(shù)據(jù)處理中將會有更廣闊的應(yīng)用前景，從而為下一步研究打下良好的基礎(chǔ)。

參考文獻：

[1]高榮強，范世福. 現(xiàn)代近紅外光譜分析技術(shù)的原理及應(yīng)用[J]. 分析儀器，2002（3）：9-12.

[2]張小超，吳靜珠，徐 云. 近紅外光譜分析技術(shù)及其在現(xiàn)代農(nóng)業(yè)中的應(yīng)用[M]. 北京：電子工業(yè)出版社，2012.

表1 預(yù)測集樣品的預(yù)測結(jié)果

樣品序號 實測值

（mg/g）

小波 EMD融合小波 EMD自適應(yīng)

預(yù)測值

（mg/g） 絕對誤差

（mg/g） 預(yù)測值

（mg/g） 絕對誤差

（mg/g） 預(yù)測值

（mg/g） 絕對誤差

（mg/g）

6 2.174 1 2.383 4 0.209 3 2.285 4 0.111 3 2.083 7 -0.090 4

9 2.184 1 2.3729 0.188 8 2.280 4 0.096 3 2.252 8 0.068 7

10 1.706 1 1.562 8 -0.143 3 1.796 7 0.090 6 1.754 4 0.048 3

17 3.298 5 3.176 2 -0.122 3 3.217 3 -0.081 2 3.241 2 -0.057 3

26 3.253 6 3.132 8 -0.120 8 3.331 4 0.077 8 3.326 1 0.072 5

32 2.276 9 2.172 7 -0.104 2 2.417 4 0.140 5 2.159 5 -0.117 4

35 2.607 3 2.726 7 0.119 4 2.691 2 0.083 9 2.696 8 0.089 5

37 2.401 7 2.265 5 2.119 7 -0.145 8 2.146 7 -0.118 8 2.208 2

39 2.265 5 1.449 5 1.339 8 -0.109 7 1.380 7 -0.068 8 1.395 5

40 1.449 5 2.972 9 2.874 9 -0.09 8 2.872 6 -0.100 3 3.035 6

r2 0.942 0.971 0.984

RMSE 0.141 0.099 0.075

[3]李 勇，魏益民，王 鋒. 影響近紅外光譜分析結(jié)果準確性的因素[J]. 核農(nóng)學(xué)報，2005，19（3）：236-240.

[4]Huang N E，Shen Z，Long S R，et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Soeiety of London，Series A，1998，454：903-995.

[5]高忠劍波，蔡宗平. 運用Hilbert-Huang變換的振動測試分析系統(tǒng)[J]. 工業(yè)儀表與自動化裝置，2009（3）：88-90.

[6]楊志華，齊東旭，楊力華. 一種基于Hilbert-Huang變換的基音周期檢測新方法[J]. 計算機學(xué)報，2006，29（1）：106-115.

[7]邵忍平，曹精明，李永龍. 基于EMD小波閾值去噪和時頻分析的齒輪故障模式識別與診斷[J]. 振動與沖擊，2012，31（8）：96-101，106.

[8] Donoho D L，Johnstone L M，Kerkyacharian G，et al. Wavelet shrinkage：asymptopia？[J]. IEEE Transaction on Information，1995，41（3）：613-627.

[9]Boudraa A O，Cexus J C. EMD-based signal filtering[J]. IEEE Transactions on Instrumentation and Measurement，2007，56（6）：2196-2202.

[10]孫偉峰，彭玉華，許建華. 基于EMD的激光超聲信號去噪方法[J]. 山東大學(xué)學(xué)報：工學(xué)版，2008，38（5）：121-126.

[11]陳鳳林. 一種新的基于EMD模態(tài)相關(guān)的信號去噪方法[J]. 西華大學(xué)學(xué)報：自然科學(xué)版，2009，28（6）：20-24.

[12]張 銀，周孟然. 近紅外光譜分析技術(shù)的數(shù)據(jù)處理方法[J]. 紅外技術(shù)，2007，29（6）：345-348.