所采用的泡桐实验样品全部来自古筝乐器厂, 进行相应的干燥和陈化处理, 含水率低于16%, 根据乐器技师按经验选出的高级品、 中级品以及普及品三个等级的木材各250块, 共750块板材用于近红外光谱谱线的筛选、 建模以及验证, 将样本标号为0— 749, 每块样品大小设置为10 cm× 10 cm× 1 cm。

采集光谱用的设备为PerkinElmer Spectrum 100 N系统(PerkinElmer Co, USA), 通过NIR积分球漫反射附件从每个样品的木块获得近红外光谱。 光谱分辨率为16 cm^-1, 收集10 000~4 000 cm^-1的光谱, 2 cm^-1为间隔, 即每个样本具有3 000点数据。 采集木材样品的两个面, 每个面上采集10个点, 共20个点最后取其平均值作为实验样本。 所使用的软件为Python3.6和TensorFlow系统, 根据Tensorboard可视化工具观察实验过程中参数的变化。

由于样品的颗粒程度不同, 光的散射等影响, 近红外光谱中往往包含一些与待测样品性质无关的信息, 因此对光谱数据进行预处理是非常必要的。 常用的预处理方法有Savitzky-Golay卷积平滑、 微分方法和标准正态变量变换(SNV)方法等。 其中微分方法可以降低或消除系统误差。 导数光谱提供比原始光谱更高的分辨率。 一阶导数可以清除不依赖于波长的基线漂移, 二阶导数可将线性漂移背景完全消除。 SNV方法用于校准基线漂移和光通行长度中存在的差异。 本文采用S-G微分方法和SNV。 方根误差反映了使用数据偏离真实值的程度, 信号平方和能够很好地反映出预处理手段提取信息的精密度, 所以光谱预处理手段选取依赖于均方根误差(RMSE)与信号平方和(SUMSQ)两个指标, 计算公式如式(1)和式(2), 其中S_i(n)为原始光谱信号, S˙i(n)为预处理后光谱信号, n为光谱数据长度。

RMSE=1n∑i=1n(Si(n)-S˙i(n))2(1)SUMSQ=∑i=1nS˙i2(n)(2)

将光谱数据进行可视化, 实验样本的原始光谱及经过预处理后的光谱信息如图1— 图4所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 原始光谱数据图Fig.1 Raw spectral image

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 一阶导数预处理图Fig.2 First derivative image

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 二阶导数预处理图Fig.3 Second derivative image

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 标准正态变量变换预处理图Fig.4 Standard normal variable transformation data image

近红外光谱定性分析中, 大量光谱数据的计算是非常复杂的, 通常进行数据压缩。 常用的数据压缩变换方法有傅里叶变换、 主成分分析和小波变换等, 本文采用主成分分析。

主成分分析是把多个特征映射为几个综合特征的一种统计分析方法, 计算出的主成分向量是各个原变量的线性组合, 主成分向量之间是正交的, 以此将数据降维, 排除化学信息中相互重叠的信息。 选择特征值较大的几个主成分作为实验数据参与分析。 主成分个数的选取很重要, 主成分个数过少不足以代表整个光谱信息, 过多则会引入噪声信息, 使模型的预测能力下降。 一般认为累计解释度达到85%以上即可代表整个光谱信息。 不同预处理方法基于实验数据进行主成分分析的累积解释度如图5所示。 可以看出3组光谱信息经过数据压缩后, 原始光谱经过4个主成分变量分析, 主成分累积解释度就已经接近99%, 一阶导数光谱取7个主成分时主成分累积解释度达到85%, 而二阶导数光谱主成分个数取56时主成分累积解释度才可以达到85%。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 主成分累积解释程度图Fig.5 Principal component cumulative interpretation image

光谱数据时常会出现异常样品的情况。 造成样品光谱数据异常的主要原因有测量环境的变化、 光谱扫描过程中的误操作和样品前处理不当等。 异常样品对于实验模型的影响是不能忽视的, 应该将其从建模集剔除。 异常样本的正确判别和剔除对于提高泡桐木材等级判别模型准确率具有十分重要的意义。

本实验所使用的剔除异常样本的方法为马氏距离方法, 它是一种有效的计算两个未知样本集的相似度的方法。 由于光谱数据计算量过大, 在进行剔除判别之前, 先对原始光谱数据进行主成分分析, 选取前两个主成分用于马氏距离计算。 马氏距离的计算公式如式(3), 其中x为待求样本点向量, μ 为实验样本均值向量, 协方差矩阵为S。

