三文鱼样本来自法罗群岛当地最大的三文鱼养殖公司Bakkafrost。 三文鱼被捕捞后, 由加工厂在3 h内将鲜活三文鱼屠宰加工为冰鲜三文鱼。 然后全程4 ℃冷链空运至上海纯尔贸易发展有限公司, 再通过冷藏车运送至镇江麦德龙超市, 真空包装后, 置放在铺满碎冰的保鲜盒中快速运送到实验室, 整个冷链过程约为4 d。

三文鱼去皮清洗干净后, 用灭菌后的刀片取中段鱼块切成长3 cm、 宽3 cm、 高2 cm的鱼片, 个体质量约(20± 5) g, 真空包装并贴上标签。 将分割好的三文鱼样本分为冰鲜组(C)和冻融组(FT), 冻融组又分为一次冻融组(FT-1)、 两次冻融组(FT-2)以及三次冻融组(FT-3), 每组30个样本, 共120个。 参考Ali^[7]等的循环冻融方法, C组真空包装后当天进行分析。 FT-1组在-30 ℃下冷冻30 d后解冻分析。 FT-2组先在-30 ℃下冷冻15 d, 后取出在4 ℃下解冻12 h, 然后继续保藏在-30 ℃条件下, 第30天取出, 解冻后分析。 FT-3组在在第10天和第20天分别完成一次冻融, 在第30天取出分析。

采用的高光谱成像系统主要包括高光谱相机(ImSpector V10E, Spectral Imaging Ltd, Oulu, Finland), 石英卤素灯(Fiber-LiteDC-950Illuminator, Dolan-Jenner Industries Inc, America), 精密自动平移装置(SC30021 A, Zolix Instruments co. Ltd., China), 以及装配了图像采集卡的电脑(P4P800-X型号, Asus computer co. Ltd., Taiwan, China), 其中高光谱相机由一台光谱仪和一台CCD相机组成, 光谱仪是整个高光谱成像采集系统最核心的部件, 能够在测试时获取样品面每个点在431~962 nm波段的光谱信息。 系统硬件部分主要组成示意图见图1。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 高光谱成像系统Fig.1 Hyperspectral imaging system

1.3.1 L^* , a^* , b^* 测定

仪器开机预热之后进行黑白板校正, 校正后开始对样品颜色进行采集, 每个样本进行三次颜色测定, 取三次结果的平均值作为样品的L^* , a^* , b^* 的真值。

1.3.2 剪切力测定

取3 cm× 3 cm× 2 cm大小的三文鱼样品, 用滤纸吸除样品表面水分后, 利用TA-XY2型物性测量仪对样品剪切力进行测定, 选用道具型号为warner-bratzler blade, 选用模式为compression, 剪切速度均为60 mm· min^-1, 目标模式设置为distance, 鱼片按肌纤维垂直于刀口方向放置于刀槽上。 仪器测定过程中刀具受到的最大剪切力值即样品的剪切力测量值。

1.3.3 K值的测定

(1) K值的计算方法

鱼死后, 机体内的三磷酸腺苷(Adenosine-5’ -triphosphate, ATP)会按顺序分解为二磷酸腺苷(Adenosine 5’ -diphosphate, ADP) 、 一磷酸腺苷(Adenosine monophosphate, AMP)、 肌苷酸(Inosincacid, inosinemonphosphate, IMP)、 肌苷(Inosine, HxR)、 次黄嘌呤(Hypoxhantine, Hx)。 K值为鱼肉中HxR和Hx总量之和与ATP, ADP, AMP, IMP, HxR和Hx总量之和的比值, K值越大说明鱼新鲜度越差。 计算公式如下^[5]

K=MHxR+MHxMATP+MADP+MAMP+MIMP+MHxR+MHx×100%

式中: M_ATP, M_ADP, M_AMP, M_IMP, M_HxR和M_Hx分别为样品中ATP, ADP, AMP, IMP, HxR和Hx的含量, 单位均为μ mol· g^-1。

(2) 核苷酸关联物的测定

采用水产标准SC/T 3048— 2014的方法对K值进行测定。

1.4.1 样品光谱提取与预处理

