采用5种东北常见的杨木、 桦木、 樟子松、 白松、 落叶松和15种其他树种分别为巴西花梨、 红花梨、 印尼菠萝格、 美国红橡、 非洲卡斯拉、 缅甸柚木、 海棠木、 克隆木、 金丝柚木、 漆木、 水煮柚木、 南美柚木、 唐木、 香樟木、 白蜡木作为实验样本。 将其加工成木块。 每个树种采集80幅高光谱图像, 则共有1 600幅高光谱图像作为实验样本。

图1为显微高光谱成像系统示意图, 该系统主要是由便携式高光谱成像仪(SOC710VP, USA)和CCD相机两部分构成, 成像波段为可见光/近红外波段, 光谱分辨率为5.2 nm, 光谱范围372.53~1 038.57 nm。 采集木材横切面的显微高光谱影像, 所采集的高光谱图像尺寸是520× 696× 128。 台式机电脑配置为i7处理器, 4G运行内存。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 木材体视显微高光谱成像系统 (a): 高光谱成像系统; (b): 采集的木材显微高光谱图像立方体(木材树种是印尼菠萝格)Fig.1 A wood stereomicroscopic hyper-spectral imaging system (a): The hyper-spectral imaging system; (b): The wood hyper-spectral image cube for wood specties “ Merbau”

这里强调一点, 考虑到仪器的价格, 我们使用了价格较低的便携式高光谱成像仪, 它的成像波段在可见光和少部分近红外波段内。 因此, 实验时应该保持外界环境例如温度、 湿度、 样本保存时间等因素尽量恒定, 尽量确保各个木材树种样本的颜色保持不变。 否则, 可能造成可见光波段内的光谱图像变化, 影响后续的特征提取和分类识别精度。

由于CCD相机暗电流及不均匀照明强度会影响高光谱图像质量, 应校准原始图像。 校正公式如式(1)

R=IO-IBIW-IB(1)

式(1)中, I_O为原始未修正的高光谱图像, I_B为扫描得到的黑板校正图像, I_W为白板校正图像, 从而得到校正后的图像R。

1.3.1 ROI的选择

采用ENVI5.1软件中region of interest工具进行高光谱图像感兴趣区域(region of interest, ROI)的提取, 以此增加实验样本数量。 本实验在样本中心区域手动选取像素的正方形区域作为ROI, 每个图像选取3个ROI, 共提取4 800个ROI, 然后计算ROI内所有像素点光谱反射率的平均值。

1.3.2 光谱预处理

为了提取高光谱图像中光谱有效信息, 提高建模稳健性和预测精度, 需要对光谱数据进行预处理来减少光谱信号中的光散射、 高频噪声等干扰信息。 常用的预处理方法有中值滤波平滑、 变量标准化校正(standard normal variate, SNV)、 多元散射校正(multiplicative scatter correction, MSC)、 Savitsky-Golay平滑算法(Savitzy-Golay smoothing, SG)、 一阶导数(first derivative, FD)等。 本工作采用SNV, MSC和SG三种预处理方法, 根据模型的预测精度选择最佳的预处理方法进行建模。

1.3.3 光谱特征选择

所获得的高光谱图像样本具有波段多、 数据量大等特点, 存在大量冗余和多重共线性问题, 在采集波段光谱中很多波长与分类预测无关, 因此删除这些不相关波长可以提高运行速度。 采用连续投影算法(successive projection algorithm, SPA)和竞争性自适应重加权算法(competitive adaptive reweighted sampling, CARS)进行特征波长选择。

1.3.4 纹理特征提取

高光谱图像具有图谱合一的特点, 除了对高光谱图像的光谱维进行特征选择外, 还需要对其空间维进行纹理计算。 这里采用经典的灰度共生矩阵(gray level co-occurrence matrix, GLCM), GLCM基于像素灰度的空间相关性表示纹理特征。 首先对高光谱图像样本进行主成分变换, 最后选择第一主成分图像, 距离参数值设为1。 对第一主成分图像基于GLCM分别提取图像中的0° , 45° , 90° 和135° 方向上的能量、 熵、 惯性矩、 相关性16个特征值作为纹理特征。

1.4.1 支持向量机

支持向量机(support vector machine, SVM)是由Vapnik^[15]等提出, 是一种以结构风险最小化原则为基础的模式识别算法, 克服了在传统机器学习中维数灾难问题, 实验表明它在小样本数据集中的分类问题上有显著优势。 SVM求解最优分类超平面问题等价求解如下方程式(C, ε _i是相应参数, ϕ 是非线性映射函数, y_i∈ [-1, 1])

