Prediction and study of strawberry hardness based on hyperspectral index segmentation

草莓为多年生草本植物，原产南美，果实色泽鲜艳，汁多味甜，具有生津利咽、美容养颜、降血脂等保健功效^[1-2]，被誉为“水果皇后”。草莓皮薄，保质期较短，难以长时间运输和存储，因此常采摘硬度较高的果实, 硬度是表征草莓成熟度和货架时间的重要指标，也是水果品质的重要体现。为了实现草莓品质和新鲜程度的判断，人们做了很多工作。

为了研究草莓的质量特性和储存过程中的变化，人们探讨了不同品种、不同贮存条件下草莓的硬度变化规律和影响因素，为草莓的生产和质量控制提供了科学依据^[3-5]。这些研究大多通过硬度仪或质构仪对水果硬度进行测试^[3-7]，操作繁琐且对样本具有破坏性，而更快速的无损检测研究较少。近年来兴起的高光谱成像技术能够同时获取目标图像和光谱信息，具有光谱分辨率高、信息量大的特点，已广泛应用于水果、农产品等的硬度和品质检测^[8-11]。

由于高光谱图像的光谱通道数多，数据量大，因此需要对原始光谱降维以消除不相关或冗余的数据^[12]。XIONG等人通过主成分分析(principal component analysis, PCA)确定敏感波长来对高光谱数据降维，再通过最小二乘支持向量机建模，用于荔枝果实的微损伤检测^[13]。WU等人利用PCA和自适应加权采样算法结合偏最小二乘法实现桑葚果实硫菌素-甲基残留物的检测^[14]。ZHAO等人通过连续投影算法(successive projections algorithm, SPA)对原光谱降维以快速识别野生黑枸杞^[15]。这些研究都表明，降维算法能有效减少高光谱信息的冗余度，但不同降维方法获得的特征波长和信息不同，因此，在提取特征时，应选择能使草莓硬度预测模型更稳定和更准确的方法。

本文作者将不同降维方法和回归算法结合起来，共构建8种硬度回归模型，比较评价指数，以确定最优硬度预测模型。通过分析背景和草莓组分(萼片、草莓籽、果肉和霉变部位)的光谱特征, 选择3个特征光谱波段，构建基于特征波段反射率的归一化指数，实现果肉组分和其它组分的分割，结合最优模型实现草莓果肉硬度的重建。

实验中采用的多光谱成像系统，主要由一个高感度相机、两个150 W卤素光源、一个步进马达移动平台、机箱和计算机等部件组成，可采集280 nm~1167 nm共616个波段范围内的光谱信息，分辨率为1.3 nm。GY-3水果硬度计主要由硬度计表和一个GY-3测头组成。

实验样本为2022-05-03采自合肥某草莓园的新鲜草莓，采摘后3 h送至实验室，18:00进行数据采集实验，间隔24 h进行相同实验，共完成5天的数据采集。在剔除掉霉变及损坏的草莓后，最终保留457颗草莓的光谱图像和硬度数据。草莓分为草莓籽、霉变果肉、果肉、萼片不同组分，如图 1所示。

Figure 1.  Schematic diagram of strawberry samples and different compositions

将样本平整放在黑色吸光绒布上，置于载物台上，与相机垂直距离为70 cm。设置合适的推扫速度和曝光时间，并进行黑白校正，以减弱电流及系统不稳定的影响。

由于高光谱信息在400 nm前信息较紊乱，因此选用波段范围为400 nm~1167 nm的高光谱图像数据。选取草莓赤道区域50 pixel×50 pixel范围共2500个像素点作为感兴趣区域(region of interest, ROI)，以ROI范围平均光谱作为该草莓图像的光谱信息。为了对高光谱图像中草莓不同组分进行分割，需提取其高光谱数据，每个类型随机提取126个5 pixel×5 pixel范围ROI的平均光谱信息。

根据中国现行标准^[16]，以硬度计测头对草莓垂直施压, 当探头进入草莓内部10 cm刻度线时停止，读取并记录外圈数值，将赤道硬度值作为整颗草莓硬度, 如图 2a所示。每次测量结束后将指针归零，再进行下次测量。经实测，草莓样本硬度分布于0.89×10⁵ Pa~3.39×10⁵ Pa范围的区域内, 如图 2b所示。

Figure 2.  Sample hardness measurement

