基于高光谱成像技术和力学参数对贡梨冲击损伤的定量研究

李 斌 邹吉平 张 烽 苏成涛 刘燕德 肖毅华

(1.华东交通大学 机电与车辆工程学院，南昌 330013； 2.华东交通大学 水果智能光电检测技术与设备国家与地方联合工程研究中心，南昌 330013)

贡梨富含人体所需的微量元素，深受消费者喜爱。但贡梨在采摘、分拣、包装、运输的过程容易受到冲击或挤压等因素的影响，使得贡梨产生机械损伤，而在机械损伤类型中，冲击损伤最容易产生[1]。损伤后的贡梨品质降低，进而造成大量的经济损失。为了减少经济损失，亟需一种准确评估和预测贡梨冲击损伤的方法。

目前，对于水果损伤检测的方法很多：1)电子鼻技术[2]；2)X射线计算机断层扫描[3]；3)核磁共振[4]；4)光谱检测技术[5]，其中光谱检测技术在水果品质检测方面具有很大的潜力。近年来，随着光谱技术的不断发展，作为新一代光电检测技术的高光谱成像技术广泛应用于农产品质量检测领域[6]。高光谱成像技术结合了光谱学、机器视觉、计算机图像学、近红外光谱等多学科知识，将光谱技术和传统的二维成像技术进行有机融合，可同时提供光谱信息和空间信息[7-8]。沈宇等[9]获取400～1 000 nm波长范围内完好和轻微损伤后不同时间段的苹果高光谱图像，并提取完好样本以及损伤样本的光谱数据。利用载荷系数法(x-LW)、连续投影法(Successive projections algorithm，SPA)以及二阶导数法提取特征波长，基于特征波长分别建立遗传算法优化的BP神经网络(GA-BP)和支持向量机(Support vector machines，SVM)损伤识别模型，对比后发现SVM是最佳判别苹果轻微损伤的模型。孟庆龙等[10]采集400～100 nm波长范围内无损伤猕猴桃和表面有疤痕的猕猴桃高光谱图像，利用多种图像处理方法实现对无损和表面有疤痕的猕猴桃进行识别，正确识别率分别为98.3%和95.0%。Yuan等[11]利用定量损伤仪获得5种不同损伤程度的灵武长枣，采用多种预处理方法建立偏最小二乘分析(Partial least squares discriminant analysis，PLS-DA)模型，结果发现MSC-PLS-DA模型是最优分类模型。又利用多种波段筛选算法对原始光谱和预处理后光谱进行特征波长的提取，建立PLS-DA分类模型，结果表明MSC-CARS-PLS-DA是最优分类模型，准确率可达89.52%。以上对水果的损伤的研究大多数是定性分析，且所建立的模型也只能区分水果是否受损。目前，也有关于水果损伤的定量研究。Zhang等[12]利用高光谱成像技术采集900～1 700 nm波长范围内高光谱图像以及压敏薄膜技术测得苹果跌落过程中的机械参数，对苹果冲击损伤进行了定量研究。Xu等[13]将900～1 700 nm波长范围内所采集的光谱数据与利用压敏薄膜技术所测得的机械参数相结合建立偏最小二乘回归(Partial least squares regression，PLSR)模型，实现了对机械参数的定量预测。文献[12]和[13]中的研究仅用了一种模型对力学参数进行预测，未进行多种模型对比分析，针对贡梨冲击损伤的定量研究未见报道。

本研究拟采用397.5～1 014.0 nm波长范围内的高光谱数据与力学参数相结合的方法，建立力学参数的不同预测模型，比较分析后选取最优预测模型，实现对贡梨冲击损伤的准确评估和定量预测，以期为贡梨的品质分选以及包装运输提供理论支持。

1 材料与方法

1.1 试验样本

贡梨购于江西省南昌市某水果市场，于-1～2 ℃的温度下贮藏。试验前将贡梨在室温为20 ℃，相对湿度(Relative humidity，RH)为40%的条件下放置24 h，从而降低试验中温度对贡梨损伤造成的影响。为了降低贡梨质量和曲率半径对贡梨冲击损伤试验的影响，共选择240个表面没有明显缺陷、形状规则的贡梨，赤道直径约为70～80 mm，质量约(240±10) g。

1.2 试验装置及损伤试验