采集的高光谱图像是一个三维数据块, 不仅包含着样品的光谱信息、 图像信息, 同时也含有很多背景光谱。 因此合理的选取一块有效的区域进行研究以及对提取的光谱进行预处理对减少数据处理运算量和提高模型的正确率有至关重要的作用^[8]。 采用ENVI 4.5软件中的矩形工具选取样品中心附近200 pixel× 200 pixel的范围作为感兴趣区域(ROI), ROI选定后, 对ROI范围内所有像素点的光谱值进行平均, 作为该样本的光谱值, 每个样本对应一条光谱值。 采用一阶导数(1^st Der)、 二阶导数(2^nd Der)、 多元散射校正(MSC)、 标准正态变量变换(SNVT)、 归一化(N)和均值中心化(MC)来减少光谱重现性差、 噪声高等不良影响。

1.4.2 特征波长选择

采集的的光谱信息包含618个波长的光谱值, 数据繁多, 其中可能会存在冗余信息和噪声, 若将每个光谱值都输入模型进行分析, 不仅影响模型的准确率, 也会增大模型的运算量, 降低模型运算速度, 难以满足光谱检测的快捷性^[9], 因此有必要通过变量筛选的方法对光谱进行降维和特征波长筛选。 故采用竞争性自适应重加权算法(competitive adaptive reweighted sampling, CARS)^[10]、 区间变量迭代空间收缩法(interval variable iterative space shrinkage approach, iVISSA)^{[11, 12]}以及iVISSA-CARS^[13]筛选特征变量, 并比较结果, 从而确定三个品质指标最优的变量选择方法。

对冰鲜与冻融三文鱼的颜色测定如表1所示, 冰鲜三文鱼具有较高的亮度值L^* 和红度值a^* , 这是因为冰鲜的三文鱼表面偏红并具有光泽, 肉色较好。 冻融三文鱼随着冻融次数的增加, 亮度值和红度值均逐渐降低, 可能是由于三文鱼肉中的虾青素和胡萝卜素等物质的降解以及随汁液流失造成的。 三文鱼的黄度值在冻融后呈现上升趋势, 可能是由于冷冻期间鱼肉表面的脂肪成分发生了酸败与氧化, 从而呈现了淡黄色^[14]。 根据颜色差异显著性分析得到不同冻融次数三文鱼的红度值a^* 之间差异较显著, 因此在后续的颜色品质检测中选用红度值a^* 作为颜色指标进行检测。

表1

Table 1

表1(Table 1)

表1 冰鲜与冻融三文鱼的L* , a* , b* 值 Table 1 L* , a* , b* values of chilled and frozen-thawed salmon 样品组 L* a* b* C 53.66± 2.72a 16.96± 0.96a 17.64± 0.65c FT-1 52.61± 2.53a 15.24± 0.80b 18.01± 1.19bc FT-2 48.63± 2.99b 14.42± 0.97c 18.41± 0.82b FT-3 48.10± 2.13b 12.09± 0.48d 18.92± 0.65a 注: 小写字母不同表示在p< 0.05水平差异显著 Note: different lower case letters show significant differences (p< 0.05)

表1 冰鲜与冻融三文鱼的L* , a* , b* 值 Table 1 L* , a* , b* values of chilled and frozen-thawed salmon

剪切力的测定结果如图2黑色柱状图所示, 从图中可以看出, 冰鲜三文鱼的剪切力高于三组冻融样本, 并且冻融样本的剪切力随着冻融次数的增加呈现下降趋势。 其原因一方面可能是在冷冻过程中, 鱼肉组织受到自身机体内冰晶的破坏, 造成肌肉纤维间隙增大、 膜组织破裂等机械损伤; 另一方面, 低温冻结过程中肌肉组织可能存在冷冻不均匀现象从而产生内应力。 受内应力的影响, 肌纤维会发生低温断裂, 从而造成肉品的机械强度降低, 剪切力下降。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 冰鲜与冻融三文鱼的剪切力和K值变化Fig.2 Changes in shear force and K values of chilled and frozen-thawed salmon

冰鲜与不同冻融次数三文鱼K值测量结果如图2红色柱状图所示, 由图中可以看出, 从市场上购入的冰鲜三文鱼即C组的初始K值在11%左右, 处于一级新鲜范围; 样本冷冻30 d后, K值有所上升, 说明在冷冻过程中样本的新鲜度在下降; 样本经过三次冻融循环后, 其K值已经超过了20%, 变为二级新鲜度, 说明多次冻融会加速样品的变质。

