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

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

引 言

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

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

1 實(shí)驗(yàn)部分

1.1 儀器與軟件

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

1.2 樣本

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

1.3 光譜的獲取

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

1.4 光譜的預(yù)處理方法

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

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

2 結(jié)果與討論

2.1 三個產(chǎn)地的大米原始拉曼光譜

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

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

2.2 大米典型拉曼峰值指認(rèn)

大米光譜特征峰[6，10]對應(yīng)著內(nèi)部化學(xué)鍵振動方式及大米中營養(yǎng)成分的差異，如圖2所示，采用多項(xiàng)式擬合去除背景后的大米拉曼光譜的明顯峰值出現(xiàn)在480，866，942，1 088，1 130，1 263，1 344，1 385，1 458，1 822和2 911 cm-1處，峰值對應(yīng)大米內(nèi)部的主要營養(yǎng)物質(zhì)，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為蛋白質(zhì)的酰胺Ⅲ帶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鍵伸縮振動；由此可見，主要特征峰出現(xiàn)在200～1 900和2 800～3 000 cm-1這兩個位置區(qū)間，根據(jù)主要特征峰值出現(xiàn)的波段，選擇200～3 100 cm-1的全波段進(jìn)行建模分析。

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

2.3 大米拉曼光譜預(yù)處理方法

當(dāng)前常用的預(yù)處理方法包括一階導(dǎo)數(shù)、二階導(dǎo)數(shù)、平移平滑、小波變換、多項(xiàng)式擬合等，結(jié)合大米光譜特征拉曼峰值的特點(diǎn)，下面選擇四類九種預(yù)處理方法對光譜數(shù)據(jù)進(jìn)行處理。

2.3.1 一階導(dǎo)數(shù)+平移平滑的預(yù)處理方法

一階導(dǎo)數(shù)的數(shù)學(xué)表達(dá)式為

(1)

其中xi為第i個樣品的光譜峰值的縱坐標(biāo)，g為步長，試驗(yàn)中是離散點(diǎn)求導(dǎo)，采用步長為1。

再對一階導(dǎo)數(shù)后的光譜用移動平均法進(jìn)行平滑處理，數(shù)學(xué)表達(dá)式為

(2)

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

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

圖3 一階導(dǎo)數(shù)+平移平滑的預(yù)處理方法Fig.3 Pre-processing method of first derivative+translation smoothing

2.3.2 二階導(dǎo)數(shù)+平移平滑的預(yù)處理方法

二階導(dǎo)數(shù)的數(shù)學(xué)表達(dá)式為

(3)

其中xi為第i個樣品的光譜峰值的縱坐標(biāo)，g為步長，采用步長為1，再對二階導(dǎo)數(shù)后的光譜用移動平均法進(jìn)行平滑處理。

常規(guī)的二階導(dǎo)數(shù)+平移平滑的預(yù)處理方法，在一階導(dǎo)數(shù)基礎(chǔ)上進(jìn)行二階導(dǎo)數(shù)并結(jié)合平滑濾波處理，如圖4所示，因?yàn)槎A導(dǎo)數(shù)是對一階導(dǎo)數(shù)處理后曲線再求拉曼強(qiáng)度的變化率，導(dǎo)致結(jié)果曲線峰值變小，特征譜峰不明顯甚至消失。

圖4 二階導(dǎo)數(shù)+平移平滑的預(yù)處理方法Fig.4 Pre-processing method of second derivative+translation smoothing

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

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

(4)

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

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

2.3.4 分段多項(xiàng)式擬合+去除基線的預(yù)處理方法

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

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

(5)

(6)

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

yi

(7)

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

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

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

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

2.4 基于偏最小二乘法的不同預(yù)處理方法分類結(jié)果分析

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

表1是對不同預(yù)處理方法所作的統(tǒng)計(jì)結(jié)果，從表中可見，在訓(xùn)練集中一階導(dǎo)數(shù)+平移平滑的預(yù)處理方法相關(guān)系數(shù)值最大、均方誤差和均方根誤差最小，3點(diǎn)2次擬合+去除基線的預(yù)處理方法相關(guān)系數(shù)值稍差，但與一階導(dǎo)數(shù)+平移平滑差距不明顯，小波變換+去除基線的預(yù)處理方法相關(guān)系數(shù)值最小、均方誤差和均方根誤差最大；在測試集中采用3點(diǎn)2次擬合+去除基線的預(yù)處理方法的相關(guān)系數(shù)值最大、均方誤差和均方根誤差最小，3點(diǎn)3次擬合+去除基線的預(yù)處理方法稍差，二階導(dǎo)數(shù)+平移平滑的預(yù)處理方法最差。經(jīng)過綜合比較，采用3點(diǎn)2次擬合+去除基線的預(yù)處理方法在訓(xùn)練集和測試集中都是比較理想的預(yù)處理方法。

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

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

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

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

圖7 真值與預(yù)測值關(guān)系圖Fig.7 Relationship between true value and predicted value

3 結(jié) 論

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