本研究基于单摆原理设计贡梨碰撞试验装置，该装置主要由底座、垂直支撑杆、水平支撑杆、鱼线以及水果夹具等组成(图1)。通过对水平支撑杆、压力传感器(HZC-H1，测量范围0～100 kg，灵敏度2.00 mV/V，金属材质)以及量角器装置进行调节，从而确保贡梨撞击在压力传感器的中心。每次贡梨与压力传感器碰撞后，立即用手抓住贡梨，防止产生二次碰撞。在整个试验过程中，智能数据分析采集器(Data acquisition & signal processing, DASP)与压力传感器相连，可及时记录力响应过程，并将记录的数据传输至计算机中储存。利用相关软件对采集的数据进行计算分析，可获得峰值力、平均接触力等力学参数。

利用贡梨碰撞试验装置获取不同损伤程度的贡梨：将贡梨样本随机分为6组，每组40个。6组贡梨样本分别以释放角度α为35°、40°、45°、50°、55°和60°释放，释放角度α为图1中量角器装置中轴线和鱼线之间的夹角。由于受损的贡梨样本难以通过肉眼辨别，还需将损伤的贡梨样本在室温条件下放置24 h，使其产生褐变，以便后续数据的测量与计算。

α，释放角度，(°)。 Release angle. 1.底座；2.压力传感器；3.导轨；4.垂直支撑杆；5.量角器装置；6.水平支撑杆；7.鱼线；8.水果夹具；9.智能数据分析采集器；10.计算机 1.Baseplate; 2.Force transducer; 3. Guide rail; 4.Vertical support rod; 5.Angle measuring device; 6.Horizontal supporting rod; 7.Fishingline; 8.Fruit fixture; 9.Data acquisition & Signal processing; 10.Computer图1 贡梨碰撞试验装置Fig.1 Gongli impact test device

1.3 高光谱图像的采集及提取

高光谱成像系统的整体结构见图2，其硬件组成包括2盏20 W卤素灯(DECOSTAR51 MR16型，OSRAM，德国)、成像光谱仪(ImSpector V10E型，Specim，芬兰)、CCD相机(C8484-05G型，Hamamastu,日本)、位移平台、配有图像采集卡的计算机、暗箱等[14]。位移平台的移动通过步进电机来实现；光源、成像光谱仪、CCD相机、位移平台等硬件都安装在一个暗箱内。

1.步进电机；2.暗箱；3.CCD相机；4.成像光谱仪；5.镜头；6.卤素灯；7.位移平台；8.贡梨样本；9.计算机 1.Stepper motor; 2.Camera obscura; 3.CCD camera; 4.Imaging spectrometer; 5.Lens; 6.Halogen lamp; 7.Displacement platform; 8.Gongli sample; 9.Computer图2 高光谱成像系统示意图Fig.2 Diagram of the hyperspectral imaging system

在进行图像采集时，首先使高光谱成像系统预热30 min左右，以防止基线漂移影响数据精度。然后利用电脑和成像光谱仪配套的SpectraVIEW软件对高光谱成像系统的一些参数进行设置，将光谱分辨率设置为2.8 nm，CCD相机的曝光时间设置为20 ms，位移平台前进和后退的速度设置为15 mm/s[15]。最后将在室温下放置24 h后产生褐变的贡梨放置在移动载物台上进行高光谱图像采集。本研究中所采集的光谱数据的波长范围为397.5～1 014.0 nm，共176个波段。

由于高光谱图像易受暗电流和光照不均匀等因素的影响，在处理数据之前必须对高光谱图像进行黑白校准[16]。在校准时，首先利用镜头盖将镜头遮住，进行扫描得到全黑的校准图像，然后将镜头盖取下，扫描由聚四氟乙烯材料制成的白板从而获得全白的校准图像。校准公式[16]为：

(1)

式中：Rλ为校准后图像数据；RW为全白校准图像数据；RB为全黑校准图像数据；Iλ为原始图像数据。

利用ENVI4.5软件提取光谱数据。首先，通过SpectraVIEW软件将采集的高光谱图像进行校正，然后，将校正后的高光谱图像导入至ENVI4.5软件中，以贡梨受损区域为感兴趣区域(ROI)，提取从不同释放角度下贡梨损伤区域的光谱数据。

1.4 贡梨力学参数的测量及计算

1.4.1损伤面积

由于损伤过后贡梨在室温下放置24 h，损伤区域产生褐变。通过肉眼可以辨别损伤区域，方便对贡梨损伤区域进行测量。本试验中，将所有贡梨样本的损伤区域近似为圆形。用数显游标卡尺(精度为0.01 mm)测量损伤区域的直径，为确保试验数据的精度，测量损伤区域的直径时，每隔120°测量1次直径，共测量3次损伤区域直径，最后取平均值作为贡梨的损伤直径。损伤面积计算公式为：

(2)

式中：A为贡梨的损伤面积，mm2；D为贡梨损伤区域直径，mm。

