数据来源于2021年中国大学生数学建模竞赛E题的附件2。 其中, 总样本量673条, 分别来自于11个产地(数据类别标签), 缺失数据为5条, 没有重复样本。 原始光谱如图2所示, 其中横轴表示波长(波频段), 纵轴表示波吸收系数。 从图2中可以看出, 不同波段的原始数据横向变化大, 这说明数据独立特征属性多。 同一波段的不同样品特征高度相似, 表明建立产地鉴别模型难度大。 经过异常值分析, 发现第7号样本在波段0~1 500 nm范围内的测量值和第410号样本在波段2 500~3 000 nm范围内的测量值比其他样本的对应值偏大(如图2所示), 但是数值比较发现这两个样本的测量值在正常范围内, 所以最终使用的有效数据样本量为658条, 波长范围从551~3 998 nm。 原始数据中不同产地的样本量如表1所示, 可以看出样本分布相对均衡。 如果用矩阵格式表示该光谱数据(称作光谱数据矩阵), 则该数据有658行(样本个数, 表示来自不同产地的药材), 3 448列(数据属性或维度, 表示不同光谱波长的药材吸收系数), 还有1列为产地标签数据。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 原始光谱数据Fig.2 Raw spectral data

表1

Table 1

表1(Table 1)

表1 实验数据的样本产地分布 Table 1 Origin distribution of experimental data 产地 1 2 3 4 5 6 7 8 9 10 11 67 59 47 67 29 87 50 59 31 66 55

表1 实验数据的样本产地分布 Table 1 Origin distribution of experimental data

红外光谱数据包含了样品相关信息和与待测样品性质无关的干扰, 如样品的背景、 噪声、 杂散光以及仪器的响应等。 背景噪音影响待测样品的预测效果。 因此, 对红外光谱数据进行合适的预处理是建立稳健预测模型的关键和基础^{[1, 4]}。 文献[1]给出了光谱分析的预处理方法原理。 用于中药材产地鉴别的光谱数据, 按两种处理方式及组合展开:

(1)基于数据样本的预处理方式(记作光谱预处理)。 对于光谱数据矩阵, 第1类预处理方法基于每一行数据(样本), 即以每一个样品的不同频谱段的数据为处理单位开展。 考察了10种数据预处理方法(含不做数据处理), 其基本公式如下, 详细可参考文献[1]和[5]。

注: 下面各公式中, 令x={x_i}, i=1, 2, …, n表示行向量(样本)数据, x_i是其特征分量。 x^*={xi*}, i=1, 2, …, n表示x经过处理后的数据。

①最大最小值归一化处理

x*=x-xminxmax-xmin(1)

式(1)中: x_max和x_min分别为样本序列x的属性最大值和最小值。 此方法利用数据中的最大值和最小值实现对原始数据的等比例缩放至[0, 1]之间, 使得不同数量级之间的数据在数值上有可比性。

②中心化处理

x*=x-x¯(2)

式(1)中:x¯为样本序列x的平均值。 中心化也叫零均值化, 即原数据减去该组数据的平均值。 经过中心化处理后, 原数据的坐标平移至中心点(0, 0), 且中心化的数据均值变为0。 该方法的目的是消除数据量纲的影响。

③标准化处理

x*=x-x¯σ(3)

式(3)中:x¯为样本序列x的平均值; σ 为标准差。 此方法将数据按照比例进行缩放, 但不改变数据的原始概率分布。 缩放后的数据均值为0, 方差为1。 处理后的数据不一定符合标准正态分布, 但由于其更好地保持了样本间距, 且更符合统计学假设, 使得该方法提升了基于距离类型的算法模型精度。

④移动平均平滑(moving average, MA)处理

xi*=1N∑j=0N-1xi-j(4)

式(4)中:xi*的含义同式(3)。 移动平均平滑是最常用的滤波器, 其主要作用是降低数据噪声影响。

⑤SG平滑滤波(savitzky-golay filter)处理

xi*=1H∑j=-ω+ωxi+jhj(5)

式中:

H=∑i=-ω+ωhi(6)

式(6)中, H为归一化因子, h_i为平滑系数。 SG平滑滤波器是一种数字滤波器, 其能够在不改变信号趋势的前提下提高数据精度。

⑥正则化处理

xi*=xix(7)

