趙肖宇，方一鳴，譚 峰，佟 亮，翟 哲

1. 黑龍江八一農(nóng)墾大學(xué)信息技術(shù)學(xué)院，黑龍江 大慶 163319 2. 燕山大學(xué)電氣工程學(xué)院，河北 秦皇島 066004 3. 齊齊哈爾大學(xué)通信與電子工程學(xué)院，黑龍江 齊齊哈爾 161006 3. 中國(guó)林業(yè)科學(xué)研究院華北林業(yè)實(shí)驗(yàn)中心，北京 102300

EMD時(shí)頻分析拉曼光譜和近紅外光譜

趙肖宇1，方一鳴2，譚 峰1，佟 亮3，翟 哲4

1. 黑龍江八一農(nóng)墾大學(xué)信息技術(shù)學(xué)院，黑龍江 大慶 163319 2. 燕山大學(xué)電氣工程學(xué)院，河北 秦皇島 066004 3. 齊齊哈爾大學(xué)通信與電子工程學(xué)院，黑龍江 齊齊哈爾 161006 3. 中國(guó)林業(yè)科學(xué)研究院華北林業(yè)實(shí)驗(yàn)中心，北京 102300

用時(shí)頻方法分析拉曼光譜和近紅外光譜。經(jīng)驗(yàn)?zāi)B(tài)分解光譜成為特征模態(tài)分量，模態(tài)分量比重計(jì)算顯示拉曼光譜能量均勻分布于各個(gè)分量，而近紅外光譜的低階特征模態(tài)分量只承載了較少的原光譜有效信息。真實(shí)光譜和數(shù)值實(shí)驗(yàn)均顯示，經(jīng)驗(yàn)?zāi)B(tài)分解視拉曼光譜為調(diào)幅信號(hào)，具有高頻能量吸附特性； 視近紅外光譜為調(diào)頻信號(hào)，在一階特征模態(tài)分量中可以較好實(shí)現(xiàn)高頻窄帶解調(diào)。一階特征模態(tài)分量希爾伯特變換顯示，經(jīng)驗(yàn)?zāi)B(tài)分解拉曼光譜時(shí)易出現(xiàn)模態(tài)混疊現(xiàn)象。進(jìn)一步在時(shí)頻域分析玉米葉片近紅外光譜，經(jīng)驗(yàn)?zāi)B(tài)分解后截掉低能量的一、二階分量，用剩余特征模態(tài)分量重構(gòu)光譜信號(hào)，均方根誤差為1.001 1，相關(guān)系數(shù)為0.981 3，兩個(gè)指標(biāo)反映出重構(gòu)精度較高； 分解趨勢(shì)項(xiàng)表明在近紅外光波段，吸光度隨著波長(zhǎng)的減小呈現(xiàn)遞增趨勢(shì)； 特征模態(tài)分量的希爾伯特變換顯示，657 cm-1是堿脅迫光譜特有頻率，可作為堿脅迫光譜特征頻率來(lái)辨識(shí)。

時(shí)頻分析； 經(jīng)驗(yàn)?zāi)B(tài)分解； 希爾伯特變換； 近紅外光譜； 拉曼光譜

引 言

透射率(T)、反射率(R)或者吸光度(A)相對(duì)于波長(zhǎng)或波數(shù)的變化構(gòu)成了近紅外光譜； 拉曼光譜橫坐標(biāo)是散射光相對(duì)于入射光的波數(shù)差(拉曼位移)，縱坐標(biāo)是光子計(jì)數(shù)。事實(shí)上，根據(jù)波速、波長(zhǎng)與頻率關(guān)系，無(wú)論是拉曼還是近紅外光譜都是頻率譜圖，可以借助傅里葉反變換得到對(duì)應(yīng)時(shí)頻譜，那么現(xiàn)有時(shí)頻分析方法即可以用到光譜分析，操作細(xì)節(jié)有待進(jìn)一步研究。近紅外或拉曼光譜通過(guò)特征峰定性辨識(shí)某種分子結(jié)構(gòu)或定量分析該成分，特征峰定位過(guò)程與尋找時(shí)間函數(shù)變異點(diǎn)過(guò)程相似，所以從譜圖分析角度，可以直接將光譜視作時(shí)變函數(shù)處理，這樣除了信號(hào)原有“時(shí)間”參量，通過(guò)時(shí)頻方法處理，光譜還具備了“頻率”參數(shù)，以及可以設(shè)計(jì)時(shí)間和頻率聯(lián)合函數(shù)，用它們表達(dá)光譜更多方面的信息。

