玻璃文物特征分析及风化预测模型——基于机器学习与统计分析

Characteristic Analysis and Weathering Prediction Model of Glass Cultural Relics—Based on Machine Learning and Statistical Analysis

摘要: 本文基于2022年高教杯C题，首先将数据集行one-hot编码，对化学成分数据进行填补缺失值，并使用Person相关系数对数据相关性进行分析。后基于奇异值分解的PCA算法对数据降维。通过随机森林、支持向量机、Xgboost、Logistic回归对降维后的数据进行分类并求得决策边界，基于软分类器预测出了各个文物的风化程度，得到在文物信息数据的验证集分类准确度达到了86.7%，化学成分数据达到了94.1%。之后利用召回率、f1值、精度、ROC曲线等评价指标对模型进行了评价与选择。最终得到了各特征的相关性以及各文物的风化概率，同时得出了Xgboost在该数据集中预测的优越性。

Abstract: Based on the 2022 Higher Education Cup A question, this paper first codes the data set line one-hot, fills in the missing values for the chemical composition data, and analyzes the data correlation using the Person correlation coefficient. After the PCA algorithm based on singular value decomposition reduces the dimensionality of the data. Through random forest, support vector machine, Xgboost, logistic regression to classify the dimensionality reduction data and find the decision boundary, based on the soft classifier, the degree of weathering of each cultural relics is predicted, and the classification accuracy of the verification set of the cultural relics information data reaches 86.7%, and the chemical composition data reaches 94.1%. After that, the model was evaluated and selected using evaluation indicators such as recall rate, f1 value, accuracy, and ROC curve. Finally, the correlation of each feature and the weathering probability of each artifact are obtained, and the superiority of Xgboost’s prediction in this data set is obtained.

1. 引言

古代玻璃玻璃制品是研究早期丝绸之路各国贸易往来的重要文物，具有非常宝贵的研究意义 [1] 。但由于各国制造玻璃制品的取材不同致其化学成分各不相同，且受掩埋环境影响大，容易风化，内部与环境元素进行大量交换，玻璃制品的成分比例发生变化，故影响对其真实类别的判断 [2] ，因此，如何对风化后的玻璃制品进行成分分析与鉴别，从而还原其研究价值是较有意义的研究课题。姜中宏、张勤远 [3] 利用同位素标记法分析铅钡玻璃的化学成分；赵凤燕、陈斌等人 [4] 利用PXRF分析研究文物成分；李鹏艳、谢承利 [5] 等人利用激光诱导击穿光谱法分析玻璃成分。本文提出使用Person相关系数对数据相关性进行分析，后基于奇异值分解的主成分分析算法对数据降维。通过随机森林、支持向量机、Xgboost、Logistic回归对降维后的数据进行分类并求得决策边界，通过一系列指标得出最优模型，并基于软分类器预测出了各个文物的风化概率。

2. 数据处理

2.1. One-Hot编码

One-Hot 编码，又称为一位有效编码，主要是采用N位状态寄存器来对N个状态进行编码，每个状态都由他独立的寄存器位，并且在任意时候只有一位有效，其数学的本质在于将文本或字符串映射到 n 维二进制空间上words→Rⁿ [6] 。

One-Hot编码是分类变量作为二进制向量的表示。这首先要求将分类值映射到整数值。然后，每个整数值被表示为二进制向量，除了整数的索引之外，它都是零值，它被标记为1。则：对于纹饰A、B、C有： $y\in \left\{\left(1,0,0\right),\left(0,1,0\right),\left(0,0,1\right)\right\}$ 其中A、B、C分别对应集合当中的三行向量。对于玻璃类型高钾类、铅钡类，可以把特征映射到以下集合： $y\in \left\{\left(1,0\right),\left(0,1\right)\right\}$ 其中铅钡玻璃对应，高钾玻璃对应。同理，对于玻璃颜色，一共有浅绿、浅蓝、深绿、深蓝、紫、绿、蓝绿、黑共8种颜色，可以将其映射到对应编码为：