采用SPXY方法对样本集划分, 划分结果见表2。 由表中可以看出, 基于SPXY方法划分的校正集和预测集的品质指标较均匀, 并且预测集的各品质指标值均在校正集的品质指标值范围内, 有助于定量模型精度提升^[15]。

表2

Table 2

表2(Table 2)

表2 校正集和测试集的划分结果 Table 2 Division results of calibration set and test set 品质指标 校正集 预测集 含量范围 平均 值 标准 偏差 含量范围 平均 值 标准 偏差 a* 11.21~18.52 14.94 1.95 11.51~18.03 14.16 1.84 剪切力/N 4.79~10.98 7.46 1.58 5.21~8.66 6.63 1.00 K值 9.71~25.89 16.66 4.49 10.48~25.88 16.86 4.03

表2 校正集和测试集的划分结果 Table 2 Division results of calibration set and test set

样本集划分后, 建立了基于全光谱波段范围内的原始光谱和不同预处理后光谱的a^* 、 剪切力、 K值PLS预测模型, 并对模型结果进行了分析和比较, 结果见表3。 由表中结果可以看出, 采用一阶导数预处理, 基于全光谱变量的PLS模型对a^* 预测效果最好, 此时的R_c和R_p值分别为0.926 9和0.913 5, RMSECV和RMSEP分别为0.72和1.00。 经过二阶导数预处理后的全光谱变量所建模型对剪切力的预测效果最优, 其R_c和R_p分别为0.874 6和0.842 1, RMSECV和RMSEP分别为0.68 N和0.82 N。 在K值的预测模型中, 经过归一化(N)预处理的光谱所建模型效果最优, 其R_c和R_p分别为0.881 7和0.846 3, RMSECV和RMSEP分别为1.75和2.27。 为尝试提高模型的预测效果和降低模型输入变量的复杂性, 分别采用CARS, iVISSA以及iVISSA-CARS方法对光谱变量进行筛选, 并将筛选后的光谱值作为变量建立a^* 值、 剪切力和K值的PLS定量模型, 比较模型的预测效果。

表3

Table 3

表3(Table 3)

表3 基于全光谱变量的三文鱼a* 、 剪切力、 K值预测结果 Table 3 Prediction results of a* , shear force and K value based on full spectrum variables 品质 指标 预处理 方法 Rc RMSECV Rp RMSEP a* 值 Raw 0.914 1 0.94 0.865 6 1.23 1st Der 0.926 9 0.72 0.913 5 1.00 2nd Der 0.789 8 1.35 0.755 9 1.42 MSC 0.846 6 1.27 0.792 6 1.31 SNVT 0.761 7 1.42 0.707 6 1.43 N 0.919 9 0.89 0.852 5 1.20 MC 0.899 7 1.14 0.882 5 1.18 剪切力 Raw 0.862 6 0.69 N 0.796 6 0.73 N 1st Der 0.851 6 0.73 N 0.826 2 0.85 N 2nd Der 0.874 6 0.68 N 0.842 1 0.82 N MSC 0.703 2 1.12 N 0.569 2 1.55 N SNVT 0.767 4 1.02 N 0.618 3 1.46 N N 0.713 0 1.04 N 0.676 1 1.23 N MC 0.536 2 1.94 N 0.478 9 2.02 N K值 Raw 0.851 4 2.17 0.810 2 2.65 1st Der 0.853 8 2.02 0.840 7 2.38 2nd Der 0.859 3 1.94 0.774 2 2.83 MSC 0.860 4 1.93 0.831 2 2.45 SNVT 0.851 0 2.24 0.802 5 2.74 N 0.881 7 1.75 0.846 3 2.27 MC 0.852 6 2.13 0.823 1 2.54

表3 基于全光谱变量的三文鱼a* 、 剪切力、 K值预测结果 Table 3 Prediction results of a* , shear force and K value based on full spectrum variables

