李路，黃漢英，趙思明，楊素仙

1(華中農(nóng)業(yè)大學(xué) 工學(xué)院，湖北 武漢，430070)2(農(nóng)業(yè)部長江中下游農(nóng)業(yè)裝備重點(diǎn)實(shí)驗(yàn)室，湖北 武漢，430070) 3(華中農(nóng)業(yè)大學(xué) 食品科技學(xué)院，湖北 武漢，430070)

稻谷的谷殼率和整精米率是2個(gè)重要的品質(zhì)指標(biāo)[1-2]，與農(nóng)民和相關(guān)企業(yè)的經(jīng)濟(jì)效益密切相關(guān)。目前，稻谷谷殼率和整精米率的檢測方法主要有人工抽樣法[3-4]和圖像識(shí)別法[5-8]。這兩類方法都需要對(duì)稻谷樣本進(jìn)行礱谷、碾米等加工，樣本的準(zhǔn)備過程繁瑣，且不適用于快速無損檢測。因此，在稻谷的收獲、收獲、存儲(chǔ)等環(huán)節(jié)，急需一種操作簡單的稻谷加工品質(zhì)無損檢測方法。

相關(guān)文獻(xiàn)[9-12]表明，稻谷的谷殼率、整精米率等品質(zhì)很大程度上受到品種、生長環(huán)境、水分、直鏈淀粉含量等因素的影響，可見稻谷顆粒自身的特性，尤其是化學(xué)成分與上述指標(biāo)間存在一定的相關(guān)性。而近年來，國內(nèi)外有關(guān)稻谷化學(xué)成分的近紅外檢測研究較多[13-16]，相關(guān)技術(shù)也較成熟。

因此，本文針對(duì)稻谷谷殼率和整精米率的預(yù)測問題，采用近紅外漫反射光譜信息，根據(jù)篩選出的相關(guān)特征波長，建立稻谷谷殼率和整精米率的近紅外預(yù)測模型，為實(shí)現(xiàn)稻谷加工品質(zhì)的無損檢測提供一種有效的手段。

1 材料與方法

1.1 試驗(yàn)材料與設(shè)備

試驗(yàn)材料包含：中9A/R591、A4A/R326、廣占S/R166、岡紅1A/R15等46個(gè)稻谷樣本，產(chǎn)自海南省，收獲至試驗(yàn)在60 d內(nèi)完成。

試驗(yàn)設(shè)備為漫反射式Supnir-2720近紅外光譜儀,聚光科技杭州股份有限公司;JLGJ4.5型檢驗(yàn)礱谷機(jī),臺(tái)州市糧儀廠;JNMJ3型檢驗(yàn)?zāi)朊讬C(jī),臺(tái)州市糧儀廠。

1.2 近紅外光譜采集

近紅外光譜儀參數(shù)如下：儀器帶寬1 nm，掃描間隔1 nm，波長范圍1 000～1 799 nm，光譜數(shù)據(jù)點(diǎn)數(shù)800，掃描溫度15～25 ℃。

近紅外光譜儀預(yù)熱30 min，經(jīng)性能測試和參比后進(jìn)行光譜掃描。將整粒稻谷放入樣本盤中，裝滿壓實(shí)，每個(gè)樣本掃描3次，最終的光譜數(shù)據(jù)取平均值。

1.3 谷殼率與整精米率測定

將試驗(yàn)材料經(jīng)礱谷、去殼和碾米等工序制成精米，收集加工過程中的谷殼與整精米。

1.3.1 谷殼率測定

設(shè)待加工的稻谷質(zhì)量記為m1，礱谷加工時(shí)得到的谷殼質(zhì)量記為m2，單位為g。稻谷谷殼率具體計(jì)算公式為：

(1)

1.3.2 整精米率測定

設(shè)稻谷樣本加工后的整精米總質(zhì)量為m3，單位為g，按照GB/T 21719—2008[4]中的計(jì)算公式計(jì)算稻谷整精米率：

