本实验所用巴旦木为市场上购买, 实验中的巴旦木样品包含三类: 饱满的巴旦木、 内部缺损的巴旦木和空壳巴旦木。 选取饱满的巴旦木若干, 饱满的巴旦木无需任何处理, 内部缺损的巴旦木是用小刀截取饱满的巴旦木内部果仁来模拟不同饱满度的巴旦木, 空壳巴旦木是使用小刀和镊子取出饱满的巴旦木内部果仁来制备, 制备的巴旦木样品图片如图1(a)所示, 去除部分外壳巴旦木样品的图片如图1(b)所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 样品准备 (a): 巴旦木样品图; (b): 去除部分外壳巴旦木样品图Fig.1 Sample preparation (a) Picture of almond samples; (b) Picture of almond with partial shell removed

本实验采用太赫兹成像系统(中国青源峰达太赫兹科技有限公司QT-TO1000)对巴旦木样品进行扫描, 频率范围为0.1~4.0 THz, 最大扫描面积为100 mm×100 mm, 最大检测厚度为50 mm, 成像速度为60像素·s^{-1 [15]}。 图2(a)为QT-TO1000透射成像系统原理图, 飞秒激光器发射的飞秒脉冲经过分束器后分为泵浦光和探测光。 泵浦光入射到太赫兹发射器上以激发太赫兹脉冲。 然后太赫兹脉冲被聚焦并穿过样品。 携带样品信息的太赫兹脉冲与探测光共线, 到达探测器晶体并将信号传输到计算机, 在X-Y二维平移台上按0.2 mm的步长逐点扫描, 完成图像采集。 图2(b)为QT-TO1000仪器测试样品图, QT-TO1000系统在测量样品前需预热30分钟, 环境温度设置为(23± 0.5) ℃, 环境湿度设置为10%以下, 考虑到聚乙烯板(PE板)在THz频段基本无吸收^[16], 实验中使用厚度为1 mm的PE板作为样品的承托平台。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 QT-TO1000透射成像原理图(a)和仪器测试样品图(b)Fig.2 Transmission imaging schematic diagram (a) and whole machine (b) of QT-TO1000

由于提取无样品、 空壳巴旦木和满仁巴旦木感兴趣区域(regions of interest)的THz光谱, 根据Dorney等^[17]和Duvillaret等^[18]提出的光学参数提取模型。 采用快速傅里叶变换(FFT)获取了THz脉冲在频率上的频谱分布。 可表述为式(1), 其中, A(ω )表示电场幅值, φ (ω )为参考信号和样本信号的相位差, E(t)为太赫兹时域波形。 由式(2)和式(3)得到样品的折射率和吸收系数。

E(ω)=A(ω)e-iφ(ω)=∫E(t)e-iφ(ω)dt(1)

n(ω)=φ(ω)cωd+1(2)

α(ω)=2k(ω)ωc=2dln4n(ω)ρ(ω)(n(ω)+1)2(3)

其中, ω 是频率, k(ω )为消光系数、 ρ(ω )为幅值比函数, d为样品厚度, c为真空中的光速。

数据处理主要包括数据提取、 光谱分析和图像处理三部分, 数据提取是利用QT-TO1000太赫兹系统获取巴旦木的太赫兹光谱图像, 提取感兴趣区域的光谱信息。 光谱分析通过KS算法将光谱按照3∶1比例划分为建模集和预测集, 波段筛选出有效的波段点, 建立定性分析模型。 最后, 利用模型评价指标对所建立的模型进行评价。 利用matlab2018a对光谱数据进行处理和分析。 图像处理是对波段筛选后太赫兹特征图像和原始图像进行图像二值化并计算像素点, 证明太赫兹成像技术对检测巴旦木饱满度的可行性。 数据处理流程图如图3所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 数据处理流程图Fig.3 Data processing flowchart

太赫兹光谱数据量较大, 若要提高建模效率, 减少计算量, 则需要剔除光谱中的冗余信息, 提取与建模相关性高的特征波长信息。 本文采用的波段筛选算法为CARS、 无信息变量消除(uninformative variable elimination, UVE)、 连续投影算法(successive projections algorithm, SPA)、 蒙特卡罗无信息变量消除法(Monte carlo uninformation variable elimination, MCUVE)和遗传算法(genetic algorithm, GA)。