1.4.2峰值力、平均接触力和平均压强

贡梨与压力传感器碰撞时，压力传感器的数据传输至DASP中，从而记录了力响应过程。峰值力为贡梨与压力传感器碰撞过程中力传感器所记录的最大值。贡梨与压力传感器碰撞过程中，从碰撞开始至碰撞结束这一时间段内，力变化的平均值为平均接触力，通过Coinv DASP V11软件对采集的数据进行时域分析，获得峰值力和平均接触力。

通过智能数据采集系统所获得的平均接触力以及所测量的损伤面积计算平均压强，公式为：

(3)

式中：σ为平均压强，MPa；F为平均接触力，N。

1.5 光谱预处理及特征波长的选择

由于仪器、样本、以及测量环境等条件的影响，所采集的原始光谱信息存在噪声和背景等干扰信息，这些干扰信息会造成定量分析出现误差、组分预测错误等问题[17-18]。因此，有必要对原始光谱进行预处理。光谱预处理方法有多种，如平滑、导数、多元散射校正(Multiplicative scatter correction，MSC)和标准正态变换(Standard normal variate,SNV)等[19]。本研究中采用Gap-segment求导、SG求导和基线校准(Baseline)3种预处理方法。

在提取的全光谱中存在部分无效或干扰信息，这些信息会降低模型的预测精度，同时全光谱中所含的数据量多，使得模型的运行时间增加。使用特征波长提取算法从全光谱中选择最具代表性的特征波长可以有效提高模型的检测效率。本研究采用竞争性自适应重加权(Competitive adapative reweighted sampling，CARS)和无信息变量消除(Uninformative variable elimination，UVE)2种算法选取特征波长。

CARS算法是基于蒙特卡洛采样与偏最小二乘(Partial least squares，PLS)回归系数的特征变量选择方法，该算法主要依据达尔文进化论中“适者生存”的原则。其核心思想是通过在采样过程中优先选择PLS回归过程中回归系数绝对值较大的变量。经过多次采样后，最终找出交叉验证均方根误差(RMSECV)最小的光谱数据子集作为最优变量组合[20]。

UVE算法的基本思想是将回归系数作为波长重要性的衡量指标。首先，将一定变量数目的随机变量矩阵引入光谱矩阵；然后，利用蒙卡罗交叉验证建立PLS模型；最后，计算PLS回归系数平均值与标准偏差的比值，进而选取有效光谱信息[21-22]。

1.6 模型建立及评价

试验样本共240个，均分成6组，利用Kennard-Stone(KS)方法，将样本按3∶1的比例分为建模集和预测集，即建模集180个样本，预测集60个样本[23]。分别建立峰值力、平均接触力、损伤面积和平均压强的PLSR和主成分回归(principal component regression, PCR)模型。

为检验模型精度，需要相关参数对所建立的模型进行评价。通过比较建模集相关系数(RC)、预测集相关系数(RP)、建模集均方根误差(RMSEC)和预测集均方根误差(RMSEP)来评价[24]。如果模型具有高的建模集相关系数(RC)与预测集相关系数(RP)，低的建模集均方根误差(RMSEC)和预测集均方根误差(RMSEP)，则表明该模型具有较好的性能。

2 结果与分析

2.1 力学参数测量结果分析

图3示出6种不同释放角度(α)下，贡梨碰撞力学参数的测量结果：随着贡梨释放角度的增大，峰值力、平均接触力、损伤面积和平均压强也在逐渐的增大。峰值力平均值从138.40 N增大至335.86 N，平均接触力平均值从77.13 N增大至188.20 N，损伤面积平均值从208.07 mm2增大至544.42 mm2，平均压强平均值从0.34 MPa增大至0.42 MPa。损伤面积可以直观的体现贡梨的损伤程度，随着贡梨释放角度的增大，损伤面积越大，冲击损伤越严重。损伤面积增大时，峰值力、平均接触力和平均压强也随之增大。可以看出，本研究中的力学参数存在一定的波动性，这可能是由于贡梨的果肉硬度、曲率半径和贡梨的内部结构差异等自身生理特性所导致。总体来说，随着释放角度的增大，力学参数呈现一定的线性变化。在一定程度上表明力学参数表征贡梨的冲击损伤程度是可行的。

柱上竖线为力学参数的标准偏差。 The vertical lines on the columns are the standard deviations of the mechanical parameters.图3 不同释放角度(α)下贡梨力学参数的变化Fig.3 Changes of mechanical parameters of Gongli under different release angles

2.2 光谱分析

