大米拉曼光譜不同預處理方法的相近產地鑒別研究

王亞軒，譚 峰，辛元明，李 歡，趙肖宇，鹿保鑫

1.黑龍江八一農墾大學土木水利學院，黑龍江 大慶 163319 2.黑龍江八一農墾大學電氣與信息學院，黑龍江 大慶 163319 3.黑龍江八一農墾大學食品學院，黑龍江 大慶 163319

引 言

大米是我國主要的主食來源[1-2]，全國大米種植區域廣、種類多，土壤、環境和水質等差異形成地域因素會導致大米的品質發生變化。如五常大米、牡丹江響水大米，獨特的地理環境形成特有的口感和營養價值，使其成為具有鮮明地理標識的大米產品。但一些商家為了追求更高的利潤，用相近產地的大米代替地域品牌大米，購買時僅通過消費者肉眼判斷很難區分，這不僅損害了糧農的利益，也不利于品牌產業鏈的健康發展。因此，研究相近產地大米的快速準確無損鑒別的方法能為鑒別地理標識大米提供理論和技術支持。

拉曼光譜通過物質內部分子對可見單色光的散射強度不同來識別分子結構，從而對物質內部官能團進行特定指紋標定，光譜譜峰強度與分子濃度有關。目前已廣泛應用在食品、藥材、化工、寶石等多個領域進行定性或定量的檢測。拉曼光譜用于農產品檢測方面，主要集中在對糧食、奶制品、果蔬類、食用油等的研究上，通過拉曼光譜分析產品內部是否摻雜其他物質，進行農產品質量和年份的鑒定。拉曼光譜應用于大米檢測方面[3-6]，主流做法是通過光譜采集樣本的原始特征光譜，再去除熒光和噪聲，將樣本分為訓練集和測試集，結合主成分分析和偏最小二乘法進行數學模型的建立，來判別大米的產地、品種、新陳度等指標。黃嘉榮[3]等對廣東大米、東北大米及糯米進行分類，識別率是97.9%，孫娟[4]等對大米進行產地分類，選擇黑龍江大米、江蘇大米、湖南大米三地大米識別率為94%以上。趙迎[5]等對儲存三年以上和當年大米進行新陳大米進行分類，識別率為95%。從以上分析可以看出，研究主要是集中在不同品種大米的種類區分、對南方和北方產地大米的產地區分、不同年份大米的新陳度區分，而基于相近產地對大米進行分類鮮有研究。因為光譜鑒別不可避免要引入機器的噪聲和熒光背景等干擾因素，因地域相近大米內部的淀粉、糖類等物質含量差異不大，從光譜中提取這些結構特征性片段難度很大，需要通過有效的預處理算法去除干擾，提取真實準確的拉曼特征峰。本文研究比較四類九種不同的預處理方法結合偏最小二乘法建模，提出一種鑒別相近產地大米的預處理方法，為大米產地鑒別提供新的理論依據。

1 實驗部分

1.1 儀器與軟件

光譜采集使用廈門奧譜天成光電有限公司制造的波長785 nm便攜式拉曼光譜儀1臺，檢測范圍在124.79～3 324.66 cm-1，在最佳測量條件下，測量標準峰的位移值偏差為零，符合位移準確度不超過±4 cm-1的使用要求；大米脫殼采用上海超星LJJM精米機1臺，脫殼率≥99%，工作電壓220 V，試驗用量50～170 g；數據處理軟件為matlab2010b。

1.2 樣本

大米樣本于2018年11月采自黑龍江依安縣田間，分別是富饒鄉黎明村(北緯47.389 49、東經125.406 00)、新興鎮東萊村(北緯47.752 09、東經125.187 28)、上游鄉紅五月村(北緯47.933 883、東經125.322 755)相同品種的粳米，依次用A，B和C表示上述的三個產地大米，每個產地均隨機選取50個脫殼后的表面完好的大米作為試驗樣本，3個產地共計150個樣本。其中選擇每個產地樣本數的2/3即33個樣本用作訓練集，剩余的1/3即17個樣本用作測試集，共計51個樣本用于測試。

1.3 光譜的獲取