CARS^[19]是通过自适应重加权采样(ARS)技术选择出PLS模型中回归系数绝对值大的波长点, 去掉权重小的波长点, 利用交互验证选出RMSECV值最低的子集, 可有效寻出最优变量组合。

UVE^[20]是基于分析PLS回归系数b的算法, 核心是利用噪声的无关变量信息统计去选择光谱自身的特征变量, 得到最优变量组合。

SPA^[21]是一种前向循环选择方法, 利用向量的投影分析, 选取含有最少冗余度和最小共线性的有效波长。 连续投影算法对光谱数据进行初步压缩, 将优选出的波长按其对贡献值的大小进一步筛选, 输出按贡献值大小进行排列。

MCUVE^[22]是蒙特卡罗采样和无信息变量消除法的结合, 其变量重要性衡量依据是回归系数的稳定性值, 通过设定阈值, 去除稳定性小于阈值的波长。

GA^[23]是一种通过模拟进化过程来寻找最优解的算法。 遗传算法的基本原理是模拟自然选择、 遗传和突变等生物进化过程, 通过对种群的不断进化, 不断逼近最优解。

本研究对巴旦木样品所建立的模型通过分析建模集精度和预测集准确率对检测模型进行评价, 检测模型的建模集和预测集准确率越高, 则模型的精度越高。 准确率计算公式如式(4)所示

Accuracy=XY×100%(4)

式(4)中, X为正确分类的光谱条数, Y为总的光谱条数。

饱满度计算原理: 通过计算太赫兹图像壳仁像素点之比来计算巴旦木饱满度。 分为太赫兹图像检测饱满度P_detection(detected plumpness)和实际饱满度P_actual(actual plumpness):

检测饱满度原理如式(5)所示

Pdetection=AB×100%(5)

式(5)中, A为太赫兹图像中外壳所占的像素点个数, B为太赫兹图像中果仁所占的像素点个数。

实际饱满度原理如式(6)所示

Pactual=CD×100%(6)

式(6)中, C为实际照片中外壳所占的像素点个数, D为实际照片中果仁所占的像素点个数。

太赫兹光谱的信息十分丰富, 具有透过率、 吸收系数、 消光系数、 折射率、 介电常数、 相位角、 功率等太赫兹光学参数, 可以多维度的反映物质内部信息。 采用太赫兹成像技术采集巴旦木样品THz图像, 为了保证光谱数据的有效性, 分别提取在巴旦木满仁、 空壳和无样品感兴趣区域的时域光谱。 为了提高模型的鲁棒性, 选取的每类感兴趣点尽可能在该类区域图像上均匀分布。 图4(a)为巴旦木样品THz时域信号光谱图, 为了降低噪声干扰, 选择20~28 ps范围内的时域光谱进行分析。 从图4(a)中可以看出, 无样品区域、 巴旦木空壳和满仁区域的时域信号存在一定的时间延迟, 无样品区域的幅值最大, 空壳区域的幅值次之, 满仁区域的幅值最小, 振幅存在明显的衰减。 越厚的样品对太赫兹波的吸收越强, 这也是区分巴旦木满仁和空壳区域的有效依据。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 无样品区域、 空壳区域、 满仁区域的THz光谱 (a): 时域信号; (b): 频域信号Fig.4 Terahertz spectra of the sample free region, empty shell region and full almond region (a): Time domain signal; (b): Frequency domain signal

时域信号通过FFT(快速傅里叶变换)得到频域信号。 图4(b)为巴旦木样品频域信号光谱图, 在0.3~0.8 THz频率范围内, 幅值强度由高到低依次为: 无样品区域、 巴旦木空壳和满仁区域。 样品厚度越大, 其对太赫兹波的吸收能力越强, 也是从频域上验证了太赫兹技术检测巴旦木空壳区域和满仁区域的可行性。 图4(a)和图4(b)显示, 无论是从时域还是频域, 无样品区域、 巴旦木空壳和满仁区域太赫兹信号都有明显的差异, 文章中主要是通过太赫兹时域光谱中时域信号进行后续的特征波段提取, 这为太赫兹光谱结合模式识别方法对巴旦木的分级和分选提供了理论基础。