經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition，EMD)方法是美國(guó)國(guó)家宇航局美籍華人黃鍔(N. E. Huang)等于1998年創(chuàng)造性提出新型的信號(hào)時(shí)頻處理方法[1]，該方法一經(jīng)提出即得到了學(xué)術(shù)界的最廣泛的關(guān)注。2009年出現(xiàn)基于EMD方法分解拉曼光譜并結(jié)合小波算法對(duì)特征模態(tài)分量(intrinsic mode function，IMF)去噪的研究，文獻(xiàn)[2]對(duì)近紅外導(dǎo)數(shù)光譜經(jīng)驗(yàn)?zāi)B(tài)分解后針對(duì)低階IMF閾值濾波，也可以直接去掉近紅外原光譜中低階IMF分量以實(shí)現(xiàn)濾波[3]。目前已有EMD方法在光譜預(yù)處理方面的應(yīng)用，但就低階IMF處理是否可以實(shí)現(xiàn)去噪，時(shí)頻方法處理近紅外光譜和拉曼光譜的區(qū)別以及EMD方法是否可以用于光譜定性和定量分析，這些涉及到光譜IMF屬性的研究及其時(shí)頻分析，至今依然未見(jiàn)報(bào)道。本研究的主要目的就是通過(guò)數(shù)值實(shí)驗(yàn)的手段，研究光譜EMD的一階IMF時(shí)頻屬性，揭示EMD方法分解拉曼光譜和近紅外光譜的差異性，并試圖用光譜的希爾伯特變化(hilbert transform, HT)辨識(shí)光譜的特異性。

1 光譜采集

實(shí)驗(yàn)采用WQF-600N傅立葉變換近紅外光譜分析儀，光源電壓： 5 V，在10 000～3 500 cm-1光譜范圍內(nèi)以8 cm-1分辨率掃描32次，本底掃描64次。采用美國(guó)DeltaNu拉曼光譜儀Advandge532系列，激光激發(fā)波長(zhǎng)為785 nm，連續(xù)輸出功率120 mW，分辨率為8 cm-1，光譜采集范圍200～2 000 cm-1，陣列CCD2048個(gè)像素。

人工溫室環(huán)境下，水基培育玉米。玉米三葉一芯時(shí)采樣，距離葉頂端5 cm處截取葉樣(此處葉脈不明顯并且平整)，采集同一樣品近紅外光譜和拉曼光譜，如圖1和圖2所示。另一組樣品采用堿脅迫培養(yǎng)(氫氧化鈉)用來(lái)模擬大慶地區(qū)堿性生長(zhǎng)環(huán)境，目測(cè)堿脅迫光譜圖同圖1和圖2無(wú)差異。

Fig.1 Raman spectrum of corn leaf

Fig.2 Near infrared spectrum of corn leaf

2 光譜EMD的高頻能量吸附屬性

記拉曼光譜和近紅外光譜為l(ν)和i(ν)，在安裝有EMD Tool-box的Matlab2011b環(huán)境下經(jīng)驗(yàn)?zāi)B(tài)分解拉曼光譜和近紅外光譜，分別得到9和7個(gè)IMF分量，即

其中IMF為分解出來(lái)的調(diào)制窄帶信號(hào)。

定義重量和比重函數(shù)來(lái)衡量分解分量的重要程度： 任意光譜分量IMFj重量

分量比重

分別計(jì)算拉曼光譜和近紅外光譜的比重，如表1所示。

Table 1 Weight distribution of Raman spectrum and near infrared spectrum components

表中q(l)為拉曼光譜分量比重，比重分布范圍3.3%～20.5%，能量分布均勻。

q(i)為近紅外光譜分量比重，可以看出光譜中有效成分集中在IMF3～I(xiàn)MF7。

重點(diǎn)考察一階IMF分量，如圖3和圖4所示，圖中(a)為一階IMF分量，(b)為扣除一階IMF后剩余信號(hào)。

Fig.3 IMF1 and remaining components of Raman spectrum decomposed by EMD

Fig.4 IMF1 and remaining components of near infrared spectrum decomposed by EMD

通過(guò)圖1和圖3的剩余分量部分對(duì)比觀察，EMD對(duì)拉曼光譜幅值改變很大，即EMD基本上將拉曼光譜視作AM(amplitude modulated)調(diào)制信號(hào)，拉曼光譜中的大部分能量吸附到第一階IMF中； 從圖2和圖4看到，EMD對(duì)近紅外光譜幅值改變不大； 圖4中，IMF1是絕對(duì)高頻信號(hào)，剩余分量是相對(duì)低頻信號(hào)，EMD將近紅外光譜信號(hào)視作FM(frequency modulation)信號(hào)。這樣的結(jié)論是否準(zhǔn)確呢？