(2)

1.4 近紅外光譜消噪

采集到的近紅外光譜，除含有樣本的信息外，還包含了一些噪聲，運(yùn)用原始近紅外光譜數(shù)據(jù)進(jìn)行建模，會(huì)使模型的準(zhǔn)確性降低。因此，對(duì)采集到的近紅外原始光譜進(jìn)行預(yù)處理十分必要。

小波變換可以將光譜信息分解為背景信息(低頻)、組分信息(中頻)和噪音(高頻)部分，據(jù)此可進(jìn)行去噪[17]。為確定小波分解尺度，同時(shí)考慮到預(yù)處理過程的運(yùn)算量不宜過大，本文運(yùn)用2階Daubechies小波(db2)為母小波，分解尺度為2和3，將樣品1 000～1 799 nm的近紅外光譜作為原始信號(hào)進(jìn)行小波消噪處理，根據(jù)定標(biāo)標(biāo)準(zhǔn)差的最小值來挑選一個(gè)最佳分解尺度。

1.5 光譜預(yù)處理

小波消噪后，為了有效放大光譜信息，將光譜壓縮在相同的范圍進(jìn)行比較，消除基線漂移、樣品不均勻、光散射、光程變化等對(duì)光譜的影響，提高模型的預(yù)測能力和精度，需要對(duì)光譜進(jìn)行預(yù)處理。本文對(duì)Z-score歸一化[18]一階導(dǎo)、二階導(dǎo)和一階導(dǎo)+Z-score歸一化等4種預(yù)處理方法進(jìn)行試驗(yàn)對(duì)比，然后根據(jù)不同預(yù)處理方法所建立模型的決定系數(shù)和定標(biāo)標(biāo)準(zhǔn)差來確定最終的光譜處理方案。

1.6 樣本集劃分

為建立模型和對(duì)模型性能進(jìn)行驗(yàn)證，需要將樣本集劃分為訓(xùn)練集和驗(yàn)證集，理想的訓(xùn)練集的數(shù)據(jù)范圍應(yīng)該包含驗(yàn)證集的數(shù)據(jù)范圍，且2集的平均值和標(biāo)準(zhǔn)差應(yīng)該相差不大，這樣建立的模型預(yù)測的準(zhǔn)確率更高。

本文用Kennard-Stone法[19]對(duì)樣本集按照4∶1的比例進(jìn)行劃分，該方法首先將所有的樣本都看作訓(xùn)練集，然后從所有的樣本中依次挑選出光譜差距大的樣本進(jìn)入訓(xùn)練集，直到訓(xùn)練集的樣本數(shù)量達(dá)到要求，剩余的樣本則自動(dòng)歸為驗(yàn)證集。具體計(jì)算方法為：(1)計(jì)算兩兩稻谷樣本的近紅外光譜數(shù)據(jù)之間的歐式距離，將歐式距離最大的2個(gè)稻谷樣本選入訓(xùn)練集；(2)根據(jù)剩余稻谷樣本與已選入訓(xùn)練集的稻谷樣本的歐式距離，選取距離最短和最長的樣本加入訓(xùn)練集；(3)重復(fù)(2)的操作，直到訓(xùn)練集的樣本數(shù)量達(dá)到要求為止。稻谷樣本近紅外光譜數(shù)據(jù)的歐式距離計(jì)算公式為：

(3)

式中：xi，稻谷樣本近紅外光譜吸光值；a、b，分代表2個(gè)樣本。

1.7 特征波長篩選

稻谷樣本的近紅外光譜由800個(gè)數(shù)據(jù)點(diǎn)組成，而樣本數(shù)為46個(gè)，自變量的個(gè)數(shù)遠(yuǎn)多于因變量，在建模的過程中，共線性非常嚴(yán)重。且采集到的近紅外光譜含有很多冗余信息，若采用全光譜建立模型，計(jì)算工作量大。因此，若能挑選出稻谷近紅外光譜的特征波長，運(yùn)用特征波長進(jìn)行建模可有效提高模型的預(yù)測精度和穩(wěn)定性。本文運(yùn)用競爭自適應(yīng)重加權(quán)采樣(competitive adaptive reweighted sampling, CARS)[20]方法篩選特征波長。