图 3为本文中方法的框图。主要由数据获取、建模、基于指数阈值的光谱图像分割和硬度图像重建4个部分组成。其中，硬度和光谱信息获取已在前面讨论过。建模部分主要包括在特征降维基础上，采用回归算法建立硬度预测模型，并在分析评价后确定最优模型。基于指数阈值的光谱图像分割是利用建立的新指数对光谱图像进行分割。最后，利用确定的最优硬度回归模型对分割出的草莓果肉部分进行硬度图像重建。

Figure 3.  Block diagram of the methodology for this article

高光谱相机获取的图像有616个通道光谱数据，为了减少数据冗余，提高处理效率，需对采集到的光谱信息降维^[17]。选取SPA和PCA及这两种算法组合的方式对采集到的高光谱图像进行降维处理。

SPA是一种前向迭代搜索算法，利用波长投影原理和交叉验证均方根误差(root mean square error of cross-validation, RMSECV)选出不共线的特征波长^[18]。以草莓光谱数据作为自变量，硬度数据作为因变量，选出与硬度相关的不共线光谱特征。

PCA将原来的波段信息通过协方差矩阵重组，产生新的特征向量，尽可能多地凸显原光谱信息^[19]。利用PCA算法将草莓光谱数据降维，选择贡献度高的特征作为输出结果。

采用SPA对裁剪后的光谱降维后，再利用PCA法对选择的光谱信息进行2次降维，即SPA-PCA的2次降维方法，进一步降低数据冗余，提取光谱域特征，以提高模型的稳定性。

将降维后高光谱数据与回归算法结合建立草莓硬度模型，不同的回归算法会直接影响模型的优劣^[20]。将457个草莓样本以4 ∶1的比例随机划分为训练集和测试集，其中训练集366个样本，测试集91个样本。选择偏最小二乘(partial least squares, PLS)法和随机森林(random forest, RF)法对草莓硬度建模。

PLS是一种多自变量对单因变量或多因变量的回归建模方法，已应用于高光谱技术对目标的无损检测中^[21]。RF使用决策树作为弱学习器，在样本特征维度很高的时候，具有较好的泛化能力。利用高光谱遥感数据和RF可对农作物生物量与产量进行估测^[22]。

为了评估模型的性能，引入5个评价指标来衡量模型效果^[23]，分别是校正集相关系数(correlation coe-fficient of calibration, CCC)、校正集均方根误差(root mean square error of calibration, RMSEC)、预测集相关系数(correlation coefficient of prediction, CCP)、预测集均方根误差(root mean square error of prediction, RMSEP)以及校正集相对分析误差(relative percentage difference, RPD)。

硬度分布重建需要分割掉草莓果肉以外的组成部分，在分析各组分光谱分布特征的基础上提取特征波段，依据特征波段建立光谱指数用于高光谱图像分割。

图 4为各组成部分的平均反射率曲线，不同类型组分的光谱差异明显。背景部分在全光谱范围反射率最低，低于0.08。正常果肉部分总体反射率较高，在656 nm和693 nm处存在波峰和波谷。萼片和草莓籽在波长693 nm处均存在波谷，且在693 nm~760 nm范围内反射率急速上升，具有明显的红边效应，可见两个部分均含有叶绿素，两者的波峰存在差异，草莓籽的波峰位于640 nm，而萼片波峰处于515 nm。霉变果肉部分的光谱反射率总体上随着波长的增加而增加，且整体斜率较低，并失去了波峰和波谷。因此，选择将萼片波峰、草莓籽波峰及红边波谷的位置作为光谱特征波段，分别是515 nm、640 nm和693 nm。

Figure 4.  Spectral feature curves of different components

为了增强特征波段之间的差异，提高识别效果，利用归一化植被指数^[24-27]来对草莓果肉进行提取。归一化指数可消除光照不均对反射率的影响，根据选出的草莓特征波长，构建新的归一化指数I₁、I₂和I₃：

式中: R₅₁₅、R₆₄₀和R₆₉₃分别代表515 nm、640 nm和693 nm处的光谱反射率。

由于5种组成部分反射率曲线在900 nm~1000 nm范围内分布平坦稳定，且各组分间分布差异较大，背景反射率明显低于其它组分，因此，选择960 nm波段为背景分割波长。利用阈值逐层分割的方法对不同组分进行分割，其分割阈值由960 nm反射率及3个归一化指数确定。

首先，建立像素数据集Q，记录各像素点位置信息并计算出用于分割的各参数值。

式中：n为像素总数; q_i记录第i个像素点位置、参数和标签信息; (x_i, y_i)和V_i分别为第i个像素点的坐标和用于组分分类的参数值; l为此像素点待标标签，取值为0、1、2、3和4，分别代表背景、萼片、果肉、发霉果肉和草莓籽。