兩類(lèi)光譜的區(qū)別是近紅外為低頻信號(hào)，拉曼為高頻信號(hào)。信號(hào)采集過(guò)程，光譜中都會(huì)引入噪聲信號(hào)?？紤]使用兩個(gè)單頻率信號(hào)分別模擬理想光譜信號(hào)和噪聲信號(hào)，其表達(dá)形式如下：

S(n)=A1cos(2πf1n)，N(n)=A2cos(2πf2n)，X1(n)=S(n)+N(n)，X2(n)=S(n)×N(n)，為了更好地模擬光譜與噪聲，設(shè)定A1/A2=1 000；X1可以表示一部分幅度調(diào)制信號(hào)，用來(lái)模擬拉曼信號(hào)，其中光譜信號(hào)與噪聲信號(hào)頻率相對(duì)接近f1/f2=0.8； 用X2表達(dá)近紅外信號(hào)，是一個(gè)調(diào)頻信號(hào)，其中f1/f2=0.005。拉曼和近紅外模擬信號(hào)如圖5和圖6所示。

Fig.5 Ideal Raman spectrum

Fig.6 Ideal near infrared spectrum

EMD分解如圖7和圖8所示。

圖7中的幅值情況再一次說(shuō)明了拉曼模擬信號(hào)中大部分能量被吸附到一階IMF中，稱(chēng)為高頻能量吸附特性，這時(shí)EMD更多的視拉曼光譜為AM信號(hào)； 圖8中IMF1頻率比剩余分量高，IMF1基本保持了噪聲性狀，EMD視近紅外光譜為FM信號(hào)。

Fig.7 IMF1 and remaining components of ideal Raman spectrum decomposed by EMD

Fig.8 IMF1 and remaining components of ideal near infrared spectrum decomposed by EMD

通過(guò)以上實(shí)驗(yàn)得出，對(duì)拉曼光譜的低階IMF低通濾波、閾值濾波或者直接截掉低階IMF以實(shí)現(xiàn)去噪，會(huì)使拉曼光譜能量損失很大，光譜形變； 用上述方法處理近紅外光譜結(jié)果差異很大，可以實(shí)現(xiàn)簡(jiǎn)單快速去噪。

3 光譜模態(tài)混疊現(xiàn)象

EMD的實(shí)質(zhì)是按照信號(hào)的頻率，以特征尺度為依據(jù)解析信號(hào)為IMF分量，這意味著IMF為短時(shí)單頻率組分信號(hào)。拉曼和近紅外光譜EMD分解后，對(duì)一階IMF希爾伯特變換，瞬時(shí)頻率如圖9和圖10所示。

由圖9和圖10可見(jiàn)，拉曼光譜的一階IMF中頻率組分比近紅外光譜復(fù)雜，這是否意味著EMD分解拉曼光譜較近紅外光譜更易出現(xiàn)模態(tài)混疊現(xiàn)象，即IMF不能?chē)?yán)格滿(mǎn)足單頻率組分的條件，或不同時(shí)間尺度的信號(hào)被強(qiáng)制分解到同一IMF中，如圖11所示為拉曼光譜EMD分解情況。

圖11中IMF2和IMF3中均出現(xiàn)了差異很大的特征時(shí)間尺度，以及IMF4和IMF5中存在相近特征尺度，圖中存在模態(tài)混疊現(xiàn)象。經(jīng)過(guò)多組不同種類(lèi)樣品(豆油、水稻葉片、油脂、土壤等)拉曼光譜分解實(shí)驗(yàn)表明，拉曼光譜極易出現(xiàn)EMD模態(tài)混疊問(wèn)題[4]。

Fig.9 Instantaneous frequency of Raman spectrum IMF1

Fig.10 Instantaneous frequency of near infrared spectrum IMF1

Fig.11 Raman spectrum decomposed by EMD

拉曼光譜的間歇斷續(xù)使得應(yīng)該被分解到不同時(shí)間尺度的分量被強(qiáng)制分解到同一IMF中，即在信號(hào)斷續(xù)處產(chǎn)生了瞬時(shí)頻率的畸變，極大影響拉曼光譜瞬時(shí)頻率的正確表達(dá)，為拉曼光譜的進(jìn)一步時(shí)頻分析造成了障礙。有兩種方法可以解決模態(tài)混疊問(wèn)題： (1)改變EMD篩分規(guī)則以克服模態(tài)混疊[5-6]； (2)在待分解光譜中加入屏蔽信號(hào)，利用兩信號(hào)之間相互作用解決模態(tài)混疊問(wèn)題，總體均值經(jīng)驗(yàn)?zāi)７纸?ensemble empirical mode decomposition，EEMD)就是一個(gè)有效的方法[7]。