$\begin{array}{l}y\in \{\left(1,0,0,0,0,0,0,0\right),\left(0,1,0,0,0,0,0,0\right),\left(0,0,1,0,0,0,0,0\right),\left(0,0,0,1,0,0,0,0\right)\\ \left(0,0,0,0,1,0,0,0\right),\left(0,0,0,0,0,1,0,0\right)\left(0,0,0,0,0,0,1,0\right),\left(0,0,0,0,0,0,0,1\right)\}\end{array}$

其中，上述行向量为1的位置依次对应上述各颜色，这里不再赘述。通过One-Hot编码，将离散特征数字化，为后文相关性分析以及机器学习模型的预测做准备。

2.2. 特征相关性分析

Person相关系数是用来判断两个变量是否存在线性相关性的一种数学模型，其原理是通过求两个变量的协方差从而得到相关系数 [7] ，本文一一将两个变量做线性相关性分析，由于只考虑各特征与是否风化的相关性，得到各个变量与风化与否的相关系数矩阵相关矩阵如表1所示。协方差：

$Cov\left(X,Y\right)=\frac{{\displaystyle {\sum }_{i=1}^n\left(X_i-E\left(X\right)\right)\left(Y_i-E\left(Y\right)\right)}}{n}$ (1)

其中，E(X)、E(Y)表示两个变量的期望，由于只需要研究各个变量同风化与否的离散变量的关系，因此对于相关系数矩阵这个方阵，只需要求出最后两行所代表的值进行具体分析。

Table 1. A matrix of correlation coefficients for each variable and weathering

表1. 各个变量与风化的相关系数矩阵

由以上相关系数矩阵可以看出，纹饰A与纹饰B分别与风华呈现弱负相关性和弱正相关性，纹饰C则无影响，铅钡与风化呈现弱正相关性，但高钾恰好相反，在颜色方面，只有浅绿色、深蓝色、绿色和黑色对风华具有影响，推测是因为化学成分的不同，因此对风化于否的影响不同，而又由于颜色由化学成分决定，纹饰又和颜色关联度高，因此推测各个特征之间具有线性相关性，因此该数据还需要进行降维处理。

2.3. 基于奇异值分解的主成分分析降维模型

对于我们所建立的稀疏矩阵特征，由于变量间存在线性相关，并且计算成本较大，因此需要算法来对高维的离散数据进行降维，本文引入主成分分析法对数据进行处理。

2.3.1. 奇异值分解

对于一个m × n的任意矩阵A，在本数据集的划分当中，离散特征转化后数据集一m = 58, n = 13，对于数据集二来说，m = 67, n = 16，也就是说一个存在13个离散特征，另一个数据集存在16个连续特征，并且特征与目标变量之间线性相关性普遍不强，特征之间存在更强的线性相关，因此对于这个数据矩阵A，可以做奇异值分解处理 [8] ，来提取它的主成分：

$A=X_m{\delta }_n{\left(X_n\right)}^{\text{T}}$ (2)

其中 $X_m$ 为m阶方阵， ${\delta }_n$ 为对角线为 ${\delta }_1,{\delta }_2,\cdots ,{\delta }_n$ 的对角矩阵， ${\left(X_n\right)}^{\text{T}}$ 为n阶方阵转置，上述矩阵可化简为：

$A=U\Sigma V^{\text{T}}$ (3)

其中U，Σ是由A^TA的特征值和特征向量得到，V是由AA^T的特征向量得到。

2.3.2. 主成分分析降维模型

针对奇异值分解算法当中的V^T是一个n矩阵，其中前i个列，对应 $\delta$ 矩阵当中前i个奇异值，也分别可以作为该数据集的主成分，可以认为主成分分析的本质就是奇异值分解。一、二数据集的奇异值分解如图1所示。

利用sklearn中的PCA求解器进行计算，发现在数据集一当中得到了六个主成分所包含整个数据集的信息量分别为36.7%、26.2%、13.2%、6.51%、5.00%、3.71%，共计能够包含91.5%的信息；数据集二当中能够提取到三个主成分的信息量分别为84.2%、9.38%、3.73%，共计97.3%的信息。