將從田間采集的帶殼稻米裝入尼龍網兜，在實驗室晾曬10 d后，采用統一加工對其用精米機進行兩次脫殼、每次脫殼50 s，再用100目篩子過篩，篩選出其中表面光滑完整的大米胚乳(去除胚芽)作為樣本。拉曼光譜檢測參數設置為：激光功率300 mW，激發波長785 nm，分辨率為6.58 cm-1，積分時間為5 000 ms，掃描范圍為200～3 300 cm-1的波段，測試條件為室溫，相對濕度為55%。每個樣本選擇米粒中間區域的背部或腹部采集數據，連續進行4次采集，取其平均值作為每個樣本的存儲數據。

1.4 光譜的預處理方法

光譜中普遍存在著熒光和背景噪聲，僅靠儀器的精度和準確度來消除檢測干擾受到儀器自身的限制，需要結合數學處理原始光譜數據來去除噪聲和基線漂移，常用的方法有導數處理、平移平滑、多項式擬合、歸一化等。導數處理主要是扣除儀器背景或漂移(散射)對信號的影響；平移平滑、多項式擬合能夠非常有效的提高譜圖信噪比，降低隨機噪聲的影響；歸一化可以消除尺度差異過大帶來的不良影響。

用大米樣本的原始光譜進行數據分析時，雖然可以用現有方法進行光譜數據預處理[7]，但其精度和準確度都達不到近地大米光譜鑒別的要求，試驗對比四類九種不同預處理方法進行原始數據分析，包括一階導數+平移平滑、二階導數+平移平滑、小波變換+去除基線三種常用的預處理方法，另外提出一種改進的分段多項式擬合+去除基線共四種預處理方法進行平滑去噪和去除基線漂移，再用極差歸一的方法進行單位統一，預處理后的數據分別采用偏最小二乘(pratial least squares,PLS)方法[8-9]進行建模分析，旨在探尋研究一種適合近地大米光譜的預處理方法。

2 結果與討論

2.1 三個產地的大米原始拉曼光譜

不同產地大米的營養成分基本一致，但各自的含量差異導致強度不同。圖1所示為200～3 300 cm-1范圍內三個產地的典型大米原始拉曼光譜，可見不同產地的大米峰值強度不同，但產生峰值位置基本相同。

圖1 三個產地大米原始光譜圖Fig.1 Raw Raman spectrum of three producing area of rice

2.2 大米典型拉曼峰值指認

大米光譜特征峰[6，10]對應著內部化學鍵振動方式及大米中營養成分的差異，如圖2所示，采用多項式擬合去除背景后的大米拉曼光譜的明顯峰值出現在480，866，942，1 088，1 130，1 263，1 344，1 385，1 458，1 822和2 911 cm-1處，峰值對應大米內部的主要營養物質，480 cm-1為淀粉的骨架振動；866和942 cm-1為支鏈淀粉的C—O—H和C—O—H變形振動；1 088 cm-1為直鏈淀粉的C—O—H鍵彎曲振動；1 130 cm-1為糖的C—O鍵伸縮振動和C—O—H鍵彎曲變形振動；1 263 cm-1為蛋白質的酰胺Ⅲ帶C—N鍵伸縮振動；1 344 cm-1為糖的C—C鍵伸縮振動和C—O—H鍵彎曲變形振動；1 385 cm-1為淀粉的C—C鍵伸縮振動；1 458 cm-1為糖的C—H鍵彎曲振動；1 822 cm-1為淀粉的O—C—O鍵伸縮振動；2 911 cm-1為淀粉的H—C—H鍵和H—N—H鍵伸縮振動；由此可見，主要特征峰出現在200～1 900和2 800～3 000 cm-1這兩個位置區間，根據主要特征峰值出現的波段，選擇200～3 100 cm-1的全波段進行建模分析。

圖2 大米拉曼光譜主要特征峰Fig.2 Main characteristic peaks of rice Raman spectrum

2.3 大米拉曼光譜預處理方法

當前常用的預處理方法包括一階導數、二階導數、平移平滑、小波變換、多項式擬合等，結合大米光譜特征拉曼峰值的特點，下面選擇四類九種預處理方法對光譜數據進行處理。

2.3.1 一階導數+平移平滑的預處理方法

一階導數的數學表達式為

(1)

其中xi為第i個樣品的光譜峰值的縱坐標，g為步長，試驗中是離散點求導，采用步長為1。

