蒙古栎(Quercus mongolica)盛产于黑龙江省大、 小兴安岭一带, 这种木材密度高、 质地坚硬、 纹理直, 具有良好的抗弯曲性能, 且切削面光滑、 耐磨损, 是一种很好的家具材料。

试验用的蒙古栎样本采自黑龙江带岭林业局林场。 在蒙古栎人工林内, 取3株样木, 树龄20 a以上, 在每株标准木的胸高(1.3 m)附近连续截取50 cm圆盘。 实验室内对每个圆盘去皮, 按照木材物理力学试验方法(GB/T 1928— 2009)制取300 mm× 20 mm× 20 mm的抗弯力学试样并标记横切面。 挑选出切面标准、 表面平整、 无缺陷试样126个, 编号后, 将试样放入恒温恒湿箱内调节含水率至12%, 进行近红外光谱扫描和力学性质测试。

按照《木材抗弯强度试验方法》(GB/T 1936— 2009)测定蒙古栎无疵试样抗弯强度, 将采集的126个样本数据按照2:1的比例随机分为校正集和验证集, 其中校正集共84个样品, 验证集共42个样品。 校正集用来建立校正模型, 预测集对模型进行验证。

光谱仪采用INSION公司近红外光纤光谱仪对木材样本进行光谱扫描。 光谱仪采集波长范围在900~1 700 nm, 共117个波段, 数据采样间隔小于16 nm光谱分辨率。 实验环境温度20 ℃、 平均相对湿度50%、 光线均匀。 使用两分叉光纤探头采集试样表面的近红外光谱, 光源使用钨光源。

光纤探头垂直固定在金属支架上, 试样放置于支架底面, 非接触测量距离1 mm, 光斑直径5 mm。 软件设置为积分时间为20 ms, 30次均值滤波, 用聚四氟乙烯白板进行标定后, 对横切面进行光谱采集, 并记录保存。

由于光谱采集过程中的抖动、 光线散射等因素, 会产生基线漂移、 噪声干扰。 为提高模型性能, 需要对近红外光谱数据进行预数据处理后再进行建模分析。

采用一阶导数与S-G卷积平滑结合算法, 对数据进行预处理。 首先对光谱数据进行一阶导数处理, 以消除基线和其他背景干扰。 再用S-G多项式卷积对一阶导数光谱数据进行平滑去噪处理, 确保信号形状、 宽度不变, 保证数据不失真。

式(1)为对相邻光谱点相减的一阶导数处理

dxdλ=xi+1-xiΔλ(1)

式中: x为吸光度; λ 为波长。

S-G平滑把光谱区间的连续点作为一个窗口, 用多项式对窗口内的光谱数据做最小二乘拟合, 得到相应的多项式系数来计算该窗口中心波长点的导数平滑值。

S-G平滑多项式可以定义为

f(x)=∑k=0nbnkxk(2)

利用光谱数据拟合多项式系数b_nk, 计算窗口中心点k阶导数光谱值a_nk, 如式(3)

dkfidiki=0=k!bnk=ank(3)

k为导数的阶数; n为多项式最高次数; b_nk为多项式系数; a_nk为窗口中心点S-G平滑光谱值。

卡尔曼滤波通过状态预测结合观测来估计当前系统的状态, 通过上一时刻状态最优值来预测当前状态, 通过状态预测值和状态观测值来计算状态最优估计值。

在光谱特征选择中, 预测模型选用PLS, 在此, 以PLS模型系数矩阵作为系统状态, 抗弯强度作为观测值, 即可定义线性系统和量测方程, 系统在迭代中不断更新状态方程以预测系统状态, 即PLS模型系数。

首先定义光谱-抗弯强度模型系统方程:

状态方程: b_k+1=b_k+w_k(4)

观测方程: y_k=x_kb_k+v_k(5)

