基于近紅外漫反射光譜的稻谷谷殼率和整精米率預測

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

1(華中農業大學 工學院，湖北 武漢，430070)2(農業部長江中下游農業裝備重點實驗室，湖北 武漢，430070) 3(華中農業大學 食品科技學院，湖北 武漢，430070)

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

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

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

1 材料與方法

1.1 試驗材料與設備

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

試驗設備為漫反射式Supnir-2720近紅外光譜儀,聚光科技杭州股份有限公司;JLGJ4.5型檢驗礱谷機,臺州市糧儀廠;JNMJ3型檢驗碾米機,臺州市糧儀廠。

1.2 近紅外光譜采集

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

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

1.3 谷殼率與整精米率測定

將試驗材料經礱谷、去殼和碾米等工序制成精米，收集加工過程中的谷殼與整精米。

1.3.1 谷殼率測定

設待加工的稻谷質量記為m1，礱谷加工時得到的谷殼質量記為m2，單位為g。稻谷谷殼率具體計算公式為：

(1)

1.3.2 整精米率測定

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

(2)

1.4 近紅外光譜消噪

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

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

1.5 光譜預處理

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

1.6 樣本集劃分

為建立模型和對模型性能進行驗證，需要將樣本集劃分為訓練集和驗證集，理想的訓練集的數據范圍應該包含驗證集的數據范圍，且2集的平均值和標準差應該相差不大，這樣建立的模型預測的準確率更高。

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

(3)

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

1.7 特征波長篩選

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

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

T=XW

(4)

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

(5)

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

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

(6)

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

1.8 建模

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

2 結果與分析

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

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

表1為4種預處理方法的效果比較，運用Z-score歸一化處理后所建模型的決定系數最大、定標標準差最小，因此選用Z-score歸一化作為光譜數據預處理方法，其結果如圖2所示。

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

圖2 經過Z-score歸一化后的光譜圖Fig.2 Spectra after Z-score normalize

表2為Kennard-Stone法選取訓練集與驗證集的結果見。可見，訓練集和驗證集的平均值和標準差相差不大，驗證集的數據范圍包含在訓練集的數據范圍內，說明樣本集劃分均勻、合理。

表2 Kennard-Stone 分組結果Table 2 Results of Kennard-Stone

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

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

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

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

序號xim偏回歸系數aitp序號xi偏回歸系數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為稻谷整精米率預測模型的參數，其回歸常數項b=-10 065。模型的R2為0.998 7，RMSEC為0.982 1，相對偏差為2.34%，說明模型的穩定性和預測能力較好。由表4可知，在1 114、1 257、1 659、1 680 nm等4個特征波長處，偏回歸系數的絕對值最大，t值相對較大，p值相對較小，說明這些特征波長對整精米率預測模型的影響較顯著。

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

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

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

表5 驗證集預測結果 單位：%

3 結論

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

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

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

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

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