桦木(Betula spp.)木材为淡褐色至红褐色, 纹理直且明显, 材质结构细腻, 质地较软; 其力学强度性能优秀, 富有弹性, 吸湿性大, 加工性能好, 切面光滑。 桦木样本来自黑龙江带岭林业局林场, 共取3株, 树龄均在20年以上, 在样木胸高(1.3 m)附近截取约50 cm的圆盘后去皮, 气干后锯解成力学试材毛坯条。 参照GB/T 1938— 2009《木材顺纹抗拉强度试验方法》, 制取厚度为4 mm, 宽度为15 mm尺寸的顺纹抗拉力学试样, 挑选无疵试样150条并编号, 将试样放入恒温恒湿箱内调节至含水率为12%, 保持室温为(20± 2)℃、 相对湿度为65%± 3%, 然后进行近红外光谱采集和力学性能测试。

采用INSION公司近红外光纤光谱仪进行光谱采集, 光谱仪波长范围900~1 700 nm, 在环境温度20 ℃, 平均相对湿度50%条件下, 在每个试件的径、 弦切面上分别均匀扫描采集4个样本点, 每点扫描30次自动平均为1个光谱, 并记录保存, 试件原始光谱如图1所示。 由于原始光谱存在基线漂移、 样本颗粒大小不一和光散射等噪声的影响, 采用多元散射校正(MSC)、 一阶求导和卷积平滑(SG)相结合的方法进行预处理。 经MSC-1ST-SG预处理后的光谱曲线如图2所示, 由图可知, 经预处理后的光谱图像更加平滑且各吸收峰更突出。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 原始光谱Fig.1 Original spectrum

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 MSC-1ST-SG预处理Fig.2 MSC-1ST-SG pretreatment

根据《木材顺纹抗拉强度试验方法》(GB/T1935— 2009), 测定桦木的无疵试样抗拉强度真值。 对150组桦木样本, 按8:2的比例划分校正集与预测集。 采用光谱-理化值共生距离(SPXY)算法对其进行分割, 保证数据集间的一致性, 得到校正集样本120个, 预测集样本30个。 由表1可知, 抗拉强度范围在49.97~220.66 MPa内, 且各数据集的平均值与标准差相差不大, 预测集信息被校正集信息所覆盖。

表1

Table 1

表1(Table 1)

表1 样本校正集和预测集抗拉强度的测定结果 Table 1 Test results of compressive strength of sample correction set and prediction set 样本 数量 抗拉强度/MPa 最大值 最小值 平均值 标准差 训练集 120 220.57 49.97 130.10 39.43 预测集 30 204.66 51.50 133.48 39.12 数据集 150 220.57 49.97 131.12 39.17

表1 样本校正集和预测集抗拉强度的测定结果 Table 1 Test results of compressive strength of sample correction set and prediction set

蒙特卡罗无信息变量消除法(MC-UVE): 通过添加与原变量数相同的白噪声变量, 再采用交叉留一法得到每个变量的回归系数, 采用模型的稳定性值剔除无信息变量。 通过蒙特卡洛采样可以从校正集中随机选取一定比例的样本建立PLS回归子模型, 比例设置为P_MCUVE, 得到回归系数b, 重复n_MCUVE次, 得到回归系数矩阵B(n× m), 则样本每个变量稳定性值S_i根据式(1)可得

$\begin{array}{c} S_{i}=\frac{\operatorname{mean}\left(\boldsymbol{B}_{i}\right)}{\operatorname{std}\left(\boldsymbol{B}_{i}\right)} \end{array}$(1)

式(1)中, mean(B_i)和std(B_i)为回归系数矩阵B的平均值和标准差。 稳定性值S_i越大, 则表示该特征波长越重要。

迭代变量子集优化算法(IVSO): 通过加权二进制采样(WBMS)^[12]建立子模型, 在每一轮迭代中, 对子模型的回归系数归一化处理后求和, 以评价各特征波长重要性, 逐步优选出最佳波长子集。 其具体过程如下:

(1)利用WBMS创建维数为l× p波长变量矩阵, 并记录WBMS抽样后的变量数L₁。

(2)建立PLS子模型计算得到回归系数矩阵W, 如式(2)

$\begin{array}{c} \boldsymbol{W}=\left[w_{1}, w_{2}, w_{3}, \cdots, w_{l}\right]^{\mathrm{T}} \end{array}$(2)

式(2)中, w_j(1≤ j≤ l)表示第j个子模型的回归系数绝对值。 将矩阵W进行如式(3)归一化矩阵C处理。

$\begin{array}{c} \boldsymbol{C}_{i j}=w_{i j} / \max \left(w_{j}\right) \end{array}$(3)

式(3)中, w_ij表示w_j中第i个变量的回归系数绝对值, max(w_j)表示行向量w_j的最大值。

(3)将矩阵C的每一列求和, 记为向量s=[s₁, s₂, …, s_p]。 将s中的元素从小到大排序, 值越大表示对回归方程贡献越大。