minw, b, εi12‖w‖2+C∑iεis.t. yi(< ϕ(xi), w> +b)≥1-εi, εi≥0, ∀i=1, …, n(2)

引入Lagrange乘子α _i可以推导得出二次规划问题

maxαi∑iαi-12∑i, jαiαjyiyjK(xi, xj)s.t.∑iαiyi=0, 0≤αi≤C, i=1, 2, …, n(3)

其中K(· )是满足Mercer条件的核函数, 相应的SVM判别函数为

f(x)=sgn∑i=1nyiαiK(xi, x)+b(4)

1.4.2 复合核函数

为了提升SVM分类器的性能, 选择合适的核函数至关重要。 因此, Gustavo等提出了基于复合核函数的新型SVM分类器^[16], 将其应用于高光谱遥感图像的地物分类中。 在遥感图像的地物分类中, 处理的基本数据单元是像素, Gustavo等重新定义一个新的像素点特征矢量x_i={ xis, xiw}, 它是由光谱维特征 xiw∈ RNw和空间维特征 xis∈ RNs组成, 其中N_w代表光谱维特征个数, N_s代表空间维特征个数。 相应的, K_s代表空间核函数, K_w代表光谱核函数, 用它们或者复合核函数替换式(3)和式(4)中原始的核函数K, 具体的4种复合核函数表示如下^[16]。

(1) Stacked Features Approach。 这里重新定义一个像素点特征向量x_i={ xis, xiw}, 它是由该像素点空间维特征和光谱维特征直接连结, 相应复合核函数表达式为

K{s, w}≡K(xi, xj)=< ϕ(xi), ϕ(xj)>(5)

(2) Direct Summation Kernel。 在希尔伯特空间H中定义φ ₁(· )和φ ₂(· )两个非线性变换, 再生成新的变换组合ϕ (x_i)={φ ₁( xis), φ ₂( xiw)}, 进而构造复合核函数, 这里dim(K)=dim(K_s)=dim(K_w)=n× n。

K(xi, xj)=< ϕ(xi), ϕ(xj)> =< {φ1(xis), φ2(xiw)}, {φ1(xjs), φ2(xjw)}> =Ks(xis, xjs)+Kw(xiw, xjw)(6)

(3) Weighted Summation Kernel。 这种情况可以表示为B的加权形式, 其中0≤ μ ≤ 1, μ 为权因子

K(xi, xj)=μKs(xis, xjs)+(1-μ)Kw(xiw, xjw)（7）

(4) Cross-Information Kernel。 它考虑光谱维信息和空间维信息的交叉关系, 在希尔伯特空间H中定义一个非线性映射φ (· )和从H到H_k三个线性变换A_k(k=1, 2, 3), 从而对光谱维和空间维的交互信息进行融合, 得到如下基于交互信息的复合核函数表达式。

ϕ(xi)={A1φ(xis), A2φ(xiw), A3(φ(xis)+φ(xiw))}（8）

K(xi, xj)=< ϕ(xi), ϕ(xj)> =φ(xis)TR1φ(xjs)+φ(xiw)TR2φ(xjw)+φ(xis)TR3φ(xjw)+φ(xiw)TR3φ(xjs)（9）

R1=AT1A1+AT3A3, R2=AT2A2+AT3A3, R3=AT3A3（10）

写成复合核函数的形式为

K(xi, xj)=Ks(xis, xjs)+Kw(xiw, xjw)+Ksw(xis, xjw)+Kws(xiw, xjs)（11）

这里的复合核函数需要光谱维信息和空间维信息同维度才能使用, 即N_w=N_s。 根据选取的非线性映射函数ϕ (· )的不同, 一般有以下3种常见的核即线性核、 多项式核及径向基函数核。

K(xi, xj)=< xi, xj> K(xi, xj)=(< xi, xj> +1)dK(xi, xj)=exp(-‖xi-xj‖2/2σ2)（12）

SVM采用这些复合核函数以后, 可以明显提高高光谱遥感图像的地物分类精度^[10]。 但是, 这类方法并不能直接应用到高光谱显微图像的木材树种分类识别中。 这是因为, 在遥感影像地物分类中, 基本任务是将每个像素分类到合适的地物类别中(例如沙地、 河流、 耕地、 林地、 建筑物等类别), 它处理的基本数据单元是像素。 在这种情况下, 单个像素的光谱维特征是各个波长下的光谱反射率值, 而其空间维特征是该像素邻域窗口(例如3× 3, 9× 9邻域)的均值和标准差。