使用iVISSA在431~962 nm波段范围内筛选出14个波段[见图3(a)]共217个与a^* 值相关的特征变量, 筛选出13个与剪切力、 K值相关的波段以及对应的443个和324个特征波长点[见图3(b, c)]。 将iVISSA筛选出的变量与对应的三文鱼品质指标建立预测模型。 iVISSA-PLS对a^* 值预测模型的R_c和R_p分别为0.906 6和0.888 3, RMSECV和RMSEP分别为0.80和0.99。 预测结果稍低于基于全光谱变量的PLS模型, 但变量简化了, 降低了模型的运算复杂度。 iVISSA-PLS对剪切力预测模型的R_c和R_p分别为0.650 5和0.644 4, RMSECV和RMSEP分别为1.10 N和1.22 N, 虽然相对于全光谱建立的模型来说变量数减少了28.31%, 但是R_p也降低了23.47%。 iVISSA-PLS对K值预测模型的R_c和R_p分别为0.878 4, 0.878 0, RMSECV和RMSEP分别为2.07, 2.29, 相对于全光谱建立的模型, iVISSA-PLS对K值的预测效果得到提高。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 iVISSA挑选的与a* (a)、 剪切力(b)以及K值(c)相关的波长区间位置Fig.3 Wavelength interval positions selected by iVISSA for a* (a), shear force (b) and K value (c)

在简化a^* 值预测PLS模型过程中, 使用CARS变量筛选方法在431~962 nm波段范围内共筛选得到51个与a^* 值相关的特征波长, 将120个样本的51个波长组合为新的120× 51的变量矩阵, 建立三文鱼a^* 值的PLS预测模型, 模型结果如图4(a)所示。 由图4(a)结果可以看出, PLS对a^* 值预测模型的R_c和R_p分别为0.931 6和0.929 7, RMSECV和RMSEP分别为0.72和0.74。 CARS挑选出61个与剪切力相关的特征变量, 用这61个特征变量建立的PLS预测模型得到较好的剪切力预测效果, R_c和R_p分别为0.892 1, 0.887 3, RMSECV和RMSEP分别为0.67 N和0.80 N, 见图4(b)。 在简化K值预测模型中, CARS从全光谱中筛选出了51个特征变量。 建立的CARS-PLS模型对K值预测结果如图4(c)所示, R_c=0.951 4, R_p=0.950 0, RMSECV=1.33, RMSEP=1.53。 CARS-PLS对三个品质指标的预测结果与各自的基于全光谱变量的PLS模型效果相比有较明显的提升。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 CARS-PLS模型对a* (a)、 剪切力(b)以及K值(c)最佳预测结果Fig.4 The best prediction results of a* (a), shear force (b) and K value (c) obtained by CARS-PLS model

为了更简化模型数据量, 提高模型的预测性能, 将波段选择方法与波长点选择方法结合, 先通过iVISSA波段选择方法对特征变量进行粗选, 然后使用CARS精选波长点。 图5分别展示了CARS对三个品质指标特征变量筛选过程, 在iVISSA筛选的217个变量的基础上, 筛选了31个与a^* 相关的特征变量[见图5(a)], 以这些变量建立的对a^* 最佳PLS预测模型的R_c和R_p分别为0.915 4和0.906 3, RMSECV和RMSEP分别为0.79, 0.83。 使用该变量选择方法提取的变量数最少, 极大的简化了模型的数据量, 并且预测效果高于iVISSA-PLS, 稍低于全光谱建立的PLS及CARS-PLS。 在iVISSA提取了443个与剪切力相关的变量后, CARS经过波长点精选出38个变量[见图5(b)], 以这些变量建立的对剪切力最佳PLS预测模型的R_c和R_p分别为0.885 3和0.860 9, RMSECV和RMSEP分别为0.69 N和0.90 N。 相对于全光谱建立的PLS模型, iVISSA-CARS-PLS模型的变量数从618减少到38, 而R_p从0.842 1增加到0.860 9; 相对于iVISSA-PLS, iVISSA-CARS-PLS模型的预测性能得到显著提高; 但相对于CARS-PLS, iVISSA-CARS-PLS模型的变量数减少37.70%, R_p为CARS-PLS的97.02%。 经过iVISSA-CARS筛选后, 31个变量用于建立预测K值的PLS模型[见图5(c)], 以该变量建立的对K值的最佳PLS预测模型的R_c和R_p分别为0.950 7和0.940 4, RMSECV和RMSEP分别为1.35和1.56。 相对于全光谱建立的PLS模型以及iVISSA-PLS模型, iVISSA-CARS-PLS采用的变量数大大较少, 并且R_c和R_p也分别提高了11.1%与7.11%; 虽然iVISSA-CARS-PLS得到的R_p为CARS-PLS的98.98%, 但是采用的变量数减少了39.21%。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 CARS选择与a* (a)、 剪切力(b)以及K值(c)相关特征变量的过程Fig.5 The processes of CARS selecting characteristic variables related to a* (a), shear force (b) and K value (c)