式中: R₉₆₀和I₁、I₂、I₃分别代表第i像素点960 nm处的反射率及3个归一化指数。

利用多阈值方法对各像素点进行类别划分，并根据分类结果将标签赋值。阈值为：

式中: C为阈值的集合; θ₀为用于960 nm处的反射率阈值，由前面各组分光谱分析可知，背景光谱反射率很低，因此θ₀取为0.1;θ₁~θ₃为各新建的归一化指数阈值，其值根据各组分在箱线图上的分布情况由统计学方法确定，具体计算方法如下：

式中: Q_{l, 1}为类别1下四分位数; Q_{u, 2}为类别2上四分位数; I_{IQR, 1}、I_{IQR, 2}为两个类别四分位距(inter quartile range, IQR)。

具体逐层分割方法如下：首先通过判断像素点q_i的R₉₆₀与阈值θ₀的关系，判断是否为背景，即R₉₆₀≤θ₀则l=0为背景，否则为草莓的有效组成部分; 然后通过I₁~I₃及其阈值θ₁~θ₃逐层确定像素点的类别，即通过判断I₁≤θ₁、I₂≥θ₂、I₃≥θ₃分别给l赋值为1、2、3，若上述条件均不满足，则将l赋值为4，这样对数据集Q中每个像素点完成标签赋值，实现不同组分的分割。

根据阈值，将背景之外的其它各组分数据单独提取，生成组分-背景二值化图，分离不同组分。

在对高光谱图像完成分割后，剔除掉无关组分，利用已建的最佳硬度预测模型, 将草莓果肉部分每一个像素点的光谱数据进行硬度拟合，得到其灰度图，再利用伪彩色处理将不同硬度附上不同的颜色，绘制出草莓的硬度分布图像。

SPA根据RMSECV确定选择的特征波长数和光谱具体波段，当RMSECV达到明显拐点时，此时所对应的特征波长数即为降维后的结果。由图 5a可知，RMSECV随着选择特征波长数量的增加整体上呈下降趋势，当达到拐点时，值为0.1168，对应选择的特征波长数为11。实验中所选400 mm~1167 nm共524个波段，图 5b为SPA降维后11个特征波长数在此524个波段中的具体位置。

Figure 5.  Feature wavelength selection based on SPA

作为一种无监督的机器学习算法，PCA可由原始信息产生新的特征，新特征对于原光谱数据的总解释贡献度如表 1所示。当新特征数为9时，贡献度达到0.9997，因此，选择前9个特征作为降维后的新特征。

利用SPA-PCA 2次降维后产生的新特征贡献度如表 2所示。当特征数达到7时，贡献度为0.9997，因此，选择前7个特征作为SPA-PCA降维后的特征波段。

3种降维方法提取的特征个数占原特征数量的百分比分别为1.79%、2.46%和2.14%，有效地减少了数据冗余。

将原始光谱(full)和3种降维后的数据分别与PLS和RF结合建立硬度预测模型，构建8种模型的评价指标。

由表 3和表 4可知，在使用3种降维方法对特征进行提取后，PLS在校正集上的CCP均有不同程度的提升，其中基于SPA-PCA的2次降维效果最好，由全光谱的0.8636提高到0.9099。但是使用降维算法后，RF在预测集上的CCP有不同程度的下降，其中SPA-PCA的2次降维下降最多，由全光谱的0.8889降低到0.8158。因此，选择SPA-PCA-PLS和full-RF两个回归算法模型做进一步比较。

两者的校正集和预测集样本模型预测值和硬度计测量值分布散点图如图 6所示。

Figure 6.  Hardness model based on optimal PLS and RF

可以看出，SPA-PCA-PLS模型的预测效果更好。其次，full-RF模型未对光谱特征进行降维，存在着大量冗余的信息，而SPA-PCA-PLS模型确定的特征个数仅有7个，其中RPD也是所有模型中最高的，达到了3.2623，因此，确定SPA-PCA-PLS为硬度预测模型，建立硬度预测公式为：

式中: H_p为硬度的预测值; F_i为2次降维后的特征变量(i=1, 2, 3, 4, 5, 6, 7)。

将4个类型504个样本(草莓果肉、发霉果肉、萼片、草莓籽，每类126个)的3个参数(I₁、I₂和I₃)计算出来并分别绘制箱线图。箱线图用来比较各组分在3个参数中的分布情况，阈值由上下四分位数及四分位距来确定，结果如图 7所示。

Figure 7.  Distribution of the sample by index