4 光譜EMD時(shí)頻分析

綜合考慮EMD分解時(shí)，拉曼光譜重量分布情況、高頻吸附特性和易出現(xiàn)模態(tài)混疊現(xiàn)象，以下僅針對(duì)近紅外光譜作出EMD時(shí)頻分析。

Fig.12 Reconstructed spectrum signal on proportion and original spectrum signal

目測(cè)兩信號(hào)擬合良好。只是在光譜開(kāi)始和結(jié)束位置，i′(ν)信號(hào)相對(duì)光滑，意味著IMF1和IMF2集中承載了原光譜中噪聲分量。

可以用以下兩個(gè)指標(biāo)度量光譜重構(gòu)擬合情況。

均方根誤差

相關(guān)系數(shù)

兩個(gè)指標(biāo)反映出來(lái)重構(gòu)精度很好，i′(ν)可以表達(dá)i(ν)。

下面通過(guò)IMF3～I(xiàn)MF7時(shí)頻分析近紅外光譜。首先分析IMF7，一般稱(chēng)之為趨勢(shì)項(xiàng)，如圖13所示，該項(xiàng)總體呈緩慢下降的趨勢(shì)，這表示隨著近紅外波長(zhǎng)的減小，葉片對(duì)光線的吸收呈現(xiàn)遞增趨勢(shì)，而從原始光譜很難得到這種趨勢(shì)規(guī)律。另外，計(jì)算IMF7的導(dǎo)數(shù)，從趨勢(shì)項(xiàng)變化情況看出，波數(shù)600～1 400(為繪圖方便，文中波數(shù)從0開(kāi)始，實(shí)際光譜的波數(shù)從4 000開(kāi)始)區(qū)間變化加劇，對(duì)應(yīng)圖12，光譜的主要透射峰剛好落在該區(qū)間。所以縮小分析范圍，在600～1 400區(qū)間對(duì)比辨識(shí)非堿脅迫與堿脅迫兩類(lèi)光譜。

Fig.13 Components of IMF3～I(xiàn)MF7

分量IMF3～I(xiàn)MF6經(jīng)過(guò)Hilbert變換，得到波長(zhǎng)頻率分譜，將波長(zhǎng)頻率譜求和，如圖14所示。

Fig.14 Wavenumber and frequency spectrum by HT

圖14(a)為堿脅迫光譜，圖14(b)為常規(guī)水基培養(yǎng)光譜。對(duì)比觀察，兩幅圖譜略有差異，相似度很高。在600～1 400 cm-1波數(shù)考察范圍內(nèi)，圖14(a)657 cm-1處多出一個(gè)特征頻率。實(shí)驗(yàn)中全部四個(gè)濃度堿基樣品共100幅光譜，均獲得了該特征頻率，說(shuō)明通過(guò)657 cm-1特征頻率可以辨識(shí)出堿脅迫近紅外光譜。

5 結(jié) 論

(1) EMD分解光譜，近紅外和拉曼兩種類(lèi)型光譜存在如下差異：

①一階IMF能量大，EMD對(duì)拉曼光譜幅值影響大，拉曼光譜可視作調(diào)幅信號(hào)，拉曼光譜具有高頻能量吸附特性；

②一階IMF頻率高，EMD對(duì)近紅外光譜幅值影響不大，近紅外光譜可視作調(diào)頻信號(hào)。

(2) EMD分解兩類(lèi)光譜時(shí)，拉曼光譜較近紅外光譜更容易出現(xiàn)模態(tài)混疊現(xiàn)象，無(wú)法實(shí)現(xiàn)窄帶單頻分解。

(3) 在近紅外光波段，近紅外光譜EMD趨勢(shì)項(xiàng)可以表達(dá)隨著波長(zhǎng)減小，葉綠素吸光度呈現(xiàn)遞增趨勢(shì)；

(4) 通過(guò)IMF分量的希爾伯特變換發(fā)現(xiàn)，堿脅迫葉片的近紅外光譜具有657 cm-1特征頻率，而非堿脅迫光譜不具備此特征。

進(jìn)一步研究方向：

