EMD時頻分析拉曼光譜和近紅外光譜

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

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

EMD時頻分析拉曼光譜和近紅外光譜

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

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

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

時頻分析； 經(jīng)驗模態(tài)分解； 希爾伯特變換； 近紅外光譜； 拉曼光譜

引 言

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

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

1 光譜采集

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

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

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)驗模態(tài)分解拉曼光譜和近紅外光譜，分別得到9和7個IMF分量，即

其中IMF為分解出來的調制窄帶信號。

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

分量比重

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

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

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

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

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

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

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

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

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

Fig.5 Ideal Raman spectrum

Fig.6 Ideal near infrared spectrum

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

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

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

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

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

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

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

圖11中IMF2和IMF3中均出現(xiàn)了差異很大的特征時間尺度，以及IMF4和IMF5中存在相近特征尺度，圖中存在模態(tài)混疊現(xiàn)象。經(jīng)過多組不同種類樣品(豆油、水稻葉片、油脂、土壤等)拉曼光譜分解實驗表明，拉曼光譜極易出現(xiàn)EMD模態(tài)混疊問題[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ù)使得應該被分解到不同時間尺度的分量被強制分解到同一IMF中，即在信號斷續(xù)處產生了瞬時頻率的畸變，極大影響拉曼光譜瞬時頻率的正確表達，為拉曼光譜的進一步時頻分析造成了障礙。有兩種方法可以解決模態(tài)混疊問題： (1)改變EMD篩分規(guī)則以克服模態(tài)混疊[5-6]； (2)在待分解光譜中加入屏蔽信號，利用兩信號之間相互作用解決模態(tài)混疊問題，總體均值經(jīng)驗模分解(ensemble empirical mode decomposition，EEMD)就是一個有效的方法[7]。

4 光譜EMD時頻分析

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

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

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

可以用以下兩個指標度量光譜重構擬合情況。

均方根誤差

相關系數(shù)

兩個指標反映出來重構精度很好，i′(ν)可以表達i(ν)。

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

Fig.13 Components of IMF3～IMF7

分量IMF3～IMF6經(jīng)過Hilbert變換，得到波長頻率分譜，將波長頻率譜求和，如圖14所示。

Fig.14 Wavenumber and frequency spectrum by HT

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

5 結 論

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

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

②一階IMF頻率高，EMD對近紅外光譜幅值影響不大，近紅外光譜可視作調頻信號。

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

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

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

進一步研究方向：

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

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

[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(農業(yè)機械學報), 2010, 41(9): 182.

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

[4] ZHAO Xiao-yu, FANG Yi-ming, WANG Zhi-gang(趙肖宇, 方一鳴, 王志剛). Spectroscopy and Spectral Analysis(光譜學與光譜分析), 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

黑龍江省科學基金項目(QC2015071, F201329)，國家科技支撐項目(2014BAD06B01)，黑龍江省教育廳科學技術研究項目(12521378，12541897)資助

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

O657.3

A

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