再對一階導數后的光譜用移動平均法進行平滑處理，數學表達式為

(2)

式(2)中，2n+1為窗口大小、試驗中n取2；i從第3點開始，對xi-2，xi-1，xi，xi+1，xi+2五點求平均，然后賦值給xi，之后移動窗口，使i點遍歷整個光譜到3 098點結束，即完成了移動平均法的平滑處理。

通過一階導數消除了原始光譜曲線的平移和漂移，但同時曲線噪聲被放大，原有多處波峰消失，并改變了拉曼光譜的形狀。從圖3中可知，采用常規的一階導數+平移平滑的預處理方法，需要再結合平移平滑對每個樣本數據進行校正，消除數據中的噪音，突出顯示光譜特征。

圖3 一階導數+平移平滑的預處理方法Fig.3 Pre-processing method of first derivative+translation smoothing

2.3.2 二階導數+平移平滑的預處理方法

二階導數的數學表達式為

(3)

其中xi為第i個樣品的光譜峰值的縱坐標，g為步長，采用步長為1，再對二階導數后的光譜用移動平均法進行平滑處理。

常規的二階導數+平移平滑的預處理方法，在一階導數基礎上進行二階導數并結合平滑濾波處理，如圖4所示，因為二階導數是對一階導數處理后曲線再求拉曼強度的變化率，導致結果曲線峰值變小，特征譜峰不明顯甚至消失。

圖4 二階導數+平移平滑的預處理方法Fig.4 Pre-processing method of second derivative+translation smoothing

2.3.3 小波變換+去除基線的預處理方法

小波變換改善了傅里葉變換不能進行局部分析的缺陷，將信號用母小波函數ψ(t)經過不同的平移和壓縮分解成一系列小波，因為小波變換可以精細的對時域和頻域的細節進行放大，使其具有很好的自適應性，但母小波函數不具有唯一性又使得分析時需要不斷嘗試，往往依靠經驗和不斷試驗才能達到去噪和去除基線的目的。母小波函數的數學公式為[11]

(4)

其中a為壓縮因子，b為平移因子。大米光譜屬于離散光譜經過多次對比分析選擇效果最佳的信號進行處理，選取小波高通濾波采用db9小波基函數對原始光譜棱角8級分解，濾掉低頻背景信號，選擇硬閾值去噪，如圖5所示，經小波變換處理后的光譜基線得到了校正，但基線仍有一定程度的漂移現象，主要產生在波段[1 800，3 100]這段背景噪聲較大的區間。

圖5 小波變換+去除基線的預處理方法Fig.5 Pre-processing method of wavelet transform+baseline removal

2.3.4 分段多項式擬合+去除基線的預處理方法

在相近產地大米鑒別中，因大米內部物質成分相似度極高，必要的預處理可以去除噪聲，增強特征峰的強度，上述的三種預處理方法對熒光背景進行去除后，存在不能保持原有波峰的形狀或基線漂移去除的不徹底的現象，為了改善以上缺點，提出一種分段多項式擬合+去除基線的預處理方法，這種預處理方法能保證擬合曲線恰到好處的通過原始波形下方，改進了傳統的多項式擬合方法，對光譜區間進行分段，校正后的波形與原始波形最大限度的保持相似性。

(1)窗口半寬為w，各測點i對應值為yi，在(w+1，n-w)區間取yi的平均值，記為式(5)

(5)

(6)

(3)將迭代后的yi值連接成線，找出曲線所有區間的最小值，記為

yi

(7)

(5)將每個區間的yi連接起來，形成分段多點擬合方法基線，如圖6(a)所示；

(6)在相同的拉曼位移上，用原始光譜曲線的數值對應減掉用分段多點擬合法的yi數值，形成去除基線后的光譜。

圖6 分段式多項式擬合+去除基線的預處理方法Fig.6 Pre-processing method of piecewise polynomial fitting+baseline removal

再討論擬合的階數對波形的影響，如圖6(b)中分別進行3點2次擬合、3點3次擬合、3點5次擬合，可見擬合的次方越大，會使擬合曲線震蕩的越劇烈，如圖6(b)中[200，600]區間，階數越高偏移越大，而在600 cm-1以后，幾乎沒有影響，分析原因可能是[200，600]區間分峰的大小和波形所致，如圖6(c)所示為采用3點2次擬合去除基線后的光譜，更好的保持了原有的特征峰面積和特定值，為實現光譜定量分析打下理論基礎。