式(4)和式(5)中, x为样本的光谱数据; b为系统模型系数矩阵; y为样本抗弯强度; 随机信号w_k和v_k表示过程激励噪声和观测噪声, 且p(w)~N(0, Q), p(v)~N(0, R)。 代入离散卡尔曼滤波器方程

b˙k-=b˙k-1, Pk-=Pk-1+Q, Kk=Pk-XT(XPk-XT+R)-1, b˙k=b˙k-+Kk(yk-Xb˙k-), Pk=(I-KkX)Pk-（6）

其中 b˙k-表示先验状态估计, b˙k为后验状态估计, X为迭代光谱数据, Q为过程激励噪声协方差矩阵, R为观测噪声协方差矩阵, K_k为卡尔曼增益矩阵, Pk-为先验估计误差的协方差, P_k后验估计误差的协方差。

卡尔曼滤波器不断更新模型系数矩阵b和估计误差协方差矩阵P, 迭代完成后, 使用模型系数矩阵与后验估计误差标准差的比值R_kf表征波长对构建模型的重要性, R_kf定义如式(7)所示。

Rkf=b/P(7)

波长点相关程度越高时, 模型系数矩阵b对应权值越高, 反之越低; 后验估计误差标准差越小表示模型系数矩阵b对应权值越准确, 反之越高。 依据R_kf的大小对光谱波长进行排序, 选择R_kf较高光谱维度进行建模。

卡尔曼滤波的特征选择算法流程如图1所示。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 基于卡尔曼滤波的波长选择流程图Fig.1 Flow char of KF based feature selection

对126个蒙古栎试样进行光谱采集后, 进行光谱数据预处理, 使用一阶导数与S-G平滑预处理方法。 一阶导数采用差分法, S-G平滑的多项式阶数为5、 窗口宽度为9时, 吸收峰集中, 光谱平滑, 效果较好。 原光谱、 预处理结果分别如图2、 图3所示, 一阶导数与S-G平滑可以实现在滤除噪声的同时突出信号的主要特征。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 原始光谱Fig.2 Original spectra

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 一阶导数和S-G平滑预处理光谱Fig.3 First order derivative and S-G smoothing preprocessing

近红外区谱带复杂、 重叠多、 数据量大, 其中含有大量与分类无关的冗余信息, 降低了分类器的分类性能。 变异系数为标准差与均值的比值, 可以通过计算光谱预处理后的变异系数观察光谱数据的冗余程度。 图4为光谱波长的变异系数, 可以看出光谱波段中存在大量冗余信息, 近红外光谱在1 400~1 500 nm波段、 1 600~1 700 nm变异系数值较高, 说明光谱在此波段有效信息多, 冗余信息较少。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 各波长变异系数Fig.4 The coefficient of variation for different wavelengths

以经预处理后116个波长点作为选择对象。 以校正集82个样本作为迭代样本, 状态方程初值均设为0, 卡尔曼过程激励噪声方差Q设为0.05, 观测噪声方差R设为0.01。 经校正集样本84次迭代后, 以R_kf进行评价。 图5为不同波长点对应的R_kf值。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 卡尔曼特征选择结果Fig.5 Based KF feature selection results

按照大小对R_kf进行排序, 得到光谱影响程度的排序表, 表1为光谱影响程度排序的前15组数据。

对光谱进行KF处理后, 运用PLS对不同特征数量的近红外光谱波长进行偏最小二乘回归, 并选择预测误差均方根(RMSEP)作为参照对比。 图6为卡尔曼特征选择后进行的偏最小二乘回归, 由图可得, 随着建模的变量数增加, RMSEP在不断减小; 当选择变量数为18时, RMSEP最小; 之后, 随着变量数的增加, RMSEP逐渐增加。 故选择建模变量数为18, 其预测相关系数r=0.81, RMSEP=6.59。 图7中圆点为经特征选择的特征波长点位置。

表1

Table 1

表1(Table 1)