Figure 1. The amount of retention of the dimensionality reduction difference between the two datasets obtained by principal component analysis

图1. 主成分分析所得到的两个数据集的降维后方差的保留量

显然，对于数据集二的降维的效果更好，原因是问题一当中数据分布较为稀疏，且矩阵值为离散变量，数据二当中的数值较为连续，且较为密集，因此可以更好地利用奇异值分解提取主成分。图2是数据集二的数据降维可视化，可以看出数据二中92.6%特征降维到PC1和PC2上。

3. 基于不同机器学习算法的风化概率预测模型

机器学习因其优秀的数据分析能力以及模型的泛化能力在回归预测、分类等问题上表现出优越性。由于本文的研究是基于离散数据得出文物是否风化的二分类问题，为了得到连续化数据，本文将二分类问题转化为求解风化概率的连续性问题，通过将是否风化映射到[0, 1]并设置0.5作为分界线，利用回归模型预测出一系列风化概率。

Figure 2. Visualization of PCA reduction results in dataset two

图2. 数据集二降维可视化

针对文物风化的概率，本文分别使用Logistic回归、极致梯度提升(Xgboost)、随机森林、软间隔支持向量机模型对降维后的数据分类并求出决策边界，基于软分类器求出各文物风化的概率。然后利用召回值、F1值、精度、ROC曲线等评价指标选出最好的预测模型。机器学习使用经降维后的数据集，两个数据集分别提取到了6、3个主成分。对于数据集的划分，本文选择75%的数据作为训练集，25%的数据作为测试集。

3.1. Logistic回归模型

Logistic回归是统计学中重要的经典分类方法，属于对数线性模型，所以也被称为对数几率回归。将对数概率函数取反后得到sigmoid曲线，首先利用极大似然函数拟合出参数的最佳取值，然后将降维后的特征输入分类模型得到各文物被风化的概率，由求出的参数拟合决策边界从而将连续的风化概率离散为是否风化两种情况。

设P为某件文物发生风化的概率，则1 − P为未被风化的概率。显然，概率是非负的，为了将概率记为输入特征值的线性表达式，将概率取对数，得到：

$\begin{array}{c}\ln \left(y=1|x\right)=w_0{x}_0+w_1{x}_1+\cdots +w_n{x}_n\\ ={\displaystyle \underset{i=0}{\overset{n}{\sum }}{w}_i{x}_i}\\ =W^{\text{T}}x\end{array}$ (4)

上述公式表示当输入特征x后，文物被分为风化的概率，其中为决策边界系数矩阵。将对数函数取反，令Z = X，可以得到：

$f\left(z\right)=\frac{1}{1+{\text{e}}^{-z}}$ (5)

该函数就是sigmoid函数，为了实现logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和带入sigmoid函数中。进而得到一个范围在0~1之间的数值。最后设定一个阈值，在大于阈值时判定为1，否则判定为0。

关于参数W求解，本文以极大似然函数为代价函数，分别让代价函数对W中j个权重求偏导，沿着梯度下降的方法，找到使极大似然函数取最小值时的参数W [9] 。

极大似然函数L(w)：

$\begin{array}{c}L\left(w\right)=P\left(y|x;w\right)\\ ={\displaystyle \underset{i=1}{\overset{n}{\prod }}P\left(y^{\left(i\right)}|x^{\left(i\right)};w\right)}\\ ={\left(\varphi \left(z^{\left(i\right)}\right)\right)}^{y^{\left(i\right)}}{\left(1-\varphi \left(z^{\left(i\right)}\right)\right)}^{1-y^{\left(i\right)}}\end{array}$ (6)

将极大似然函数取对数，将连乘转化为累加求和的形式，并对j个权重求偏导：

$\frac{\partial }{\partial w_j}J\left(w\right)=\left(-y\frac{1}{\varphi \left(z\right)}+\left(1-y\right)\frac{1}{1-\varphi \left(z\right)}\right)\frac{\partial }{\partial w_j}\varphi \left(z\right)$ (7)