为了将基于复合核函数的SVM应用到木材树种显微高光谱图像的分类识别任务中, 需要作以下的转换工作。 首先, 这里处理的基本数据单元是手动选取的ROI区域, 具体步骤参见1.3.1节。 其次, 对于ROI区域的光谱维特征提取, 采用了光谱降维和计算ROI区域平均光谱值的方法, 参见1.3.1节和1.3.3节。 最后, 对于ROI区域的空间维特征提取, 采用了基于GLCM的纹理特征参数计算方法, 具体步骤参见1.3.4节。

为了准确评价分类模型的精度, 从4 800个样本中, 随机挑选3 000个样本作为训练集, 1 800个样本作为测试集。 分类器评价指标是由训练集分类准确率、 测试集分类准确率、 KAPPA系数(KC)和测试时间四个指标为评价标准。

2.1.1 原始光谱

利用Matlab2016b软件对光谱进行建模和分类, 图2所示是20种木材树种样本的平均光谱曲线图。 每种木材几乎都和其他树种类别有交叉, 光谱重叠严重, 而且部分光谱走势相似, 这对树种分类造成困难, 故必须对木材光谱数据做进一步处理。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 20种木材树种样本的平均光谱曲线Fig.2 Wood spectral curves of 20 wood species

2.1.2 采集波段建模分类

首先, 对采集波段光谱数据进行建模分类, 采用SVM和BP神经网络建立判别模型, 分类准确率如表1。 由表中可知, BP神经网络的测试时间少于SVM的测试时间, 但分类准确率和SVM相比较低。 SVM方法训练集和测试集分类准确率都在90%以上, 说明采集波段范围内利用SVM模型对木材树种分类是可行的。 SVM核函数采用了径向基核函数, 最终利用交叉验证法和网格搜索法确定惩罚因子C和径向基核函数中的参数σ 的值。

表1

Table 1

表1(Table 1)

表1 采集波段光谱建模分类对比 Table 1 Spectral classification comparisons for overall spectral band Model Training set/% Testing set/% KC Testing time/s SVM 98.800 0 92.222 2 0.918 1 0.405 2 BP 90.330 0 85.780 0 0.850 3 0.013 6

表1 采集波段光谱建模分类对比 Table 1 Spectral classification comparisons for overall spectral band

2.1.3 采集波段预处理建模分类

采用3种预处理方式对原始光谱进行预处理, 然后基于各种预处理光谱建立SVM模型。 从表2中可以看出, 通过不同预处理, 分类模型的精度得到了不同程度的提升。 采用SNV预处理的训练集和测试集分类准确率最高, 精度达到98.766 7%和93.000 0%, 而且所需测试时间也最短, 因此后续建模都采用这种方式。

表2

Table 2

表2(Table 2)

表2 不同预处理下采集波段光谱SVM模型分类结果 Table 2 Spectral classification comparisons for overall spectral band and different preprocessing schemes Pretreatment SVM Training set/% Testing set/% KC Testing time/s NONE 98.800 0 92.222 2 0.918 1 0.405 2 SNV 98.766 7 93.000 0 0.926 3 0.356 5 MSC 67.866 7 52.500 0 0.500 0 1.051 3 SG 98.900 0 92.277 8 0.918 7 0.379 0

表2 不同预处理下采集波段光谱SVM模型分类结果 Table 2 Spectral classification comparisons for overall spectral band and different preprocessing schemes

2.1.4 特征波长选取及建模分类

连续投影算法(SPA)是由Araujo等在以多元校正的背景下提出^[17], 通过对向量的投影分析, 确定投影向量最大的为待选波长, 再基于校正模型选择最终特征波长。 本实验在Matlab2016b软件中对经过SNV预处理后的样本进行特征波长挑选, 根据RMSE最小原则选取特征波长的集合。 如图3所示, 波长数目为19时RMSE=0.456 56接近最小值。 因此我们挑选19个波长, 分别是372, 382, 387, 392, 397, 412, 417, 432, 457, 580, 637, 699, 725, 762, 778, 800, 816, 961和1 000 nm。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 SPA波长数量与均方根误差关系图Fig.3 The relationship between wavelength number and RMSE