表1 光谱影响程度排序 Table 1 List of spectral effects based on Rkf 序号 Rkf 光谱波长/nm 1 495.87 1 408.34 2 474.71 1 416.59 3 440.41 1 697.90 4 394.03 1 689.61 5 388.85 1 424.85 6 378.75 1 400.08 7 331.25 1 706.19 8 269.55 1 391.83 9 262.75 1 433.1 10 253.95 1 681.32 11 172.85 1 822.41 12 172.65 1 383.58 13 151.98 1 673.03 14 146.67 1 789.18 15 143.82 1 747.67

表1 光谱影响程度排序 Table 1 List of spectral effects based on Rkf

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 不同变量数的卡尔曼-偏最小二乘模型预测误差Fig.6 Predictive errors with different variable number based on KF-PLS

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 特征波长选择变量位置Fig.7 The variable positions of the selected wavelength poiots

选择无信息变量消除法(UVE)和连续投影算法(SPA)与KF优选特征的结果进行比较。

UVE算法通过对PLS模型中添加一组与原始变量数量相同的白噪声变量, 然后基于PLS模型的交叉留一法得到每个变量的回归系数, 用模型系数的稳定值剔除无关变量。 图8为UVE法筛选波长示意图, 纵轴实线左侧为波长变量稳定性分布曲线, 右侧为噪声变量稳定性分布曲线, 2条虚线外的波长变量部分被认为是有效信息, 即在2条虚线外的波长变量被保留。 经UVE选择后, 得到42个特征波长点, 进行PLS建模验证, RMSEP=6.92, 相关系数r=0.76。

SPA(连续投影算法)是从光谱变量中找出冗余信息最少的变量组, 该变量组内的变量之间的共线性达到最小。 图9为SPA方法进行光谱选择后, 不同变量数的PLS模型RMSEP曲线, 当变量数达到14时, RMSEP趋于不变, 故选择建模变量数为14, 其RMSEP=14.25, 相关系数r=0.65。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 UVE筛选波长示意图Fig.8 UVE Screening wavelength diagram

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 基于SPA的不同变量数PLS模型预测误差Fig.9 Predictive errors of PLS with different variable numbers based on SPA

图10为三种方法选择特征波长对比图, UVE选择波长数量较多, 相比之下SPA和KF选择特征波长更为精准, 且重叠部分较多。 表2计算了不同特征选择方法优选出的特征变量数, 及运用特征优选后PLS建模的相关评价指标。 使用卡尔曼光谱特征选择后, 预测相关系数相比UVE、 SPA分别提高了0.05和0.16, RMSEP降低了2.33和7.66。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 不同特征选择方法的波长对比Fig.10 The comparison of the selected wavelength points based on different methods

表2

Table 2

表2(Table 2)

表2 基于不同优选特征的PLS建模结果对比 Table 2 The comparison of different features based PLS 建模方法 变量数 RMSEP 预测相关系数r UVE-PLS 42 8.92 0.76 SPA-PLS 14 14.25 0.65 KF-PLS 18 6.59 0.81

表2 基于不同优选特征的PLS建模结果对比 Table 2 The comparison of different features based PLS

为了验证KF特征选择后, 模型预测的快速性, 在2.7 GHz Intel Core i5硬件环境下, 使用Matlab 2016b软件运行程序, 调用tic和toc函数分别记录了42个测试样本的全波

长光谱提取、 光谱预处理、 特征光谱提取与识别的时间, 各步骤的平均时间如表3所示。 通过对各步骤平均时间的累计可以计算出基于KF-PLS的蒙古栎抗弯强度预测的平均时间为65.8 ms。

表3

Table 3

表3(Table 3)

表3 基于KF-PLS的各步骤预测时间统计 Table 3 Predicted time statistics of the various steps through KF-PLS 光谱处理步骤 平均时间/ms 光谱提取 20.0 光谱预处理 5.1 KF-PLS模型预测 40.7

表3 基于KF-PLS的各步骤预测时间统计 Table 3 Predicted time statistics of the various steps through KF-PLS