CARS法在特征波長篩選過程中，首先利用蒙特卡羅采樣法采樣N次，建立偏最小二乘回歸(partial least squares regression，PLSR)模型；然后根據(jù)使用指數(shù)衰減函數(shù)強(qiáng)行去掉回歸系數(shù)相對(duì)較小的波長點(diǎn)；再用每次采樣保留的波長建立PLSR模型；最后選擇模型的交叉驗(yàn)證均方差(root mean square error of cross validation，RMSECV)最小時(shí)所對(duì)應(yīng)的波長子集作為特征波長。具體算法為：

T=XW

(4)

y=Tc+e=XWc+e=Xb+e

(5)

式中：X，46×800的光譜數(shù)據(jù)矩陣；y，46×1的稻米加工品質(zhì)指標(biāo)向量；T，X的得分矩陣，是X與組合系數(shù)W的線性組合；c，y和T建立的PLS校正模型的回歸系數(shù)向量；e，預(yù)測殘差；b， 1×800的系數(shù)向量。有下列關(guān)系式成立：

b=Wc=[b1,b2, …,bj]T

(6)

b中第j個(gè)元素的絕對(duì)值|bj| (1≤j≤800)表示第j個(gè)波長對(duì)y的貢獻(xiàn)，|bj|越大則該變量越重要，越應(yīng)該保留。

1.8 建模

運(yùn)用多元線性回歸(multiple linear regression，MLR)[21]方法建立稻谷谷殼率和整精米率的近紅外光譜預(yù)測模型。MLR方法建立的預(yù)測模型，公式含義清楚，對(duì)于自變量較少的模型回歸效果較好，并且能在建模時(shí)根據(jù)顯著性指標(biāo)對(duì)特征波長再次優(yōu)選，進(jìn)一步減少特征波長數(shù)量，降低模型復(fù)雜度。建模過程中，以決定系數(shù)R2、定標(biāo)標(biāo)準(zhǔn)差(root mean square error of calibration， RMSEC)、校驗(yàn)標(biāo)準(zhǔn)差(root mean square error of prediction，RMSEP)和相對(duì)偏差來評(píng)價(jià)模型的穩(wěn)定性和預(yù)測精度。R2越大，RMSEC和相對(duì)偏差越小，則模型的穩(wěn)定性和預(yù)測精度越好。

2 結(jié)果與分析

根據(jù)采用不同小波分解尺度進(jìn)行消噪試驗(yàn)的結(jié)果，當(dāng)分解尺度為2時(shí)RMSEC為0.003 03，當(dāng)分解尺度為3時(shí)RMSEC為0.003 01，可見分解尺度為3時(shí)性能稍好，故確定小波消噪時(shí)的分解尺度為3。圖1為46個(gè)稻谷樣本的近紅外吸收光譜經(jīng)消噪后的光譜圖。可見，不同品種的稻谷樣本，其近紅外光譜的變化趨勢大體是一致的。但由于不同樣本之間的化學(xué)成分具有微小差異，所以其光譜的吸光度略有不同。

圖1 稻谷樣本近紅外光譜圖Fig.1 Near-infrared spectra of paddy samples

表1為4種預(yù)處理方法的效果比較，運(yùn)用Z-score歸一化處理后所建模型的決定系數(shù)最大、定標(biāo)標(biāo)準(zhǔn)差最小，因此選用Z-score歸一化作為光譜數(shù)據(jù)預(yù)處理方法，其結(jié)果如圖2所示。

表1 不同預(yù)處理方法的比較Table 1 Comparison of different pretreatments

