基于近红外漫反射光谱的稻谷谷壳率和整精米率预测

李路，黄汉英，赵思明，杨素仙

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%。

综上所述，利用近红外漫反射光谱信息对稻谷谷壳率和整精米率进行预测是可行的，该研究结果能为稻谷加工品质的无损检测提供一种有效的手段。