图4 未受损及不同释放角度(α)下受损贡梨的光谱曲线Fig.4 Average spectral curves of damaged Gongli under undamaged and different release angles

未损伤和不同损伤程度的贡梨在波长397.5～1 014.0 nm范围内的光谱曲线见图4：在739 nm处存在较为明显的吸收峰，该吸收峰是由于碳水化合物中的C—H键伸缩振动引起的[25]。未损伤贡梨的光谱曲线与不同角度释放下受损贡梨的光谱曲线整体相似，但是未受损贡梨和不同释放角度下受损贡梨的光谱反射率在同一波段存在明显差异。随着释放角度的增大，贡梨损伤越严重，光谱反射率越低，产生这一变化的原因可能是贡梨与金属材质的压力传感器碰撞，使得贡梨果肉中完整的细胞壁和细胞膜被破坏，细胞破裂后，细胞中的水和酶被释放，从而造成光谱曲线上反射率的差异。从光谱曲线中反射率的差异可反映贡梨受损程度。因此本研究通过光谱数据与力学参数相结合的方法，进一步实现贡梨冲击损伤的定量预测。

2.3 光谱预处理

通过Gap-segment求导、SG求导和Baseline预处理过后的光谱数据和原始光谱数据分别建立PLSR和PCR模型，模型性能参数见表1和表2：采

表1 不同光谱预处理后力学参数的PLSR模型性能参数Table 1 Performance parameters of PLSR model with different spectral pretreatment

表2 不同光谱预处理后力学参数的PCR模型性能参数Table 2 Performance parameters of PCR model with different spectral pretreatment

用不同的预处理方法过后，力学参数的预测模型性能存在差异，其中采用Baseline预处理后，力学参数的PLSR和PCR模型性能最优。Gap-segment求导、SG求导预处理只能消除光谱数据中不重要的基线漂移，同时还会增大数据的噪声，降低信噪比。Baseline预处理消除光谱数据基线漂移的效果要比Gap-segment、SG求导预处理的效果好，因此，Baseline预处理是最佳预处理方法[26-27]。通过比较PLSR和PCR模型的RC、RP、RMSEC和RMSEP发现，对于力学参数的预测，PLSR模型的性能优于PCR模型的性能。PLSR模型对峰值力、平均接触力和损伤面积的预测效果较好，但是对于平均压强的预测效果并不理想。导致PLSR模型对平均压强预测效果差的原因可能是由于平均压强通过平均接触力除以损伤面积而获得，而贡梨本身的硬度、曲率半径、质量均存在差异，导致贡梨碰撞过程中力传感器所记录的数据以及损伤面积的测量均存在一定的误差，两者误差的累积导致平均压强不准确。综上所述，本研究采用Baseline光谱预处理并建立力学参数的PLSR模型。

2.4 特征波长的提取

利用CARS算法对Baseline预处理后的光谱数据进行特征波长的选取，以峰值力的特征波长选取为例，选取结果见图5(a)。随着采样次数的增加，被优选的波长变量个数逐渐的减少，在第一阶段变量个数减少的较快，随后逐渐减慢，表明使用CARS算法选取特征波长中有“粗选”和“精选”2个阶段(图5(a1))。当采样次数<31时，RMSECV值不断降低，表明波长中所含的干扰信息被剔除；当采样次数>31时，RMSECV值呈回升的趋势，表明波长中所含的有效信息被剔除，从而导致RMSECV值上升[28]；当采样次数为31时，RMSECV值最小(图5(a2))，即第31次采样得到的变量子集为最优的特征波长变量子集，其中该子集包含45个特征波长[29]。

图5 CARS和UVE算法选取峰值力的特征波长Fig.5 CARS and UVE algorithm selected characteristic wavelength results of peak force

利用UVE算法对Baseline预处理后的光谱数据进行特征波长的选取，以峰值力的特征波长选取为例，选取结果见图5(b)。由于UVE算法在选取特征波长时集噪声和浓度信息为一体，因此利用该算法对176个原始波长变量进行筛选时需要引入一定数量的随机变量，本研究中引入了176个随机变量，即原始波长变量和随机变量总个数为352个。分割线左右两侧分别为原始波长变量和随机变量的稳定性分布曲线，通过上、下阈值线进行特征变量的选取，阈值选定标准为随机变量稳定性最大值的99%。位于2条阈值线之间的波长变量被剔除，超出2条阈值线的波长变量被保留，用于进一步分析。经过UVE波段筛选过后，共有37个变量被选择。