圖2 經(jīng)過Z-score歸一化后的光譜圖Fig.2 Spectra after Z-score normalize

表2為Kennard-Stone法選取訓(xùn)練集與驗(yàn)證集的結(jié)果見?？梢?，訓(xùn)練集和驗(yàn)證集的平均值和標(biāo)準(zhǔn)差相差不大，驗(yàn)證集的數(shù)據(jù)范圍包含在訓(xùn)練集的數(shù)據(jù)范圍內(nèi)，說明樣本集劃分均勻、合理。

表2 Kennard-Stone 分組結(jié)果Table 2 Results of Kennard-Stone

圖3-a～圖3-b為CARS方法對(duì)稻谷谷殼率特征波長的篩選過程，圖3-d～圖3-f為整精米率特征波長的篩選過程。由圖3-a和圖3-d可知，隨著運(yùn)行次數(shù)的增加，保留的特征波長數(shù)量呈負(fù)指數(shù)規(guī)律減少。圖3-b和圖3-e為采用10折交叉驗(yàn)證得到的RMSECV的變化趨勢，RMSECV值減小，說明剔除了無關(guān)變量，RMSECV值增大，說明剔除了有效變量。在圖3-b中，1～59次蒙特卡羅采樣過程中RMSECV呈現(xiàn)遞減趨勢， 60次后開始快速增大，因此采樣次數(shù)為100次就能選出最優(yōu)的特征波長數(shù)，此時(shí)保留的特征波長數(shù)為24個(gè)，說明經(jīng)過篩選，稻谷谷殼率的特征波長由800個(gè)減少到了24個(gè)。圖3-e中，1～57次采樣過程中RMSECV呈現(xiàn)遞減趨勢， 58～130次之間緩慢上升，為保證能篩選出最優(yōu)的特征波長，采樣次數(shù)設(shè)定為200次，此時(shí)稻谷整精米率的特征波長由800個(gè)減少到了31個(gè)。圖3-c和圖3-f中各曲線表示各特征波長的偏回歸系數(shù)隨蒙特卡羅采樣次數(shù)的變化趨勢，*代表了RMSECV最小時(shí)的特征波長數(shù)。

根據(jù)篩選出的特征波長，使用MLR方法建立稻谷谷殼率和整精米率的預(yù)測模型，并根據(jù)顯著性指標(biāo)，分別去掉了3個(gè)對(duì)谷殼率和整精米率不顯著的特征波長。最終，稻谷谷殼率的特征波長為21個(gè)，整精米率為28個(gè)。

表3為稻谷谷殼率MLR預(yù)測模型的參數(shù)，其回歸常數(shù)項(xiàng)b=354.11，xi為各特征波長所對(duì)應(yīng)的經(jīng)過預(yù)處理后的吸光值，ai為各特征波長的偏回歸系數(shù)。模型的R2為0.998 3，RMSEC為0.112 9，相對(duì)偏差為0.51%，說明模型具有較好的穩(wěn)定性和預(yù)測準(zhǔn)確度。由表3可知，在1 127、1 203、1 264、1 446、1 495、1 597 nm等6個(gè)特征波長處，偏回歸系數(shù)的絕對(duì)值最大，t值相對(duì)較大，p值相對(duì)較小，說明這些特征波長對(duì)谷殼率預(yù)測模型的影響較顯著。

圖3 稻谷谷殼率和整精米率的特征波長篩選圖Fig.3 Key wavelengths selection of husk content and head rice yield paddy

序號(hào)xim偏回歸系數(shù)aitp序號(hào)xi偏回歸系數(shù)aitp11 004-80.7-7.540121 39359.631.70.002 821 01395.776.940131 396-65.06-2.130.003 231 040-86.68-4.350.000 8141 418164.9619.62041 061135.269.910151 446-242.58-15051 086-53.16-4.470.000 6161 495235.1510.71061 127237.1618.70171 597-251.58-5.390.000 171 203221.3419.30181 598205.864.330.000 881 264-335.18-9.940191 697114.935.060.000 291 291174.013.850.002201 753-172.13-11.670101 293215.185.60211 792112.479.690111 311-110.11-5.970