2.4 基于偏最小二乘法的不同預處理方法分類結果分析

為了對比上述不同預處理方法的優劣，每份樣本中隨機選取33個作為訓練集樣本、其余17個作為測試集樣本。采用偏最小二乘法進行建模分析。并采用相關系數(r)、均方誤差(MSE)、均方根誤差(RMSE)來評價預處理的效果，其中r越大、MSE和RMSE越小說明樣本的預處理效果越好。

表1是對不同預處理方法所作的統計結果，從表中可見，在訓練集中一階導數+平移平滑的預處理方法相關系數值最大、均方誤差和均方根誤差最小，3點2次擬合+去除基線的預處理方法相關系數值稍差，但與一階導數+平移平滑差距不明顯，小波變換+去除基線的預處理方法相關系數值最小、均方誤差和均方根誤差最大；在測試集中采用3點2次擬合+去除基線的預處理方法的相關系數值最大、均方誤差和均方根誤差最小，3點3次擬合+去除基線的預處理方法稍差，二階導數+平移平滑的預處理方法最差。經過綜合比較，采用3點2次擬合+去除基線的預處理方法在訓練集和測試集中都是比較理想的預處理方法。

表1 不同預處理方法的相關系數CC、均方誤差MSE、均方根誤差RMSETabel 1 Correlation coefficient，Mean square error，Root mean square error of different pretreatment methods

為了進一步驗證不同預處理效果的差異，對3個產地樣品共150份大米采用PLS進行建模分析，在訓練集中，采用表1中的9種預處理方法對A，B和C三種大米的正確判別率均為100%。在測試集中如表2所示：采用3點2次擬合+去除基線預處理方法對A，B和C三產地大米總識別率為100%，采用5點2次擬合+去除基線預處理方法對A，B和C三產地大米總識別率為52.9%，其他分段多項式擬合介于二者之間；采用一階導數+平移平滑、二階導數+平移平滑和小波變換的預處理方法總識別率分別為88.2%，86.2%和96.1%；從中發現，采用一階導數+平移平滑的方法稍好于二階導數+去除基線的方法，這是因為二階導數的噪聲使更多特征峰不能突顯出來，導數處理不如小波變換和3點2次擬合+去除基線的效果，但小波變換過程需要通過先驗知識確定的參數過多，沒有通用規律可循，分段式多項式擬合中的3點2次擬合+去除基線的預處理方法優勢明顯，與表1中r，MSE和RMSE的結果吻合，總體識別率高，鑒別效果穩定。

表2 17個測試樣本中不同預處理方法的識別個數和識別率Table 2 The number and recognition rate of different pre-processing methods in 17 test sample

采用3點2次擬合+去除基線的預處理方法進行建模，并分別將A，B和C三產地大米樣本賦值1，2和3，結果在1±0.5(不含1.5)鑒別為A大米、結果在2±0.5(不含2.5)鑒別為B大米、結果在3±0.5(不含3.5)鑒別為C大米，結果如圖7所示。A大米的測試值主要集中在0.69～1.02、B大米樣本的測試值主要集中在1.54～2.01、C大米的測試值主要集中在2.75～3.01，均具有明顯的聚類趨勢。說明該模型預測結果具有較好的精度，可以很好的實現三種近地大米的產地鑒別。

圖7 真值與預測值關系圖Fig.7 Relationship between true value and predicted value

3 結 論

拉曼光譜技術結合不同預處理方法對相近三個產地的大米進行鑒別，分別采用一階導數+平移平滑、二階導數+平移平滑、小波變換+去除基線的方法進行光譜預處理，因為這些方法存在不能保持原有波峰的形狀或基線漂移的現象，提出一種分段多項式擬合+去除基線的預處理方法，通過偏最小二乘法PLS對150個樣本三個產地大米建立拉曼模型，實驗結果表明經過分段多項式擬合+去除基線中的3點2次多項式的預處理后建立的模型精度最高，在訓練集和測試集中三個產地的識別率均為100%，聚類效果好。通過3點2次多項式+去除基線的預處理為相近產地大米鑒別分析提供了一種有效方法，同時為近地域其他農作物鑒別提供技術參考。