最小二乘支持向量机(least squares support vector machine, LS-SVM)^[24]是一种高效的有监督学习模型。 该方法适用于小样本、 高维度和非线性的数学问题, 具有较好的预测效果。 分别在太赫兹透射图像中无样品区域、 空壳区域和满仁区域三类感兴趣区域各选取400、 400、 320个像素点, 每个点提取一条太赫兹光谱。 采用KS算法^[25]将1120组太赫兹光谱数据按照3∶1的比例分成两组, 分别为建模数据和预测数据, 将无样品区域、 空壳区域和满仁区域光谱数据分别进行编号“1”、 “2”和“3”贴定标签。 为减少计算量, 剔除光谱中的冗余信息, 提取与建模相关性高的特征波长信息, 分别采用无波段筛选、 CARS、 UVE、 SPA、 MCUVE和GA分别建立LS-SVM模型, 表1为LS-SVM判别无样品区域、 空壳区域和满仁区域模型分类结果, 经过MCUVE波段筛选后建立的LS-SVM模型效果最好, 准确率达到99.05%, 但经过MCUVE筛选后波段点个数较多。 研究发现: 综合特征提取后的数据量和建模效果, 经过UVE波段筛选后的波段点个数为50, LS-SVM建模后预测集准确率为98.20%, 效果最好。

表1

Table 1

表1(Table 1)

表1 LS-SVM判别无样品区域、 空壳区域和满仁区域模型分类结果 Table 1 Classification results of sample free, empty shell and full almond region models by LS-SVM Modeling method Band screening method Number of band points RMSEP Precision of prediction set/% LS-SVM None 201 0.148 7 98.36 CARS 52 0.170 1 97.80 UVE 50 0.154 4 98.20 SPA 14 0.250 5 95.00 MCUVE 100 0.112 3 99.05 GA 36 0.234 1 95.73

表1 LS-SVM判别无样品区域、 空壳区域和满仁区域模型分类结果 Table 1 Classification results of sample free, empty shell and full almond region models by LS-SVM

为了更加直观地评价UVE-LS-SVM模型的预测效果, 通过混淆矩阵对预测集经过UVE特征提取后LS-SVM模型的预测结果进行分析。 如图5所示, 100条无样品区域谱线有2条误判为满仁区域谱线, 100条空壳区域谱线有1条被误判为无样品区域谱线、 1条被误判为满仁区域谱线, 80条满仁区域谱线有1条被误判为空壳区域谱线。 可以看出, 280条谱线中误判了5条光谱, UVE-LS-SVM模型能够有效地鉴别太赫兹图像中无样品区域、 空壳区域和满仁区域。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 UVE-LS-SVM模型预测集混淆矩阵图Fig.5 Confusion matrix of the prediction set of UVE-LS-SVM model Notes: SF=Sample free region; EM=Empty shell region; FU=Full almond region

由于K-最近邻(KNN)^[26]是一种没有确定训练阶段的分类方法, KNN算法的核心思想是如果样本在特征空间的k个最相邻的样本中大部分属于某一类, 那么该样本也属于这一类, 并且具有这类样本的特征。 与2.2处理光谱数据方法相同。 分别采用无波段筛选、 CARS、 UVE、 SPA、 MCUVE和GA分别建立KNN模型, 表2为建模结果, 综合对比下, 经过GA筛选后的波段点数为36个, 预测集准确率也为97.85%, 此时KNN模型效果最好。

表2

Table 2

表2(Table 2)

表2 KNN判别无样品区域、 空壳区域和满仁区域模型分类结果 Table 2 Classification results of sample free, empty shell and full almond region models by KNN Modeling method Band screening method Number of band points Number of prediction sets Number of correct Precision of prediction sets/% KNN None 201 280 276 98.57 CARS 52 276 98.57 UVE 50 275 97.85 SPA 14 262 93.55 MCUVE 100 276 98.57 GA 36 275 97.85

表2 KNN判别无样品区域、 空壳区域和满仁区域模型分类结果 Table 2 Classification results of sample free, empty shell and full almond region models by KNN