利用CARS算法对峰值力、平均接触力、损伤面积和平均压强选取的特征波长个数分别为45、35、62和13个，选取的波长数量分别占全波长数量的25.6%、19.9%、35.2%和7.4%(表3)。UVE算法对峰值力、平均接触力、损伤面积和平均压强选取的特征波长个数分别为37、51、34和24个，选取的波长数量分别占全波长数量的21.0%、29.0%、19.3%和13.6%(表4)。从2个表中可以看出峰值力和平均接触力中含有大多数相同的波段，而损伤面积和平均压强中大多数特征波段并不相同，表明同一属性的力学参数所含的特征信息相同，不同属性的力学参数所含的特征信息不同。

2.5 基于特征波长和力学参数的PLSR建模

本研究利用CARS和UVE算法进行特征波长的选取，通过2种算法选取的特征波长与全波长为变量分别建立PLSR模型进行比较分析，模型的性能参数见表5。相比于Baseline-PLSR模型的预测效果，Baseline-UVE-PLSR模型对力学参数的预测效果较差。可能是采用UVE算法进行特征波长的选取过程中，剔除原始光谱信息中干扰信息的同时将一些有效信息也剔除了，从而导致对力学参数的预测效果变差。Baseline-CARS-PLSR模型对力学参数的预测效果最好，在 Baseline-CARS-PLSR模型中，峰值力、平均接触力、损伤面积和平均压强的RP和RMSEP分别为0.892和31.527 N、0.883和18.861 N、0.895和54.411 mm2、0.661和0.045 MPa。表明不同的特征波长选取方法对PLSR模型性能的影响不同，选择合适的特征波长选取方法可以有效的提高模型的预测精度。

表3 CARS算法选取力学参数的特征波长Table 3 Selection of characteristic wavelength of mechanical parameters by CARS algorithm

表4 UVE算法选取力学参数的特征波长Table 4 Selection of characteristic wavelength of mechanical parameters by UVE algorithm

表5 力学参数原始光谱与特征光谱的PLSR模型性能参数Table 5 Performance parameters of PLSR model for original and characteristic spectra of mechanical parameters

基于Baseline-CARS-PLSR模型力学参数的测量值与预测值见图6。分别对建模集和预测集的力学参数的测量值与预测值进行拟合后发现，峰值力、平均接触力和损伤面积的样本点较为均匀的分布

图6 基于Baseline-CARS-PLSR模型力学参数的测量值与预测值Fig.6 Measured and predicted values of mechanical parameters based on the Baseline-CARS-PLSR model

在拟合回归线周围，且样本点与拟合回归线的距离较近。峰值力、平均接触力和损伤面积的RP和RMSEP分别为0.892和31.527 N、0.883和18.861 N、0.895和54.411 mm2。综上，Baseline-CARS-PLSR模型对贡梨峰值力、平均接触力和损伤面积等力学参数的预测较为准确。

由图6(d)可知，对于平均压强的预测效果并不理想，Baseline-CARS-PLSR模型的RP和RMSEP分别为0.661和0.045 MPa。平均压强预测不理想的原因可能是平均压强的计算误差的累积或者是平均压强并不适用于表征贡梨冲击损伤程度。

3 结 论

为实现准确评估和预测贡梨的冲击损伤，本研究采用波长为397.5～1 014.0 nm的高光谱成像技术与力学参数相结合的方法对贡梨的冲击损伤进行定量研究。首先，对峰值力、平均接触力、损伤面积和平均压强等力学参数进行统计分析。然后，采用Gap-segment求导、SG求导和Baseline对原始光谱进行预处理，并建立力学参数的PLSR和PCR模型，从而选取最佳预处理方法及建模方法。最后，利用CARS和UVE算法进一步优化力学参数的预测模型，比较各模型的预测精度并确定力学参数的最佳预测模型，探索利用高光谱成像技术与力学参数相结合的方法预测贡梨冲击损伤的可行性，主要结论如下：

1)采用3种光谱预处理后，力学参数的PLSR和PCR模型性能存在差异，其中Baseline-PLSR模型对力学参数的预测效果最佳，表明选择合适的预处理方法以及建模方法可以提高模型的性能。

2)本研究在选取最佳的光谱预处理和建模方法后，利用2种算法选取的特征波长为变量建立力学参数的预测模型，其中 Baseline-CARS-PLSR模型对力学参数的预测效果最优，该模型中峰值力、平均接触力、损伤面积和平均压强的RP和RMSEP分别为0.892和31.527 N、0.883和18.861 N、0.895和54.411 mm2、0.661和0.045 MPa，表明CARS特征提取算法结合PLSR建模方法可以有效的提高模型的预测能力。高光谱成像技术与力学参数相结合的方法对贡梨冲击损伤进行定量预测具有较高的可行性。