竞争性自适应重加权算法(CARS)^[18]是效仿达尔文进化论中“ 适者生存” 的原则提出的变量优选方法。 对采样次数反复迭代并寻找每次采样的RMSECV最小值, 所对应的变量视为优选出来的变量子集。 图4(a)表示选取变量数与运行次数的关系, 由图可知呈递减趋势。 图4(b)表示RMSECV的变化趋势, 当RMSECV变大表示剔除了有效信息; 当RMSECV变小表示剔除了无效信息。 图4(c)表示当运行次数为30时(中位线的位置), RMSECV最小。 反复进行参数选取调试, 筛选中, 蒙特卡洛采样次数设为100次。 最终CARS选择的特征波长分别是392, 397, 427, 432, 437, 447, 452, 462, 467, 549, 554, 564, 569, 590, 595, 605, 611, 616, 637, 642, 668, 673, 678, 683, 710, 741, 773, 842, 869, 885, 891, 902, 907, 912, 934, 950, 956和989 nm波长。 表3给出了采用这些特征波长进行SVM分类识别的结果。 综合对比两种特征波长选择方法, CARS的特征波长选择效果和运行速度较好一些。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 CARS方法波长筛选过程 (a): 选取变量数与运行次数的关系; (b): RMSECV的变化走势图; (c): RMSECV的最小值选取Fig.4 Spectral wavelength selection procedure with CARS scheme (a): Relationship between wavelength number and running times; (b): Varying trend chart for RMSECV; (c): Minimal value selection for RMSECV

表3

Table 3

表3(Table 3)

表3 采用特征波长SVM模型分类结果 Table 3 Classification results of SVM using feature wavelength Pretreatment Methods SVM Training set/% Testing set/% KC Testing time/s SNV SPA CARS 98.300 0 98.966 7 91.944 4 92.166 7 0.915 2 0.917 5 0.135 4 0.126 0

表3 采用特征波长SVM模型分类结果 Table 3 Classification results of SVM using feature wavelength

表4

Table 4

表4(Table 4)

表4 纹理建模分类比较 Table 4 Classification comparisons with wood texture features Model Training set/% Testing set/% KC Testing time/s SVM 99.800 0 60.333 0 0.577 2 0.128 8 BP 82.830 0 54.670 0 0.522 8 0.040 5

表4 纹理建模分类比较 Table 4 Classification comparisons with wood texture features

高光谱数据可表示为3D立方体, 不仅包含丰富的光谱维信息, 同时包含丰富的空间维/纹理信息。 采用灰度共生矩阵计算纹理特征。 先对高光谱图像进行主成分变换, 选取前3个主成分图像PC1— PC3, PC1图像保留很多原始图像信息, 而PC2和PC3的木材管孔射线等特征显示不清晰。 因此我们采用PC1对GLCM的纹理特征进行提取, 分别在0, 45, 90和135四个方向上的能量、 熵、 惯性矩、 相关性总共16个特征参数输入到SVM和BP神经网络中。 SVM的训练集和测试集分类精度最高达到99.800 0%和60.333 0%, 但测试时间却稍高于BP神经网络。 为了能准确对木材树种分类, 采用SVM对木材高光谱图像纹理进行建模分类。

将木材高光谱成像样本的38个特征波长下光谱值与16个纹理特征参数进行融合, 融合前对光谱和纹理特征做归一化处理, 融合后的特征参数作为输入数据进行复合核函数SVM建模分类, 识别结果如表5。

表5

Table 5

表5(Table 5)

表5 光谱和纹理信息融合分类结果 Table 5 Wood spectral-spatial fused classification comparisons Methods Composite kernels SVM Training set/% Testing set/% KC Testing time/s CARS+纹理 复合核1 100.000 0 93.333 3 0.929 8 0.180 0 CARS+纹理 复合核2 95.522 2 94.166 7 0.938 6 0.254 7 CARS+纹理 复合核3 93.766 7 92.610 0 0.922 2 0.260 2 PCA+纹理 复合核4 89.466 7 82.160 0 0.812 3 0.430 6

表5 光谱和纹理信息融合分类结果 Table 5 Wood spectral-spatial fused classification comparisons

从表5可知采用四种复合核函数后SVM建模分类正确率有所上升, 这是由于综合考虑了光谱维和空间维的分类特征。 第四个复合核函数分类精度较低, 这是因为第四核函数要求光谱与纹理同维限制, 为满足同维要求, 最终我们采用PCA方法, 因为通过对比CARS, SPA和PCA这3种降维方法选取16维光谱信息, PCA方法对本数据集的处理效率较高, 精度较好。 第二个复合核函数训练集和测试集分类效果最好可达到95.522 2%和94.166 7%。 综合分析, 本实验中CARS能够较好地用于可见光/近红外高光谱的选择, 该方法有效去除光谱中的无信息变量同时也能对共线性变量有效的剔除, 所以CARS+纹理模型优于其他模型, 同时也说明在高光谱图像上融合光谱和纹理特征对木材树种进行识别是可行的。 整体实验过程的流程图参见图5。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 实验过程的流程图Fig.5 Experimental flow graph