对上述三个品质指标的不同算法预测结果(见表4)比较得出, iVISSA-CARS算法筛选得到的变量数最少, 较大程度减少了模型的数据量, 但CARS筛选后的光谱作为变量建立的PLS模型预测效果最佳, 其测试集具有较高的相关系数以及最低的RMSEP值, 校正集的模型效果也高于其余三种模型。 同时R_c与R_p之间差值最小, 说明模型不仅具有较好的预测性能, 而且具有较高的一致性, 并且相对于全光谱建立的模型来说极大的减少了变量数, 说明CARS变量选择方法适合运用于三文鱼品质指标特征变量的提取。 同时从表4中可以看出, 三个指标的CARS-PLS以及iVISSA-CARS-PLS模型均高于iVISSA-PLS, 说明对于三文鱼上述三个品质指标的预测, CARS的特征变量筛选策略优于iVISSA的波段选择策略。 因此选择CARS-PLS用于后续品质指标的可视化分析中。

表4

Table 4

表4(Table 4)

表4 基于波长筛选后建立的a* 、 剪切力、 K值预测模型结果 Table 4 The model results of a* , shear force and K value prediction based on screened wavelength 品质指标 模型 变量数 Rc RMSECV Rp RMSEP CARS-PLS 51 0.931 6 0.72 0.929 7 0.74 a* iVISSA-PLS 217 0.906 6 0.80 0.888 3 0.99 iVISSA-CARS-PLS 31 0.915 4 0.79 0.906 3 0.83 CARS-PLS 61 0.892 1 0.67 N 0.887 3 0.80 N 剪切力 iVISSA-PLS 443 0.650 5 1.10 N 0.644 4 1.22 N iVISSA-CARS-PLS 38 0.885 3 0.69 N 0.860 9 0.90 N CARS-PLS 51 0.951 4 1.33 0.950 0 1.53 K值 iVISSA-PLS 324 0.878 4 2.07 0.878 0 2.29 iVISSA-CARS-PLS 31 0.950 7 1.35 0.940 4 1.56

表4 基于波长筛选后建立的a* 、 剪切力、 K值预测模型结果 Table 4 The model results of a* , shear force and K value prediction based on screened wavelength

样品生成的高光谱图像中的每个像素点对应一条光谱, 可以通过所建立的品质指标的预测模型的相关系数将光谱值转化为品质指标的预测值, 再将不同颜色表示的预测值以像素的形式叠加在原始图像中, 即可得到代表指标数值大小以及空间分布的伪彩图^[16]。

图6(A)为a^* 值的分布可视图, C组样本的a^* 值分布图中黄色像素点居多, 代表着样本具有较高的a^* 值; 冻融组FT-1, FT-2和FT-3样品的a^* 值分布图中黄色像素点个数逐渐减少, 蓝色像素点增多, 说明冻融组的a^* 值低于冰鲜组, 与表2的结果相符, 其中FT-3组样本a^* 分布图中的蓝色像素点最多, 说明三文鱼经历三次冻融后a^* 值明显下降。 图6(B)为冰鲜与冻融三文鱼的剪切力分布图, 通过剪切力可视化分布图可以直观的看出, 冰鲜三文鱼具有较高的剪切力; 经历冻融后, 三文鱼剪切力下降, 并且随冻融次数的增加, 剪切力呈下降趋势。 图6(C)为冰鲜与冻融三文鱼的K值分布图。 K值作为表征鱼类新鲜度的重要指标, 对其实时监测具有重要意义。 通过三文鱼的K值分布图, 可以直观的了解三文鱼的新鲜度情况。 从可视化图中可以看出, C组三文鱼具有较低的K值, 冻融三文鱼在冻藏期间K值有所上升, 并且随着冻融次数的增加而加快上升。 以伪彩图的形式对三文鱼品质指标的可视化分析, 可以直观的展示品质指标的大小, 还能展示其空间分布, 有利于商家对三文鱼更加精确的切割管理以及消费者对鱼体品质的快速评估。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 (A), (B), (C)分别为a* 、 剪切力、 K值可视化分布图, (a), (b), (c), (d)分别代表C组、 FT-1、 FT-2、 FT-3组Fig.6 (A), (B), (C) are the visual distribution maps of a* , shear force, and K value, respectively, and (a), (b), (c), (d) represent group C, FT-1, FT-2, and FT-3