式(7)中, 采用L1标准处理时,x=∑|xi|;当使用 L2 标准处理时,x=∑xi2;当使用max标准处理时, ‖ x‖ =max{x_i}。

⑦多元散射矫正(multiplicative scatter correction, MSC)^[8]处理:

记平均光谱计算公式

x̅=∑i=1nxin(8)

令线性回归计算公式

xi, LR=mix¯+bi(9)

则多元散射矫正的表达式

xi, MSC*=xi, LR-bimi(10)

式(10)中: m_i为线性回归方程斜率; b_i为线性回归方程截距。 通过最小二乘求解方法得到这两个参数。 MSC方法在光谱预处理中可以算作标准化方法^[5]。 通过调整m_i和b_i的大小可保留原有成分的信息, 使光谱数据的样本内特征差异变小。 该方法可以消除数据基线漂移和平移现象。

⑧一阶导数(the first derivative, 1st Der)处理

xi*=xΔx+i-xiΔx(11)

式(11)中: Δx为x_i处的变化量。 式(11)是红外光谱预处理中常用的一种方法^[5], 该种方法可以消除基线和降低背景干扰, 但是可能同时放大噪声。

⑨二阶导数(the second derivative, 2nd Der)处理

xi*=xi+Δx-2xi+xi-ΔxΔx2(12)

式(12)中, Δx的含义同式(11)。

(2) 基于数据样本特征的矩阵预处理(记作矩阵预处理): 这类(第2类)处理方式的数学公式与基于样本预处理(第1类)方法的公式相同。 只是该类数据处理方式是基于属性特征(矩阵列)的数据处理方式。 虽然基于导数类的光谱处理方法能够消除基线和背景干扰, 但会放大噪声、 降低数据信噪比^[5], 因此在矩阵预处理中不再考察^[5]。

①最大最小归一化: 式同式(1);

②中心化: 式同式(2);

③标准化: 式同式(3);

④正则化: 式同式(7)。

SVM算法是一类监督学习方法, 其基本模型是定义在特征空间上的间隔最大的线性分类器, 且其决策边界是对学习样本求解的最大边距超面。 如图3所示, SVM方法是实质上的非线性分类器。 该方法在最小化样本点误差的同时最小化结构风险, 进而提高了模型的泛化能力, 且没有限制数据维数。 因此, SVM算法在人像识别、 分类等模式识别问题中都取得了不错的应用效果。 这也是本工作选取多分类SVM方法作为检验数据预处理方法效果的主要原因。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 支持向量机模型示意图Fig.3 Schematic diagram of support vector machine model

选用了Python语言的机器学习库(Sklearn库)的SVM 方法。 对于每种数据预处理方法, 确定SVM模型中最优参数组合的K折交叉验证方法以及网格搜索策略如下:

K折交叉验证(K-fold cross validation): 在机器学习中, 当数据量不充分时采用的一种对数据进行重复使用方法。 其内容是把数据分为训练集、 验证集、 测试集, 每次实验随机选出n组数据, 训练出n个模型, 然后用测试集对n个模型进行评价, 根据得分选出最优模型。

网格搜索法^[9]: 是一种通过遍历参数组合来寻找优化模型的方法。 使用这种技术可以构建独立的模型, 并评估其性能, 从而找到产生最佳结果的模型和超参数。 比如, 在SVM参数寻优中, 以核函数RBF的SVM分类模型为例, 需要优化包括正则化常数和核函数参数在内的至少两个超参数, 并在验证集上分别评估它们的性能(或者在训练集内进行K折交叉验证), 从而得到结果。

决定系数(R²)是一种评价模型对数据拟合度的量。 见式(13)

R2=1-SSRSSTotal(13)

式中

SSR=∑(xi-xR)2(14)

其中, SS_R表示预测值和真实值的误差平方和。

SSTotal=∑(xi-x¯)2(15)

式(15)中, SS_Total表示所有样本点的平均值和真实值的误差平方和。

从式(13)可以看出: R²越接近于1, 则模型对数据的拟合度越高; 越接近于0, 则表示模型拟合精度越接近于平均值估计; 如果R²值小于0, 则拟合误差大于直接以平均值为模型估计值, 即模型拟合效果差。