表4為稻谷整精米率預(yù)測模型的參數(shù)，其回歸常數(shù)項(xiàng)b=-10 065。模型的R2為0.998 7，RMSEC為0.982 1，相對(duì)偏差為2.34%，說明模型的穩(wěn)定性和預(yù)測能力較好。由表4可知，在1 114、1 257、1 659、1 680 nm等4個(gè)特征波長處，偏回歸系數(shù)的絕對(duì)值最大，t值相對(duì)較大，p值相對(duì)較小，說明這些特征波長對(duì)整精米率預(yù)測模型的影響較顯著。

使用驗(yàn)證集稻谷樣本，對(duì)稻谷谷殼率和整精米率預(yù)測模型進(jìn)行驗(yàn)證，其結(jié)果如表5所示。谷殼率預(yù)測誤差的絕對(duì)值小于0.5%，整精米率預(yù)測誤差的絕對(duì)值最高達(dá)到了5.62%，可見前者的預(yù)測精度要高于后者。究其原因，主要是本試驗(yàn)所采集的近紅外漫反射光譜主要反映了稻谷顆粒淺表的性狀，谷殼位于稻谷顆粒的最外層，光譜特征與其相關(guān)性較大，故模型的預(yù)測精度高；而精米主要是指稻谷顆粒最里層的胚乳，近紅外漫反射光譜特征與其相關(guān)性稍小，故整精米率預(yù)測模型的精度相對(duì)較低。

稻谷谷殼率和整精米率預(yù)測模型驗(yàn)證時(shí)的R2分別為0.924 5和0.928 7，RMSEP分別為0.221 6和3.115 2，相對(duì)偏差分別為1.02%和7.90%。說明通過近紅外漫反射光譜能對(duì)稻谷的谷殼率和整精米率進(jìn)行有效的預(yù)測。

表4 稻谷整精米率預(yù)測模型的參數(shù)Table 4 Parameters of head rice yield prediction model

表5 驗(yàn)證集預(yù)測結(jié)果 單位：%

3 結(jié)論

以46個(gè)品種的稻谷樣本為研究對(duì)象，首先采集樣本的近紅外漫反射光譜，使用分解尺度為3的2階小波消噪和Z-score歸一化對(duì)光譜數(shù)據(jù)進(jìn)行預(yù)處理。然后利用Kennard-Stone法劃分了訓(xùn)練集和驗(yàn)證集，使用CARS法確定了與稻谷谷殼率和整精米率相關(guān)的特征波長。最后根據(jù)MLR理論建立了稻谷谷殼率和整精米率的近紅外光譜預(yù)測模型，并使用驗(yàn)證集樣本對(duì)模型進(jìn)行了驗(yàn)證。結(jié)果表明：

(1)稻谷谷殼率的近紅外特征波長為21個(gè)，其中最典型的特征波長為：1 127、1 203、1 264、1 446、1 495、1 597 nm；整精米率的特征波長為28個(gè)，最典型的有：1 114、1 257、1 659、1 680 nm。

(2)所建立的谷殼率和整精米率預(yù)測模型的R2分別為0.998 3和0.998 7，RMSEC分別為0.112 9和0.982 1，相對(duì)偏差分別為0.51%和2.34%。

(3)兩模型驗(yàn)證的R2分別為0.924 5和0.928 7，RMSEP分別為0.221 6和3.115 2，相對(duì)偏差分別為1.02%和7.90%。

綜上所述，利用近紅外漫反射光譜信息對(duì)稻谷谷殼率和整精米率進(jìn)行預(yù)測是可行的，該研究結(jié)果能為稻谷加工品質(zhì)的無損檢測提供一種有效的手段。