然后对sigmoid函数求偏导：

$\begin{array}{c}\frac{\partial }{\partial w{}_j}\varphi \left(z\right)=\frac{\partial }{\partial z}\frac{1}{1+{\text{e}}^{-z}}\\ =\frac{1}{{\left(1+{\text{e}}^{-z}\right)}^2}{\text{e}}^{-z}\\ =\frac{1}{1+{\text{e}}^{-z}}\left(1-\frac{1}{1+{\text{e}}^{-z}}\right)\\ =\varphi \left(z\right)\left(1-\varphi \left(z\right)\right)\end{array}$ (8)

将(8)式代入(7)式即得到了函数对j个权重的偏导，并按照式(9)的规则不断更新权重，最终求出参数W。

$w:=w-\eta \nabla J\left(w\right)$ (9)

在logistic回归中，决策边界方程定义为W^Tx = 0，上文得到假设函数：

$P\left(y=1|x;w\right)=g\left(w^{\text{T}}x\right)=\frac{1}{1+{\text{e}}^{-w^{\text{T}}x}}$ (10)

通过计算得到数据集一、二的决策边界函数分别为：

$h_{\theta 1}\left(x\right)=g\left(-0.44x_1+0.05x_2-0.48x_3-1.87x_4-1.27x_5-0.33x_6\right)$ (11)

$h_{\theta 2}\left(x\right)=g\left(-0.0319x_1-0.0447x_2-0.0743x_3\right)$ (12)

决策边界也可以是一个软边界分类器，即可以分析风化程度的概率，因此在logistic回归当中，把数据代入决策函数，可以得到该行数据是否风化的概率，由于数据集是离散的，求概率的过程相当于利用提取到的特征对离散数据连续化处理，因此可以得到每个文物的风化概率，如果风化概率小于0.5则可以说明该文物没有风化。图3是由logistic回归得到的部分文物的风化概率。

3.2. 软间隔支持向量机模型(Support Vector Machines, SVM)

支持向量机是一类通过求解最大边距超平面进而解决一个二元分类问题，而本数据集的非线性性质较高，无法做到线性可分，因此需要进行一个软间隔的支持向量机分类模型。

Figure 3. Based on Logistic regression of the weathering probability of some cultural relics

图3. 基于Logistic回归部分文物风化概率

支持向量机本身是一个约束凸二次优化问题，其目标函数为 $\frac{1}{2}{\Vert \omega \Vert }^2+C{\displaystyle {\sum }_{i=1}^N{\epsilon }_i}$。

其中 $\Vert \omega \Vert$ 是支持向量机模型中计算的超平面的L2范数；C > 0 称之为惩罚参数，当C较大的时候，惩罚更大，避免过拟合效果更好，但过大会导致分类精度较低；而C较小的时候，容易过拟合，但分类精度大，因此需要在模型验证的时候选择合理的参数C。而 ${\epsilon }_i$ 是一个松弛变量。针对线性不可分的支持向量机模型可以写作一个凸二次规划问题 [10] ：

$\begin{array}{l}\text{max}{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\alpha }_i-{\displaystyle \underset{i,j=1}{\overset{m}{\sum }}{\alpha }_i}}{\alpha }_j{y}_i{y}_j\left(x_i-x_j\right)\\ \text{Subjectto}:{\displaystyle \underset{i=1}{\overset{m}{\sum }}{\alpha }_i{y}_i=0}\\ 0\le {\alpha }_i\le C,i=1,\cdots ,m\end{array}$ (13)

3.3. 随机森林模型(Random Forest, RF)

随机森林是一个bagging法的集成学习方法，每次对数据集进行随机特征采样，再创建一个简单的决策树模型，对各个决策树的模型的预测结果进行投票，得到一个最终的结果。其思想如下所示：

Input：分类数据集X，Y

Output：投票生成的类标

1：for i=1to B do

2：对数据集进行随机采样，在数据与类标X，Y当中采样n个，称为X_b，Y_b