DM(x)=(x-μ)TS-1(x-μ)(3)

设置合适的距离阈值, 筛选实验建模集样本。 采用随机数种子机制和建模集与验证集的比例为4∶ 1原则随机选取实验样本, 构成建模集和验证集。 共剔除9个异常样本, 其中高级品6个, 中级品2个, 普及品1个, 异常样本进行显著放大标记, 具体的标记情况如图6所示, 其中横坐标为第一个主成分得分值, 纵坐标为第二主成分得分值。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 马氏距离异常建模集样本筛选图Fig.6 Mahalanobis distance method culling model set abnormal sample image

剔除异常样本后剩余的样本用于实验训练建模, 建模集与验证集具体分组情况如表1所示。

表1

Table 1

表1(Table 1)

表1 数据集划分 Table 1 Data set partition 建模集/例 验证集/例 剔除集/例 高级品 194 50 6 中级品 198 50 2 普及品 199 50 1 合计 591 150 9

表1 数据集划分 Table 1 Data set partition

Kmeans方法是典型的基于距离的分层次聚类方法, 在最小化误差函数的基础上将数据划分为预定的类数K, 采用距离作为相似性的评价指标, 即认为两个对象的距离越小, 其相似度就越大, 常用的距离度量方法包括欧几里得距离和余弦相似度等, 本文选取的度量方法为欧式距离。

建立一个三层的BP神经网络模型, 由于传统的BP网络学习率通常指定为一个固定的数值, 具有损失函数定义过于简单以及所使用的梯度下降法具有容易陷入局部最优解和训练速度慢等缺陷, 所以针对此劣势进行改进。

学习率对于模型有着显著的影响, 采用自适应学习率优化策略Adam优化方法进行实验, Adam算法集合了AdaGrad和RMSProp两个算法的优点, 能较好的处理噪声样本, 并且天然具有退火效果, 利用梯度的一阶矩估计和二阶矩估计动态调整每个参数的学习率, 对超参数的选择具有很强的鲁棒性。 交叉熵函数克服了方差代价函数更新权重过慢的问题。 本实验为多分类问题, 采用One-hot的编码方式设置样本类别, 最后一层的传递函数设置为Softmax函数。 相比于Sigmoid函数和tanh函数, Relu函数只需要一个阈值就可以得到激活值, 避免一系列复杂计算, 所以本模型输入层与隐藏层之间的传递函数选取Relu函数。

经过多次试验, 最终网络隐含层结点个数设置为15, 输出层节点数为3, 初始学习率设置为0.01, 训练迭代次数为100。

从图1可以看出, 原始光谱曲线谱峰重叠, 并且包含高频噪声信号; 在5 500~4 800和7 000 cm^-1附近不同质量的泡桐木材光谱波形具有微小的差异, 这使得区分差异品种成为可能, 同种等级质量的泡桐板材的样本光谱具有一定的聚集趋势, 但凭肉眼难以准确区分。

可以发现, 泡桐木材在波数10 000~7 100 cm^-1区吸收最少, 在6 806~5 192 cm^-1区吸收稍高, 在波数4 400~4 016 cm^-1区最高。 在8 248, 6 806, 5 804, 5 602, 5 192, 4 760, 4 400, 4 286和4 016 cm^-1处附近有明显的吸收峰; 在9 040, 8 280, 7 668, 6 084, 5 396, 4 978, 4 526, 4 354和4 166 cm^-1处有明显的肩峰。 纤维素、 半纤维素和木质素等为木材主要组成, 且所占比例为90%以上, 主要由O— H, C— H, CH₂和C— O等基团组成; 有机物在近红外谱区的吸收主要是含氢基团的倍频与合频的吸收; 5 192 cm^-1处附近由水中O— H键反对称的伸缩振动和变形振动引起, 6 806 cm^-1处附近由水中O— H键伸缩振动的第一泛音引起, 4 400 cm^-1附近由C— H键伸缩振动和变形振动引起, 8 248 cm^-1处附近为CH₂键伸缩振动的二级倍频, 8 280~7 668 cm^-1区间为C— H键伸缩振动第一泛音谱带。 根据特征峰以及吸光度大小可以发现, 与普及品相比高级品板材和中级品板材中C— H和O— H基团具有较低的吸收。