为了更加直观地评价GA-KNN模型的预测效果, 通过混淆矩阵对预测集经过GA特征提取后KNN模型的预测结果进行分析。 如图6所示, 100条无样品区域谱线没有误判, 100条空壳区域谱线有2条被误判为无样品区域谱线、 1条被误判为满仁区域谱线, 80条满仁区域谱线有2条被误判为空壳区域谱线。 可以看出, 280条谱线中误判了5条谱线, 说明GA-KNN模型能够有效地鉴别太赫兹图像中无样品区域、 空壳区域和满仁区域。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 GA-KNN模型预测集混淆矩阵图Fig.6 Confusion matrix of the prediction set of GA-KNN model Notes: SF=Sample free region; EM=Empty shell region; FU=Full almond region

随机森林(random forest)^[27]指的是利用多棵树对样本进行训练并预测的分类器, 是一种基于决策树的机器学习算法, 可以用于巴旦木太赫兹图像不同感兴趣区域的定性判别。 与2.2处理光谱数据方法相同。 分别采用无波段筛选CARS、 UVE、 SPA、 MCUVE和GA分别建立RF模型, 结果如表3所示: 综合对比下, 经过GA筛选后的波段点数为36个, 预测集准确率为98.21%, KS-GA-RF也是本文中效果最好的模型。

表3

Table 3

表3(Table 3)

表3 RF判别无样品区域、 空壳区域和满仁区域模型分类结果 Table 3 Classification results of sample free, empty shell and full almond region models by RF Modeling method Band screening method Number of band points Number of prediction sets Number of correct Precision of prediction sets/% RF None 201 280 274 97.86 CARS 52 275 98.21 UVE 50 275 98.21 SPA 14 262 93.57 MCUVE 100 275 98.21 GA 36 275 98.21

表3 RF判别无样品区域、 空壳区域和满仁区域模型分类结果 Table 3 Classification results of sample free, empty shell and full almond region models by RF

为了更加直观地评价GA-RF模型的预测效果, 通过混淆矩阵对预测集经过GA特征提取后RF模型的预测结果进行分析。 如图7所示, 100条无样品区域谱线没有误判, 100条空壳区域谱线有1条被误判为无样品区域谱线、 2条被误判为满仁区域谱线, 80条满仁区域谱线有2条被误判为空壳区域谱线。 可以看出, 280条谱线中误判了5条谱线, 说明GA-RF模型能够有效地鉴别太赫兹图像中无样品区域、 空壳区域和满仁区域。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 GA-RF模型预测集混淆矩阵图Fig.7 Confusion matrix of the prediction set of GA-RF model Notes: SF=Sample free region; EM=Empty shell region; FU=Full almond region

巴旦木图像太赫兹光谱全谱共2 250个点, 通过时域光谱区分出20~28 ps的波段点共201个点, 但是成像速度还是较慢, 采用波段筛选算法提取特征波段点, 根据特征波段点形成的图像合成波段筛选图, KS-GA-RF算法建模效果最好, 因此选取经过KS-GA处理后的36个波段点合成图像进行后续的巴旦木饱满度计算, 对应的波段点为27.97、 27.81、 27.61、 27.57、 20.77、 21.45、 22.65、 22.97、 23.01、 24.97、 26.53、 27.89、 20.85、 21.29、 22.25、 24.85、 25.29、 26.01、 26.81、 27.17、 27.21、 21.09、 22.20、 22.77、 22.89、 22.93、 23.13、 24.05、 25.33、 25.37、 25.89、 26.49、 26.57、 27.37、 27.77和27.85 ps。

巴旦木特征图像(图8)是根据GA波段筛选后, 经过36个波段筛选点所对应的特征波段点形成的单点图像叠加形成的, 未去壳的巴旦木经过太赫兹透射成像可以检测出巴旦木内部的饱满程度, 果仁的颜色为深蓝色, 果壳的颜色是浅蓝色, 区别不太明显, 经过特征提取后根据特征波段点合成的太赫兹特征图像可以看出外壳与果仁的轮廓, 通过巴旦木特征图像的壳仁像素点之比得出检测饱满度, 通过去一半壳巴旦木图片的壳仁像素点之比得出实际饱满度, 检测饱满度与实际饱满度相对比可以实现太赫兹成像技术检测巴旦木内部果仁饱满度的可视化检测。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 巴旦木特征提取对比图Fig.8 Comparative picture of feature extraction of almond