3：训练弱模型f_B

4：end for

5：对于每个弱模型，其预测的算数平均数 $\hat{h}=\frac{1}{B}{\displaystyle {\sum }_{b=1}^B{f}_b(\; x\; \prime \; )}$

3.4. 极致梯度提升模型(Xgboost)

Xgboost是一种boosting算法，其目的在于利用迭代优化，对上一次迭代的预测结果进行加权，使得误差越来越小。而Xgboost算法是一种由k个基模型进行叠加的一个带有正则化项 $\omega$ 的目标函数，利用泰勒公式将目标函数展开后得到 [11] ：

$\begin{array}{c}Obj={\displaystyle \underset{i=1}{\overset{n}{\sum }}l\left(y_i,{\hat{y}}_i\right)}+{\displaystyle \underset{i=1}{\overset{k}{\sum }}\Omega \left(f_t\right)}\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}l\left(y_i,{\hat{y}}_i^t\right)}+{\displaystyle \underset{i=1}{\overset{t}{\sum }}\omega \left(f_i\right)}\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}l\left(y_i,{\hat{y}}_i^{t-1}+g_i{f}_t\left(x_i\right)\right)}+{\displaystyle \underset{i=1}{\overset{t}{\sum }}\Omega \left(f_i\right)}\\ ={\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(l\left(y_i,{\hat{y}}_i^t\right)+g_i{f}_t\left(x_i\right)+\frac{1}{2}{h}_i{f}_t^2\left(x_i\right)\right)}+{\displaystyle \underset{i=1}{\overset{t}{\sum }}\Omega \left(f_i\right)}\end{array}$ (14)

其中g_i, f_i分别是前一项损失函数的一阶导和二阶导。在决策树模型的Xgboost当中，由于决策树的复杂度可以通过叶子数T确定，因此目标函数的正则化项可以写为：

$\Omega \left(f_t\right)=\gamma T+\frac{1}{2}\lambda {\displaystyle \underset{j=1}{\overset{T}{\sum }}{w}_j^2}$ (15)

经简化，目标函数可写为：

$Obj={\displaystyle \underset{j=1}{\overset{T}{\sum }}\left[\left({\displaystyle \underset{i\in I_j}{\sum }{g}_i}\right)w_j+\frac{1}{2}\left({\displaystyle \underset{i\in I_j}{\sum }{h}_i}+\lambda \right)w_j^2\right]}+\gamma T$ (16)

3.5. 模型评价及选择

本文选取了logistic回归、随机森林、支持向量机、Xgboost模型分别对数据进行拟合分类，通过召回率Recall、f1值，support，精确度，AUC，ROC曲线进行模型的评价，几项指标的定义如下 [12] (其中TP、FN、TN、FP代表True Positive、False Nagetive、True Nagetive、False Positive，前两者指的是预测正确的情况，后两者指的是预测反的情况)：(其中AUC指ROC曲线包围的面积)

$\text{accuracy}=\frac{TP+TN}{TP+FN+TN+FP}$ (17)

$\text{precision}=recall=\frac{TP}{TP+FP}$ (18)

$\text{F1}=\frac{2\ast precision\ast recall}{precesion+recall}$ (19)

$AUC=S\left(Roc\right)$ (20)

Table 2. Dataset 1 Learning performance of each model

表2. 数据集一各模型学习表现

上表2展示了各模型对数据集一的预测情况，可以看出：线性模型和树模型对数据一的拟合效果都比较好，最高的拟合精度在验证集上能够达到86.7%，对于机器学习任务来说，可以认为该数据集的拟合效果比较好，并且不存在过拟合。随机森林算法的精度能够达到1.00，相当于预测正类全部预测正确。

Table 3. Dataset 2 Learning performance of each model

表3. 数据集二各模型学习表现