(4)根据变量数L₁和对回归方程贡献的大小顺序排名构建L₁个子模型, 将RMSECV值最小的子模型中的变量子集作为该迭代轮的目标变量子集, 记录该值和对应模型的变量子集L₂。

在每一轮迭代中, 定义第i个变量的权值为

$\begin{array}{c} \omega_{i}=s_{i} / \max (s), i=1,2, \cdots, P \end{array}$(4)

式(4)中, max(s)表示向量s中的最大值, ω _i表示第i个变量的采样权重。

重复(1)— (4), 直到L₁=L₂, 在得到多个变量子集中, 选择RMSECV值最小的变量子集作为最优变量子集。

变量组合集群分析法(VCPA): 采用指数衰减函数(EDF)确定变量的剩余数量。 当剩余变量的下降速度和剩余变量的值成比例, 剩余变量服从指数衰减, 指数衰减微分方程为式(5)

$\begin{array}{c} \frac{\mathrm{d} \alpha(i)}{\mathrm{d} i}=-\theta \alpha(i) \end{array}$(5)

指数衰减微分方程的解如式(6)所示

$\begin{array}{c} r_{i}=\alpha \mathrm{e}^{-\theta i} \end{array}$(6)

式(6)中, i为迭代次数。 α 为迭代次数i为0时的初值。 θ 是常数, 其值由运行次数i与r_i决定。

在每次EDF运行中, 通过二进制矩阵采样BMS^[13]使变量与变量之间随机组合, 生成一个不同变量组合的子集, 建立子模型并记录RMSECV。 通过从特征变量子模型集群中提取RMSECV值较小的10%最优子模型, 计算出子模型中所有特征变量的次数频率。 频率越高, 特征变量越重要, 以此获得最佳的特征变量组合。

MC-UVE通过蒙特卡洛采样的方式建立回归子模型, 根据稳定性值判定特征的信息量, 会带来一定的随机性, 导致一定量的非信息变量存在。 特别是针对光谱特征相互差异较小的数据时会更难分辨, 从而影响最后的建模效果。 有必要对采用MC-UVE优选后的特征波长作进一步优化。 而IVSO算法与VCPA算法在面对特征变量较多的光谱数据集时, 需要更多的迭代过程, 导致计算效率低。 MC-UVE与这两种算法的有效结合可以解决上述问题, 对全光谱的特征波长根据稳定性值作粗略的选取, 再根据RMSECV值最小的变量子集作为最优变量子集, 特征波长选择组合方法将更准确剔除非信息变量, 提高计算效率, 提升回归模型的效果。

偏最小二乘回归模型(PLS)是常用的近红外光谱预测模型, 其主要将高维数据投影到由几个主元构成的低维空间, 通过主元来保留原数据中的重要信息, 每个主元间相互正交, 从而消除原数据中的冗余信息。 回归模型评价指标采用决定系数(R²), 预测均方根误差(RMSEP)和性能偏差比(RPD)。 R²被定义为模型输出向量与实际输出向量的协方差除以其标准差的乘积, RMSEP定义为模型输出量与实际之差期望值的平方根, RPD则为实际测量值的标准偏差(SD)与预测模型的RMSEP之间的比率。 R²越接近于1、 RMSEP的值越小、 RPD的值越大则回归模型预测效果越好。

将预处理后的光谱数据分别采用VCPA、 MC-UVE、 IVSO、 MC-UVE-VCPA和MC-UVE-IVSO算法进行波长筛选。 实验中设置PLS的最大主成分数为10, 交叉验证为5折, 每个特征波段选择实验各进行30次, 选取这30组中RMSECV最小的模型进行分析比较, 如图3所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同算法的波长选择过程 (a): VCPA; (b): MC-UVE-VCPA; (c): IVSO; (d): MC-UVE-IVSOFig.3 Wavelength selection process of different algorithms (a): VCPA; (b): MC-UVE-VCPA; (c): IVSO; (d): MC-UVE-IVSO

VCPA的二进制采样(BMS)采样数n_bms设为2 000, 指数递减函数(EDF)运行数i为50, 初值α 为512, 比率r_i为0.1, 如图3(a)所示。 VCPA在EDF运行到35次处, 最小交叉验证均方根误差(RMSECV)为12.09, 此时得到12个波长点, 占总波长的3%。

MC-UVE的蒙特卡洛采样次数n_mcuve设为4 000, 校正集采样比率P_mcuve为0.8, RMSECV最小为9.37, 此时模型稳定性最高, 得到波长点245个, 占总波长的48%。

MC-UVE-VCPA的参数与VCPA的参数相同, 如图3(b)(MC-UVE-VCPA), VCPA在迭代运行到34次处, 为10.93, 此时得到10个波长点, 占总波长的3%。

IVSO加权二进制采样(WBMS)采样数n_wbms设为4 000, 如图3(c)所示, 在第129轮迭代时, RMSECV为9.21最小, 此时得到67个波长点, 占总波长的13%。