由图 7a、图 7b可知，草莓果肉和萼片与另外两种组分在I₁和I₂上的分布有明显差异。I₁上发霉果肉及萼片异常值更少，对萼片的识别效果更好。I₁中草莓籽下四分位数Q_{l, 1}与四分位距I_{IQR, 1}分别为0.322、0.122;萼片上四分位数Q_{u, 2}与四分位距I_{IQR, 2}分别为-0.307、0.078，确定θ₁分割阈值为0.026。草莓果肉1.5I_{IQR, 1}和发霉果肉1.5I_{IQR, 2}在I₁上有明显重叠，在I₂上无重叠现象，确定I₂对草莓果肉的识别效果更好。I₂中草莓果肉下四分位数Q_{l, 1}与四分位距I_{IQR, 1}分别为0.807、0.118，发霉果肉上四分位数Q_{u, 2}与四分位距I_{IQR, 2}分别为0.542、0.137，则θ₂分割阈值为0.689。由图 7c可知，I₃可将4种类型分为两个大类，分别为正常果肉和发霉果肉、萼片和草莓籽, I₃中草莓籽下四分位数Q_{l, 1}与四分位距I_{IQR, 1}分别为0.118、0.075，发霉果肉上四分位数Q_{u, 2}与四分位距I_{IQR, 2}分别为-0.066、0.039，确定θ₃分割阈值为0.001。

由光谱特征分析，确定960 nm分割阈值为0.1对背景进行分割，则θ₀~θ₃分别为0.1、0.026、0.689，0.001。

在对高光谱图像中每个像素点经多阈值的逐层分割后，所分割出不同组分的二值化结果，如图 8所示。

Figure 8.  Component-background binarization results

图 8a为第3天的草莓样本RGB图; 图 8b为分割出的草莓区域图，轮廓清晰；图 8c和图 8d显示萼片和草莓籽部分，部分萼片和果梗被错分成草莓籽，原因在于萼片和果梗均含有叶绿素，两者的光谱反射率在515 nm~640 nm范围内平稳，并没有突出的波峰，光谱曲线与草莓籽相近; 图 8e为分割出的霉变果肉，在存放3天后，草莓均已经开始不同程度腐败和霉变，红色框霉变程度较明显，由点状霉变发展成片状霉变；图 8f为去除掉无关组分的草莓果肉，用于后续果肉硬度重建。

提取草莓果肉部分的光谱信息，基于SPA-PCA-PLS最优模型所确定的硬度拟合公式，对每个像素点进行硬度拟合并伪彩色化，图 9为不同贮藏天数的草莓硬度分布图。图中①~⑥为6个草莓样本，day 1~day 5表示不同的贮藏时间。硬度值最低为0 Pa，最高3×10⁵ Pa，由蓝到红硬度不断提高。总体上看，草莓硬度值经历了先下降后上升的两个过程。采摘时草莓成熟度不高，随着存放时间的增加，草莓不断成熟，细胞壁变薄，果胶质逐渐水解为果胶酸，细胞坚固程度下降，草莓硬度值不断下降，到第3天硬度达到最低。等到草莓完全成熟后保水能力变差，在室内环境不断进行蒸腾作用，水分减少，果实失水皱缩，体积变小，导致了硬度上升，从第4天开始硬度增大，且高于刚采摘后的草莓。草莓样本④在第3天出现了霉变，霉变发生在右下角区域，随着时间的变化，霉变区域不断扩大，第5天时已几近于覆盖整颗草莓。

Figure 9.  Strawberry hardness reconstruction image

利用高光谱成像系统采集草莓的高光谱图像，在SPA、PCA及2次组合降维的基础上，分别建立基于PLS和RF的硬度预测模型。构建515 nm、640 nm、693 nm波段归一化反射率指数，采用多阈值的逐层分割方法完成草莓组分分割，并实现草莓果肉部分硬度重建和直观显示。

为了提高模型的稳健性和准确性，分析了3种降维方法和两种回归算法的模型建立方法，选择最优的SPA-PCA-PLS模型作为硬度预测模型，CCP和RPD分别提升至0.9099、3.2623，RMSEP降至0.1344。除了模型建立方法外，样本量也决定了预测效果的好坏，实验中选用457颗草莓作为样本, 达到了比较好的预测效果，但由于样本个数的限制和同谱异物的影响，可能导致模型的泛化能力较差，例如在不同品种及不同产地上草莓的硬度预测效果有所差异。在本研究的基础上，更多样本数量的加入有望提高模型的泛化能力，以增加所提出方法的适用性。