针对数据集二，利用上述指标评价各个模型得到如表3所示数据。仍然能够发现数的模型表现要比线性模型要好，并且无论是随机森林还是Xgboost和SVM，表现均优于数据集一，但是发现线性模型的预测能力反而有所下降。因此可以得出结论：数据二数据的非线性性较强，由于随机森林可以叠加不同分类器，Xgboost可以自动抽取非线性特征，而SVM可以通过内核进行无穷升维，进而对非线性性有一个非常好的拟合效果，Xgboost的效果在非线性数据集当中效果是最好的。

机器学习还有另外一大评价指标：ROC曲线，也叫接收者操作特征曲线。在ROC曲线当中，曲线与横轴的面积即为AUC值，面积越大代表模型的分类效果越好。根据图4的ROC曲线，可以看出分类效果随机森林最好，Xgboost次之。支持向量机在分类当中的表现略弱。

Figure 4. ROC curves for each model, from left to right, logistic regression, support vector machine, random forest, Xgboost

图4. 各个模型的ROC曲线，从左到右依次为logistic回归、支持向量机、随机森林、Xgboost

对于各个模型的AUC曲线，把正类预测对的概率为红色曲线，预测错的曲线为蓝色曲线，评价指标是红色曲线和蓝色曲线之间的面积差，可以看出来对于各个模型来说，红色曲线数值是蓝色曲线的一个上界，证明了机器学习算法的有效性。由以下结果可以看出，基于树的模型分类效果较好，线性模型表现优于支持向量机，Xgboost表现是最好的。

从图5可以看出红色曲线与蓝色曲线之间的面积达到最大的时候，分类正确次数最多，分类错误次数最少，Xgboost的表现最好。

针对数据集二，利用上述指标评价各个模型。任然能够发现数的模型表现要比线性模型要好，并且无论是随机森林还是Xgboost和SVM，表现均优于数据集一，但是发现线性模型的预测能力反而有所下降。因此可以得出结论：数据二数据的非线性性较强，由于随机森林可以叠加不同分类器，Xgboost可以自动抽取非线性特征，而SVM可以通过内核进行无穷升维，进而对非线性性有一个非常好的拟合效果，Xgboost的效果在非线性数据集当中效果是最好的。

Figure 5. The AUC curves for each model, from left to right, logistic regression, support vector machines, random forests, Xgboost

图5. 各个模型的AUC曲线，从左到右依次为logistic回归、支持向量机、随机森林、Xgboost

3.6. Xgboost效果解释

由下图6可见，Xgboost的ROC曲线下的面积为1，并且AUC曲线当中红色曲线和蓝色曲线之间的面积差值非常大，根据概率预测和正负类划分的预测图表明，Xgboost分类的时候仅仅会把极少数“模棱两可”的数据误分类。

Figure 6. Xgboost’s ROC curve, AUC curve, and classification positive and false label curves

图6. Xgboost的ROC曲线、AUC曲线与分类正误标签曲线

对于Xgboost在数据集上表现良好的原因，可以总结为：

1) 鲁棒性强，对异常点不敏感，不需要归一化

2) Xgboost的决策流形可以看做超平面的决策边界

3) 利用贪心算法、自动特征选择，可以提取到更多的非线性特征

4) 数据量带来的增益不大，适合小样本学习

通过上文四种模型对比得出Xgboost在预测模型中的优越性，得到了如下表4所示的部分风化概率：

Table 4. Some artifact prediction data based on Xgboost

表4. 基于Xgboost的部分文物预测数据

4. 结果分析

本文首先通过person相关系数矩阵得出各特征之间可能存在相关性并使用基于奇异值分解的PCA降维将两个数据集维数降低。将处理后的数据按照训练集与测试集3:1的比例划分训练Logistic回归、随机森林、软间隔支持向量机和Xgboost四种模型，训练结果表明：数据集一线性相关较强，数据集二的非线性相关较强，三随机森林、支持向量机和Xgboost三种非线性模型中Xgboost在数据集二中表现更好。

根据本文结论，后续研究可以利用Xgboost模型对数据集二进一步分析并依据相关特征划分文物类型。然后根据划分出的类型，通过灰色关联分析、线性模型和多层感知机模型分别拟合文物中各化学元素的关联性。

NOTES

^*通讯作者。