MC-UVE-IVSO的参数与MC-UVE和IVSO相同, 如图3(d)所示, 在第52轮迭代时, RMSECV为10.55最小, 此时得到98个波段, 占总波长的19%。

近红外光谱主要是有机分子活原子振动在2 000 cm^-1以上的倍频与合频吸收, 对应的C— H, O— H, N— H等含氢基团的倍频与合频吸收带。 木材的主要成分是纤维素和木质素, 对应的二级倍频主要在1 000~1 400 nm处, 一级倍频主要在1 400~1 800 nm处^[14], 在近红外光谱区域有丰富的吸收信息。 图4为不同算法对波长的选择结果, 可以看出在980、 1 260、 1 460以及1 660 nm附近有主要吸收峰, 同时在1 360和1 570 nm处有小的吸收峰。 从各算法选择的波长点可以看出: MC-UVE选取的波长多, 该处近红外光谱数据间的差异较小时, 会保留冗余信息, 降低建模准确性。 IVSO和MC-UVE-IVSO所选特征波长位置基本相同, 同时对于1 660 nm附近存在大量噪声, 只选取了少量的特征波长, 但MC-UVE-IVSO在主要和次要的吸收峰处选取的波长数明显多于IVSO, 将提高模型的精确性, VCPA与MC-UVE-VCPA同样遵循这条规律, 但VCPA所选特征波长过少会影响模型精确性。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同算法的波长选择对比Fig.4 Comparison of wavelength selection of different algorithms

在优选特征波长基础上, 构建木材抗拉强度PLS预测模型。 基于不同的特征波长建模预测精度如表2所示。 MC-UVE-IVSO在校正集中决定系数(R²)值为0.91, 校正均方根误差(RMSECV)值为10.55。 在预测集中, 决定系数(R²)值为0.94, 预测均方根误差(RMSEP)值为7.50, 性能偏差比(RPD)值为3.17。 相比于VCPA、 MC-UVE、 MC-UVE-VCPA、 IVSO, MC-UVE-IVSO在校正集和预测集中的PLS模型效果最优。 由于木材的各向异性及不均匀性等特点, 增加了抗拉强度的预测难度, 据PLS建模的结果MC-UVE-IVSO算法相比于其他4种算法所提取出的特征波长信息与抗拉强度真值拟合度更高, 模型剔除了非信息波长, 提升了模型预测准确性。

表2

Table 2

表2(Table 2)

表2 不同算法的PLS模型对比 Table 2 PLS model comparison of different algorithms Algorithm Number of wavelengths Calibration set Prediction set R2 RMSECV R2 RMSEP RPD FULL 512 0.94 15.62 0.92 17.91 2.81 MC-UVE 245 0.96 9.37 0.93 11.70 2.91 VCPA 12 0.88 12.09 0.82 14.91 2.25 MC-UVE-VCPA 11 0.93 10.93 0.87 12.12 2.28 IVSO 67 0.98 9.21 0.93 8.47 2.78 MC-UVE-IVSO 65 0.91 10.55 0.94 7.50 3.16

表2 不同算法的PLS模型对比 Table 2 PLS model comparison of different algorithms

箱形图是用于展示数据分散情况的统计图, 由于上述特征波段选取算法具有随机性, 30次实验得到的优选波长结果存在区别, 所建立的PLS模型的RMSEP各不相同。 统计不同模型的RMSEP, 并绘制如图5所示的预测箱形图。 可以看出VCPA-PLS相对预测精度RPD最差, 且模型输出与数据集实际标签间的误差RMSEP最大。 30次训练得到MC-UVE-PLS模型的RMSEP两端极值距离过大, 说明该算法在选取波长建模的稳定性不高, 预测效果最差。 而第一、 四分位数(First Quartile, Q1), 第三、 四分位数(Third Quartile, Q3)和RMSEP均优于使用MC-UVE-VPCA和VPCA优化的PLS模型, 进一步证明MC-UVE-IVSO算法在处理多变量波长的建模稳定性。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 不同算法的稳定性比较Fig.5 Stability comparison of different algorithms

在Intel Core i7-9750H@2.6GHz的硬件环境下, 使用Matlab 2018b运行程序, 调用etime函数计算各算法特征波长选择时间, 并计算平均值。 如表3所示, MC-UVE-IVSO算法在经过MC-UVE的前期选择后, 降低了IVSO变量剔除的迭代过程, 相比于IVSO算法, 减少43%的计算时间, 优化了波长选取过程。

表3

Table 3

表3(Table 3)

表3 不同算法的波长选取时间对比 Table 3 Comparison of wavelength selection time of different algorithms Algorithm Average time/s VCPA 78.09 MC-UVE 25.48 MC-UVE-VCPA 65.91 IVSO 164.05 MC-UVE-IVSO 92.35

表3 不同算法的波长选取时间对比 Table 3 Comparison of wavelength selection time of different algorithms