一般来说, 古筝面板木材要求含水率低于16%。 根据光谱吸收峰和肩峰发现, 水分子影响仍然较大, 具有较强的吸收。 泡桐木材水分子吸收带在5 396~4 978 cm^-1和6 800~7 000 cm^-1区域附近, 分别对应于水中O— H键伸缩振动、 变形振动的混合带和水中O— H键伸缩振动的第一泛音带。 为了避开水的吸收带, 选择10 000~7 000, 6 800~5 398, 4 976~4 000 cm^-1三个波段数据, 在此近红外光谱区域具有大量的有机物信息, 为探究最佳的实验波段, 将三个波段区间进行组合, 共7种实验光谱数据区间。

图7为使用原始光谱数据前两个主成分进行Kmeans方法聚类的分析图, K的值为3, 可见不同等级的泡桐板材均具有的一定的聚类趋势, 为进一步分类不同等级的泡桐木材等级提供了可能性, 但每一类中仍存在分错现象, 对于高级品和普及品分类效果稍好, 中级品分类效果最差, 要取得精准的定性分析结果, 还需要运用有监督的模式识别分析方法。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 聚类分析图Fig.7 Cluster analysis image

由图2和图3可知, 在一阶导数光谱中, 原光谱的峰值变成了零点, 而在原峰值两侧的拐点处分别出现了正的和负的峰值。 在二阶导数光谱中, 尽管符号出现了反转, 但原峰值的波长处仍然是二阶导数光谱的峰值, 因此导数光谱可以很方便地把握原光谱的变化趋势。 对于卷积微分预处理过程涉及到窗口的考量, Savitzky-Golay卷积滤波的效果随着选取窗宽不同而不同, 过小的窗口不足以提取全部的信息, 过大的窗口会过滤掉有用信息, 随着窗口大小的不同, 均方根误差(RMSE)与信号平方和(SUMSQ)的相应变化具体如图8和图9所示, 横坐标为窗口大小。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 二阶导数预处理评价指标图Fig.8 Second derivative preprocessing evaluation index image

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 一阶导数预处理评价指标图Fig.9 First derivative preprocessing evaluation index image

光谱数据去噪的衡量标准就是既能够过滤掉过多的噪声又可以保留尽量多的有效信息, 即去噪能力越强并且去噪后信号更接近于原始光谱信息时效果最佳, 所以均方根误差值与信号平方和数值应该尽量小。 综合对比均方根误差与信号平方和随滤波窗口数值大小变化的结果, 并且根据计算量和所提取的有用信息程度, 最终选取卷积平滑与一阶导数为本实验的预处理方式, 窗口大小为15。

本实验类别为高级品, 中级品和普及品三类, 通过One-hot编码, 其中[1 0 0]表示高级品, [0 1 0]表示中级品和[0 0 1]表示普及品, 以模型对验证集的识别率作为衡量指标, 不同主成分个数和不同组合实验光谱波段数据对于模型识别率的影响如表2所示, 可以看出随着主成分个数的增加, 模型识别率递增。 在主成分个数为11和实验光谱区间为10 000~7 000和4 976~4 000 cm^-1时验证集的准确率平均值达到99.7%, 并且参与实验计算的光谱数据最少和所选光谱区间涵盖了构成古筝板材化学成分的O— H和C— H等基团特征峰的全部特征信息。 为了减少篇幅和叙述重复, 将A表示为10 000~7 000 cm^-1波段, 将B表示为6 800~5 398 cm^-1波段, 将C表示为4 976~4 000 cm^-1波段。

表2

Table 2

表2(Table 2)

表2 模型分类结果 Table 2 Model classification result 主成 分数 波段 A B C A和B A和C B和C A, B和C 5 0.893 0.793 0.873 0.926 0.953 0.883 0.900 6 0.933 0.853 0.890 0.961 0.973 0.890 0.958 7 0.946 0.906 0.913 0.960 0.993 0.933 0.993 8 0.940 0.587 0.915 0.967 0.993 0.948 0.994 9 0.947 0.940 0.920 0.980 0.995 0.953 0.995 10 0.953 0.860 0.923 0.980 0.995 0.960 0.993 11 0.953 0.694 0.933 0.980 0.997 0.978 0.997

表2 模型分类结果 Table 2 Model classification result

在模型训练过程中, 使用了自适应学习率机制, 神经网络收敛速率较快, 在主成分个数为11和参与实验计算的光谱波段组合为10 000~7 000和4 976~4 000 cm^-1时, 使用Tensorboard可视化工具来观察参数的变化, 相应的准确率和损失值变化具体如图10和图11所示, 横坐标均为迭代次数。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 损失值变化图Fig.10 Loss value change image

图11

Fig.11

	Figure Option ViewDownloadNew Window
	图11 准确率变化图Fig.11 Accuracy rate change image