2.6.1 太赫兹图像计算检测饱满度

波段筛选后的太赫兹特征图像处理包括: 图像截取(伪彩色图转为JPG格式)、 JPG转RGB格式、 G通道分离提取、 图像二值化^[28]和计算像素点个数。 本节使用pycharm运行环境下的OpenCV库实现对图像的分析和处理。

外壳的处理是通过轮廓标记的方式将外壳部分与无样品区域区分开来; 仁的处理通过提取RGB图像中的G通道图像进行图像二值化(二值化阈值: 外壳=30; 果仁=40), 实现对巴旦木壳仁的分离。 总体流程如图9所示。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 太赫兹特征图像处理过程图Fig.9 Terahertz feature image processing flowchart

接着通过Img.sum函数计算出图像的像素总数值并将数值保存。 仁和壳的二值化图像通过相同的方式进行相加计算, 使用式(5)就可以得出太赫兹图像巴旦木的检测饱满度。 具体每个巴旦木的检测饱满度如表4所示。

表4

Table 4

表4(Table 4)

表4 太赫兹特征图像巴旦木检测饱满度 Table 4 Detection plumpness of almond in terahertz feature image Almond No. Total shell pixels Total kernel pixels Detection plumpness/% No.1 59 946 42 468 70.844 No.2 62 517 39 776 63.624 No.3 56 556 31 017 54.843 No.4 61 238 31 267 51.058 No.5 61 150 22 496 36.788 No.6 66 728 10 149 15.210

表4 太赫兹特征图像巴旦木检测饱满度 Table 4 Detection plumpness of almond in terahertz feature image

2.6.2 原始图像计算实际饱满度

针对原始图像的处理主要经过以下3个步骤: 轮廓标记、 图像二值化和计算像素点个数。 本节使用labelme软件进行人工图像分割, 使用Pycharm平台调用OpenCV函数库, 通过python编程语言实现原始图像的饱满度计算。

巴旦木原始图像的轮廓标记是通过labelme软件实现, 接着图像二值化(二值化阈值: 外壳=30, 果仁=40)和计算像素点个数, 总体流程如图10所示。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 原始图像处理过程图Fig.10 Original image processing flowchart

通过计算仁和壳的二值化图像的像素值使用式(6)得出原始图像巴旦木的实际饱满度。 每个巴旦木的实际饱满度如表5所示。

表5

Table 5

表5(Table 5)

表5 原始图像巴旦木实际饱满度 Table 5 Actual plumpness of almond in original image Almond No. Total shell pixels Total kernel pixels Actual plumpness/% No.1 63 302 48 360 76.396 No.2 67 749 44 731 66.025 No.3 54 594 29 953 54.865 No.4 60 339 23 628 39.159 No.5 54 965 13 664 24.859 No.6 No calculate 0 0

表5 原始图像巴旦木实际饱满度 Table 5 Actual plumpness of almond in original image

2.6.3 饱满度计算结果对比

通过表4检测饱满度和表5的实际饱满度信息, 为方便对比, 将检测饱满度和实际饱满度结果汇总在同一表格里。 如表6所示。

表6

Table 6

表6(Table 6)

表6 检测饱满度与实际饱满度结果对比 Table 6 Comparison of detection plumpness and actual plumpness results Almond No. Detection plumpness/% Actual plumpness/% Error /% No.1 70.844 76.396 5.552 No.2 63.624 66.025 2.401 No.3 54.843 54.865 0.022 No.4 51.058 39.159 11.899 No.5 36.788 24.859 11.929 No.6 15.210 0 15.210

表6 检测饱满度与实际饱满度结果对比 Table 6 Comparison of detection plumpness and actual plumpness results

通过表6可以看出, 6个巴旦木的检测饱满度与实际饱满度相差不超过16%, 由于巴旦木的不规则, 饱满度的误差还是在有效范围内的, 说明了太赫兹成像技术结合图像二值化可以有效计算巴旦木的饱满度。