基于输入的原始光谱数据及产地标签, 分别采用上述两类数据预处理方式及其组合处理后的数据进行SVM建模。 实验流程(见图4)为: 首先对原始数据进行10种基于样本(矩阵行)的数据预处理(即光谱预处理, 见1.2节), 再按照5种数据预处理方式对数据进行基于特征或维度的数据处理(即矩阵预处理, 见1.2节), 最终得到50组预处理后的数据, 共需要创建50个SVM模型。 具体的预处理方法和组合方式见表3。 预处理后的数据集按照训练集与验证集7∶3的比例划分。 最终得460个训练样本, 198个测试样本。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 实验流程图Fig.4 Experimental flow chart

考虑到机器学习算法对数据均衡性的要求, 记录了五次实验的训练样本(460条)中的不同产地样本数量以及平均值如表2与图5所示。 从中可以得到, 产地4与产地6的样本数较多, 而产地5与产地9的样本数量较少。 尽管实验数据样本量的不均衡情况不突出, 但考虑到预测模型的准确性, 采用了K折交叉验证方法来解决样本不均衡问题。

表2

Table 2

表2(Table 2)

表2 五次实验的训练集样本分布 Table 2 Sample number of origin obtained by five experiments 试验 次数 产地 1 2 3 4 5 6 7 8 9 10 11 1 44 42 48 60 23 51 34 45 20 50 43 2 50 40 42 59 16 67 42 41 19 47 37 3 49 39 50 62 22 53 36 41 23 44 41 4 45 35 50 58 22 60 19 42 19 47 45 5 46 45 45 58 23 64 36 45 20 45 36 均值 46.8 40.2 47.0 59.4 21.2 59.0 33.4 42.8 20.2 46.6 40.4

表2 五次实验的训练集样本分布 Table 2 Sample number of origin obtained by five experiments

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 五次实验的训练集产地统计图Fig.5 Sample number of origin in each five experiments

需要说明的是, 在按照7∶3的比例训练模型时, 对50种模型中的每一种都进行了5次实验。 不同模型的每次实验输入的数据是相同的, 表3给出了每种模型的5次数值实验平均结果。

表3

Table 3

表3(Table 3)

表3 预处理方式和对应的决定系数R2 Table 3 Preprocessing methods and corresponding coefficients of determination 光谱预处理 矩阵预处理 R2/% 1不处理 2中心化 3最大最小归一化 4标准化 5正则化 均值 1 不处理 73.84 74.22 33.99 77.04 59.42 63.70 2 最大最小归一化 84.58 78.10 56.68 74.42 -16.07 58.76 3 标准化 86.27 85.69 56.80 76.75 -27.14 61.10 4 中心化 0.51 6.94 -0.25 11.49 -121.51 3.79 5 移动平均平滑 70.34 74.06 29.15 77.04 61.32 62.38 6 SG平滑滤波 77.37 74.22 33.99 77.04 59.42 64.41 7 多元散射校正 76.57 83.26 46.42 70.55 -52.55 55.36 8 正则化 77.28 73.90 56.38 74.08 -42.98 56.33 9 一阶导数 69.80 73.18 79.60 86.59 79.51 77.74 10 二阶导数 0.00 0.00 80.86 94.66 93.56 53.82 均值 63.89 64.48 48.24 72.40 32.11 56.23

表3 预处理方式和对应的决定系数R2 Table 3 Preprocessing methods and corresponding coefficients of determination

对于表3中每一种预处理方法, SVM 最优超参数组合采用网格搜索法得到(见1.4部分), 然后把测试集输入到参数优选后的模型中, 得到决定系数[R², 计算方法见式(13)]。 每个模型进行5次随机实验, 并记录决定系数R²(每大类处理方法)的平均值。 上述过程均采用Python编程语言实现。

图6—图9展示了本次试验方案中的4种典型预处理后的数据特征。 从图6—图9可以看出, 4种不同方法预处理后的数据样本属性特征差异大, 特征间的相对关系变化大。 这意味着用不同方式预处理后的数据建立的模型预测精度可能不同。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 数据只做移动平滑处理Fig.6 Data with moving average processing

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 数据移动平滑+标准化处理Fig.7 Data pretreated with moving average+normalization

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 数据二阶导数+标准化处理Fig.8 Data pretreated with 2nd Der+normalization

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 数据二阶导数+最大最小化处理Fig.9 Data pretreated with 2nd Der+max-min normalization

表3给出了50种数据预处理建模的决定系数和每大类决定系数的平均值。