(1) 用EMD、希爾伯特變換對(duì)近紅外光譜作出時(shí)頻分析，通過(guò)特征頻率辨識(shí)出堿脅迫和非堿脅迫兩類(lèi)光譜，實(shí)現(xiàn)了物質(zhì)的定性分析。進(jìn)一步通過(guò)研究特征頻率與物質(zhì)分子振動(dòng)和能級(jí)躍遷之間關(guān)系，以實(shí)現(xiàn)物質(zhì)的定量分析，期待用特征吸收峰和特征頻率雙重指標(biāo)提高近紅外檢測(cè)精度。

(2) 由于拉曼光譜高頻能量吸附特性和模態(tài)混疊問(wèn)題，EMD時(shí)頻分析方法不適合應(yīng)用到拉曼光譜解析，可以考慮應(yīng)用EEMD或其他EMD的改進(jìn)算法，相關(guān)問(wèn)題有待進(jìn)一步研究。

[1] Huang N E，Zheng Shen，Long S R，et al. Proc. R Soc. Lond A，1998，454: 903.

[2] CAI Jian-hua, WANG Xian-chun, HU Wei-wen(蔡劍華, 王先春, 胡惟文). Transactions of the Chinese Society for Agricultural Machinery(農(nóng)業(yè)機(jī)械學(xué)報(bào)), 2010, 41(9): 182.

[3] CAI Jian-hua, WANG Xian-chun(蔡劍華, 王先春). Acta Optica Sinica(光學(xué)學(xué)報(bào)), 2010, 30(1): 267.

[4] ZHAO Xiao-yu, FANG Yi-ming, WANG Zhi-gang(趙肖宇, 方一鳴, 王志剛). Spectroscopy and Spectral Analysis(光譜學(xué)與光譜分析), 2013, 33(12): 1.

[5] Huang N E, Wu M, Long S, et al. Proc. R Soc. Lond A, 2003, 459: 2317.

[6] Huang N E, Long S R. Annual Review of Fluid Mechanics, 1999, 31: 417.

[7] Wu Z H，Huang N E. Advances in Adaptive Data Analysis，2009, 1(1): 1.

EMD Time-Frequency Analysis of Raman Spectrum and NIR

ZHAO Xiao-yu1, FANG Yi-ming2, TAN Feng1, TONG Liang3, ZHAI Zhe4

1. College of Information Technology, Heilongjiang Bayi Agricultural University, Daqing 163319, China 2. College of Electrical Engineering, Yanshan University, Qinhuangdao 066004, China 3. Communication and Electronic Engineering Institute, Qiqihar University, Qiqihar 161006, China 4. Forestry Experiment Center of North China, Chinese Academy of Forestry, Beijing 102300, China

This paper analyzes the Raman spectrum and Near Infrared Spectrum (NIR) with time-frequency method. The empirical mode decomposition spectrum becomes intrinsic mode functions, which the proportion calculation reveals the Raman spectral energy is uniform distributed in each component, while the NIR’s low order intrinsic mode functions only undertakes fewer primary spectroscopic effective information. Both the real spectrum and numerical experiments show that the empirical mode decomposition (EMD) regard Raman spectrum as the amplitude-modulated signal, which possessed with high frequency adsorption property; and EMD regards NIR as the frequency-modulated signal, which could be preferably realized high frequency narrow-band demodulation during first-order intrinsic mode functions. The first-order intrinsic mode functions Hilbert transform reveals that during the period of empirical mode decomposes Raman spectrum, modal aliasing happened. Through further analysis of corn leaf’s NIR in time-frequency domain, after EMD, the first and second orders components of low energy are cut off, and reconstruct spectral signal by using the remaining intrinsic mode functions, the root-mean-square error is 1.001 1, and the correlation coefficient is 0.981 3, both of these two indexes indicated higher accuracy in re-construction; the decomposition trend term indicates the absorbency is ascending along with the decreasing to wave length in the near-infrared light wave band; and the Hilbert transform of characteristic modal component displays, 657 cm-1is the specific frequency by the corn leaf stress spectrum, which could be regarded as characteristic frequency for identification.

Time-frequency analysis; Empirical mode decomposition; Hilbert translation; Near infrared spectrum; Raman spectrum

Aug. 7, 2014; accepted Dec. 18, 2014)

2014-08-07，

2014-12-18

黑龍江省科學(xué)基金項(xiàng)目(QC2015071, F201329)，國(guó)家科技支撐項(xiàng)目(2014BAD06B01)，黑龍江省教育廳科學(xué)技術(shù)研究項(xiàng)目(12521378，12541897)資助

趙肖宇，女，1977年生，黑龍江八一農(nóng)墾大學(xué)副教授 e-mail: xy_zhao77@163.com

O657.3

A

10.3964/j.issn.1000-0593(2016)02-0424-06