(1)从表3可以看出, 不同处理方法对应的决定系数(R²)变化范围从负值到正值。 其中, 最高预测值(94.66%) 对应的方法是基于样本的二阶导数处理+按列标准化处理。二阶导数按行处理+按列正则化处理的R²系数为93.56%。 这两种方法在所有的预处理方法中表现最好。

(2)第1行、 第1列对应的不做任何预处理的数据建模的R²值为73.84%, 超过平均值50%, 说明SVM算法使用原始数据也可以取得不错的效果。

(3)从第1列不做特征预处理对应的各行R²值来看, 按样本的最大最小归一化、 标准化方法对应的R²值近86%, 但所有基于样本处理方法的R²平均值为63.89%。 说明在只做基于样本的数据预处理中, 最大最小归一化、 标准化方法效果最好。

(4)从第1行(只做基于按特征的预处理)可以看出, R²平均值为63.70%, 这与只做各种按样本处理的R²平均值(63.89%)接近。

(5)R²系数出现了负数情况, 比如第5列(正则化列)的某些行的对应值。 这意味着先按样本行进行预处理, 再按特征进行正则化处理后的数据建模效果变差(决定系数远远低于不做任何预处理的数据建模效果, 如第1行、 第1列的73.84%)。

(6)从表3的第4行可以看出, 基于样本的中心化处理+按特征数据处理后的R²平均值最低, 为3.79%。 说明这种处理方法不利于建模。

(1)合适的数据预处理是非常必要的。 尽管药材光谱数据不存在量纲不一致问题, 但是如果不做预处理, R²值最高不到74%。 这主要因为不同光谱谱段的特征差异比较大。 当然, 不合适的预处理组合训练的模型, 会出现决定系数(R²)为负的情况。

(2)不同的预处理方法及组合训练的模型预测效果差异大。 从表3中可以看出, 决定系数(R²)变化范围从负值到最高正值为94.66%, 最高预测值对应的预处理方法组合是二阶导数按样本(行)处理+按特征(列)的标准化处理。 二阶导数按样本处理+按列进行的正则化处理对应的R²次高。 在某种意义上, 数据预处理组合后的SVM预测效果甚至超过了复杂模型的预测结果, 这对于降低建模复杂度有重要的意义。

(3)对于基于样本预处理+按特征(列)的中心化处理方法, R²出现了多个负值。 对于表中所出现的R²<0情况, 分析式(13)可知, SVM这时的拟合函数预测误差大于取Y=平均值函数的预测误差。

(4)先进行基于样本的一阶导数与二阶导数处理(表3的第9、 10行), 再进行最小最大归一化、 标准化或正则化数据预处理后的预测效果显著高于其他方法。 这主要是按样本的一阶导数处理突出了不同谱段变化较大部分, 消除了不同谱段的常数项差异; 二阶导数处理突出了相同样本的曲率, 消除了不同谱段的突变和常数项差异部分, 因此, 再做按特征的列标准化(或归一化、 正则化)可以消除各谱线基线差异不同带来的影响, 从而提高了预测率。 但是, 只做基于样本的二阶求导或按样本求二阶导数+按列中心化处理效果较差, 这是因为这种预处理突出了样本的特征变化大部分, 放大了特征内的噪音, 所以降低了预测准确率。

(5)只进行基于样本标准化处理(第3行)或者按样本标准化预处理+按列中心化处理对应的R²均超过85%。 这归功于预处理后的数据满足了正态分布假设要求, 因此提高了SVM预测率。

(6)从表3可以看到, 基于样本中心化+按特征列的正则化处理后的R²值为-121.5%。 这是因为样本间的特征值差异小而特征间的属性值差异大引起。 从图2可以明显看出不同特征的吸收系数变化大, 而样本之间表现为相似的特征行为(大小变化一致且差异小)。 因此先样本中心化, 只修改了样本的中心位置, 再特征正则化(把特征值统一修正为0~1之间)修改了较大的特征间差异, 降低了算法预测效果(低于不做任何预处理的R²值73.84%)。

(7)在组合预处理方法中, 所有基于样本预处理方法+标准化处理的SVM决定系数(R²)平均值最高, 为72.40%。 基于样本的一阶导数预处理+所有按列(特征)预处理的 SVM决定系数(R²)平均值最高, 为77.74%。 说明标准化和一阶导数预处理在光谱数据建模的稳定性中有意义。