古代玻璃制品成分分析和鉴定的研究

Study on Composition Analysis and Identification of Ancient Glass Products

摘要: 古代玻璃种类繁多且易受环境影响而风化，因此需要对古代玻璃制品的化学成分数据分析，研究有无风化玻璃制品成分的变化规律，以并探索亚分类方法，进而可以根据未知分类的文物化学成分对文物进行准确的分类。本文通过使用K-means算法和BP神经网络结合的方式对玻璃制品进行亚分类划分，之后根据亚分类种类进行风化前后成分的预测；通过RUSBoost机械学习算法，70%的数据作为训练集，15%的数据作为测试集，其余部分作为预测集，来进行玻璃制品的种类鉴定。这些模型相互之间配合紧密，所得结果依次递进，使最终求解真实可靠。模型充分联系实际，具有很好的通用性和推广性。

Abstract: There are many kinds of ancient glass and it is easy to be weathered by the environment. Therefore, it is necessary to analyze the chemical composition data of ancient glass products, study the chang-ing rules of the composition of weathered glass products, and explore subclassification methods, so as to accurately classify cultural relics according to the chemical composition of cultural relics un-known. This paper uses K-means algorithm and BP neural network to subclassify glass products, and then predicts the composition before and after weathering according to the subclassification types. Through the RUSBoost mechanical learning algorithm, 70% of the data is used as the training set, 15% of the data is used as the test set, and the rest of the data is used as the prediction set to identify the types of glass products. The classification rules and identification of ancient glass prod-ucts are divided. These models cooperate closely with each other, and the results are progressive in turn, which makes the final solution true and reliable. The model is fully connected with practice and has good generality and popularization.

1. 引言

丝绸之路是古代中西方文化交流的通道，其中玻璃是早期贸易往来的宝贵物证。古代玻璃种类繁多且易受环境影响而风化，因此需要对古代玻璃制品的化学成分数据分析，研究有无风化玻璃制品成分的变化规律，以及高钾、铅钡两种玻璃类型的化学成分统计规律，并探索亚分类方法，进而可以根据未知分类的文物化学成分对文物进行准确的分类。

由于玻璃文物表面风化与其类型、纹饰和颜色具有一定关系。首先确定各因素对表面风化的影响程度。之后采用定量分析，根据化学成分将玻璃文物进行分类，并根据不同玻璃类型进行风化前后的成分分析。关于对玻璃文物风化前的成分预测，本论文采用先分类后预测的方法，提高预测准确性。通过RUBoost机械学习算法根据玻璃文物的化学成分推断玻璃文物的类别。

2. 基于K-Means算法和BP神经网络的亚分类方法

2.1. 实验数据及处理

本文研究所用的实验数据包括三个：表1：58个文物的纹饰、类型、颜色和表面风化情况的信息；表2：58个文物不同部位的化学成分；表3：8个已知化学成分和风化情况，但未知玻璃类型的文物。

将表单2中的化学成分百分比之和大于100%和小于85%的数据去除，将空白项进行填零处理。

表单1部分数据：

Table 1. Form 1 partial data

表1. 表单1部分数据

表单2部分数据：

Table 2. Form 2 partial data

表2. 表单2部分数据

表单3部分数据：

Table 3. Form 3 partial data

表3. 表单3部分数据

2.2. 相关性分析

为判断玻璃文物表面风化与其类型、纹饰和颜色是否具有关系，需要在计算皮尔逊相关系数前计算各个变量之间的显著关系p值。两者之间的相关关系可能只是偶然因素引起的，所以我们要对两个变量之间的相关关系的显著性水平进行判断。

将表1中的文字数据替换为数值矩阵形式。首先对对两变量之间是否具有统计上的显著关系(p值严格小于0.01)进行检验。通过表4可以得出纹饰对表面风化和颜色对表面风化无显著关系。之后在计算每两个变量之间的皮尔逊系数，之后画出相关性热力图。

${\rho }_{X,Y}=\frac{cov\left(X,Y\right)}{{\sigma }_X{\sigma }_Y}=\frac{E\left(\left(X-{\sigma }_Y\right)\left(Y-{\sigma }_X\right)\right)}{{\sigma }_X{\sigma }_Y}=\frac{E\left(XY\right)-E\left(X\right)E\left(Y\right)}{\sqrt{E\left(X^2\right)-E^2\left(X\right)}\sqrt{E\left(Y^2\right)-E^2\left(Y\right)}}$ (4)

Table 4. P values

表4. P值

通过图1我们不难看出：玻璃类型对于表面风化和颜色具有较强的相关性。根据古代玻璃的制作，是根据不同的玻璃类型和颜色确定纹饰进行加工的，因此纹饰应该与颜色和玻璃类型呈负相关，即颜色和玻璃类型影响纹饰。

Figure 1. Correlation thermal map

图1. 相关性热力图

2.3. 数据分析

首先我们根据表2中的数据将玻璃分为高钾和铅钡对其化学成分均值以三维柱状图形式呈现。

通过图2我们可以发现高钾玻璃主要是由大量二氧化硅元素组成；铅钡玻璃是由基本等量的二氧化硅和氧化铅元素以及部分氧化钡元素组成，因此可以通过二氧化硅和氧化铅以及氧化钡元素含量来对高钾玻璃和铅钡玻璃进行区分。

Figure 2. Three-dimensional bar chart of high potassium lead and barium content

图2. 高钾铅钡元素含量三维柱状图

2.4. 亚类划分及预测

在高钾和铅钡玻璃的前提下，进行亚类分析。通过K-Means和BP神经网络在高契合度的情况下，兼顾分类的多样性和准确性下，得到K = 8时的效果最好，BP神经网络的检测准确度为95.5%。

2.4.1. 模型建立

目前经常被人们使用的检验聚类有效性的指标有XB指标(Xie-Beni)、KL指标(Krzanowski-Lai)、Sil指标 (Silhouette)、DB指标(Davies-Bouldin)、IGP指标(In-Group Proportion)、BWP指标(Between-Within Proportion)等。其中BWP 是一种有效的聚类指标，该指标依据簇内相似与簇间相异两方面来评判聚类结果 [1] [2]。如果 $X=\left\{X_1,X_2,\cdots ,X_n\right\}$ 为聚类数据集，k = (X, R)为聚类空间，如果数据集中的n个样本分为c类，则最小类间距离、类内距离，BWP指标判别函数公式如下。

最小类间距离：

$b\left(i,j\right)=\underset{1\le m\le c,m\ne j}{\min }\left(\frac{1}{n_m}{\displaystyle \underset{p=1}{\overset{n_m}{\sum }}{\Vert x_p^m-x_i^j\Vert }^2}\right)$ (2)

类内距离：

$w\left(i,j\right)=\frac{1}{n_j-1}{\displaystyle \underset{q=1,q\ne i}{\overset{n_j}{\sum }}{\Vert x_p^m-x_i^j\Vert }^2}$ (3)

BWP指标：

$\text{BWP}\left(k\right)={\displaystyle \underset{j=1}{\overset{k}{\sum }}{\displaystyle \underset{i=1}{\overset{n_i}{\sum }}\text{bwp}\left(i,j\right)}}$ (K为类簇， $\text{bwp}\left(i,j\right)=\frac{b\left(i,j\right)-w\left(i,j\right)}{b\left(i,j\right)+w\left(i,j\right)}$ ) (4)

其中BWP指标值越高，聚类性能越好，BWP指标确定k值。

对于单个样本来说，BWP指标值越高，则意味着聚类性能就越好。对大量样本来说，取BWP的平均值来观察聚类性能，BWP平均值越高，数据集的聚类性能就会越好。其中 ${\text{avg}}_{\text{BWP}\left(k\right)}$ 为所有样本的BWP平均值， $k_{opt}$ 为聚类数，计算公式如下：

${\text{avg}}_{\text{BWP}\left(k\right)}=\frac{1}{n}{\displaystyle \underset{j=1}{\overset{k}{\sum }}{\displaystyle \underset{i=1}{\overset{n_i}{\sum }}\text{bwp}\left(i,j\right)}}$ (5)

$k_{opt}=\arg \underset{2\le k<n}{\max }\left({\text{avg}}_{\text{BWP}\left(k\right)}\right)$ (6)

本文提出一种基于最大最小距离和BWP指标相结合的K-Means算法 [1] [2]，以下命名为UPK-Means，改进后算法的基本思路如下：假如数据集 $X=\left\{X_1,X_2,X_3,\cdots ,X_n\right\}$ 可分为K类，K数目未知。首先在X中找出距离最远的两个样本点C₁、C₂，结合给定的θ值，根据最大最小距离原则确定余下的K-2个聚类中心，然后把样本点分配给距离最近的聚类中心，当数据集中的每个样本点都有确定的一个类别时，聚类结束。参数值θ不同，聚类效果也会有所不同，最后利用 ${\text{avg}}_{\text{BWP}\left(k\right)}$ 评价聚类性能， ${\text{avg}}_{\text{BWP}\left(k\right)}$ 值越高，聚类性能越好。

2.4.2. 模型求解

将表2中处理后的数据导入SPSSPRO中，使用K-Means聚类算法后得出玻璃文物的亚分类规律。之后将各分类的化学成分占比以图3形式表示。

Figure 3. Partial results of subclassification

图3. 亚分类的部分结果

在根据分化点检测的数据预测风化前的化学成分之前，需要将该风化点根据检测数据划分到不同亚类当中，之后根据再根据该亚类的风化前后化学成分变化规律来进行预测。

风化前后变化率 $\eta$

$\eta =\frac{{\alpha }_{未风化高钾}-{\alpha }_{风化高钾}}{{\alpha }_{未风化高钾}}$ (7)

由于铅钡玻璃中二氧化硅，氧化铅和氧化钡含量在各亚类中无明显差别，因此为了提高预测精度，这几项含量采用所有铅钡玻璃进行变化率的计算；在亚类中有明显区别的采用同类进行变化率计算。使用此方法得到的预测精度较传统方法有所提高。

3. 基于RUSBoost算法的玻璃制品鉴定

3.1. RUSBoost模型的建立

RUSBoost算法首先利用初始训练集训练出一个基学习器；再根据当前基学习器的训练误差调整每个训练样本的分布权重，通过增大被误分类样本的权重，使其在后续训练过程中获得更多关注；然后利用调整后的样本训练出下一个基学习器；如此反复迭代，直至生成T个基学习器；最后将根据上述T个基学习器的加权投票结果来预测未标记样本 [3] [4]。

假设训练数据集为 $S=\left\{\left(x_1,y_1\right),\left(x_2,y_2\right),\cdots ,\left(x_m,y_m\right)\right\}$，其中 $\left(x_i,y_i\right)$ 是一组样本数据， RUSBoost算法在第t次迭代中生成的基学习器为 $h_t$，每组样本 $X_i$ 对应的预测结果为 $h_t\left(x_i\right)$，并且该样本在该次迭代中对应的权重为 $D_t\left(i\right)$。因此，在第t次迭代中所有样本构成离散分布。

离散分布为 ${\Lambda }_t=\left\{D_t\left(1\right),D_t\left(2\right),\cdots ,D_t\left(m\right)\right\}$，其中 $Z_t={\displaystyle {\sum }_{i=1}^m{D}_t\left(i\right)}=1$。同时，每个基学习器 $h_t$ 的错误率定义为在分布 ${\Lambda }_t$ 下其预测错误的概率，即

${ϵ}_t=P_{x\sim {\Lambda }_t}\left[h_t\left(x\right)\ne y\right]$ (8)

其中， ${\alpha }_t$ 是基学习器 $h_t$ 的权重。而RUSBoost算法的预测结果表示为 $\text{sign}\left[H\left(x\right)\right]$。

RUSBoost算法是通过最小化指数损失函数 $l_{\exp }\left(H|\Lambda \right)$ 达到贝叶斯最优错误率：

$\begin{array}{c}{l}_{\exp }\left({\alpha }_t{h}_t|{\Lambda }_t\right)=E_{x~{\Lambda }_t}\left\{{\text{e}}^{-y{\alpha }_t{h}_t\left(x\right)}\right\}\\ ={\text{e}}^{-{\alpha }_t}{P}_{x~{\Lambda }_t}\left[h_t\left(x\right)=y\right]+{\text{e}}^{{\alpha }_t}{P}_{x~{\Lambda }_t}\left[h_t\left(x\right)\ne y\right]\\ ={\text{e}}^{-{\alpha }_t}\left(1-{ϵ}_t\right)+{\text{e}}^{{\alpha }_t}{ϵ}_t\end{array}$ (9)

令导数为零：

$\frac{\partial l_{\exp }\left({\alpha }_t{h}_t|{\Lambda }_t\right)}{\partial {\alpha }_t}=0$ (10)

则最优权重计算为：

${\alpha }_t=\frac{1}{2}\ln \left(\frac{1-{ϵ}_t}{{ϵ}_t}\right)$ (11)

RUSBoost算法在t + 1次迭代之后获得的基学习器线性组合为 $H_{t+1}=H_t\left(x\right)+a_{t+1}{h}_{t+1}$，其指数损失函数通过泰勒展开可以近似为：

$\begin{array}{c}{l}_{\exp }\left(H_{t+1}|\Lambda \right)=E_{x~\Lambda }\left\{{\text{e}}^{-y\left[H_t\left(x\right)+a_{t+1}{h}_{t+1}\left(x\right)\right]}\right\}\\ \cong E_{x~\Lambda }\left\{{\text{e}}^{-y{H}_t\left(x\right)}\left[1-y{\alpha }_{t+1}{h}_{t+1}\left(x\right)+\frac{1}{2}{y}^2{\alpha }_{t+1}^2{h}_{t+1}^2\left(x\right)\right]\right\}\\ =E_{x~\Lambda }\left\{{\text{e}}^{-y{H}_t\left(x\right)}\left[1-y{\alpha }_{t+1}{h}_{t+1}\left(x\right)+\frac{1}{2}{\alpha }_{t+1}^2\right]\right\}\end{array}$ (12)

因此，t + 1次迭代的最优基学习器为：

$\begin{array}{c}{h}_{t+1}^{*}\left(x\right)=\arg \underset{h_{t+1}}{\min }{l}_{\exp }\left(H_{t+1}|\Lambda \right)\\ =\arg \underset{h_{t+1}}{\min }{E}_{x\sim \Lambda }\left\{{\text{e}}^{-y{H}_t\left(x\right)}\left[1-y{\alpha }_{t+1}{h}_{t+1}\left(x\right)+\frac{1}{2}{\alpha }_{t+1}^2\right]\right\}\\ =\arg \underset{h_{t+1}}{\max }{E}_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_t\left(x\right)}y{\alpha }_{t+1}{h}_{t+1}\left(x\right)\right]\\ =\arg \underset{h_{t+1}}{\max }{E}_{x\sim \Lambda }\left\{\frac{{\text{e}}^{-y{H}_t\left(x\right)}}{E_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_t\left(x\right)}\right]}y{\alpha }_{t+1}{h}_{t+1}\left(x\right)\right\}\end{array}$ (13)

其中， $E_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_t\left(x\right)}\right]$ 是一个常数。令：

${\Lambda }_{t+1}\left(x\right)=\frac{\Lambda \left(x\right){\text{e}}^{-y{H}_t\left(x\right)}}{E_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_t\left(x\right)}\right]}$ (14)

${\Lambda }_{t+1}\left(x\right)$ 定义了一种新分布，故：

$h_{t+1}^{*}\left(x\right)=\arg \underset{h_{t+1}}{\max }{E}_{x\sim \Lambda }\left[y{\alpha }_{t+1}{h}_{t+1}\left(x\right)\right]$ (15)

则有：

$h_{t+1}^{*}\left(x\right)=\arg \underset{h_{t+1}}{\max }{E}_{x\sim \Lambda +1}\left\{\left[1-2I\left(h_{t+1}\left(x\right)\ne y\right)\right]{\alpha }_{t+1}\right\}=\arg \underset{h_{t+1}}{\min }{E}_{x\sim \Lambda +1}\left[I\left(h_{t+1}\left(x\right)\ne y\right)\right]$ (16)

所以， $h_{t+1}$ 是在分布 ${\Lambda }_{t+1}$ 下以最小化分类误差为优化目标而训练得到的基学习器。而分布 ${\Lambda }_{t+1}$ 可通过如下递推公式计算获得：

$\begin{array}{c}{\Lambda }_{t+1}\left(x\right)=\frac{\Lambda \left(x\right){\text{e}}^{-y{H}_t\left(x\right)}}{E_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_t\left(x\right)}\right]}=\frac{\Lambda \left(x\right){\text{e}}^{-y{H}_t\left(x\right)}{\text{e}}^{-y{\alpha }_t{H}_t\left(x\right)}}{E_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_t\left(x\right)}\right]}\\ ={\Lambda }_{t+1}\left(x\right){\text{e}}^{-y{\alpha }_t{H}_t\left(x\right)}\frac{E_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_{t-1}\left(x\right)}\right]}{E_{x\sim \Lambda }\left[{\text{e}}^{-y{H}_t\left(x\right)}\right]}\end{array}$ (17)

3.2. 模型求解

将表单2中的空缺化学成分进行填零处理，将玻璃类型为高钾和铅钡两种类型。将表单2中的数据带入Matlab工具箱中的RUSBoost模型中，通过混淆矩阵得知，训练集和测试集的平均准确度可以达到100%

由图4中的ROC曲线可知，AUC = 1.00 > 0.95，说明分类预测结果十分可靠，可以进行下一步求解。

Figure 4. Confusion matrix and ROC curve

图4. 混淆矩阵和ROC曲线图

表3中各未知玻璃类型的文物预测结果如下表5：

Table 5. Prediction result table

表5. 预测结果表

4. 结论

本文根据已有数据建立小样本数学模型，能够准确的预测出玻璃文物风化前后的化学成分变化，并且能够对玻璃文物类别进行亚类划分以及鉴别。在预测风化成分方面，相较于传统方法，采用先分类后预测的方法，能够大幅提高预测准确性。在亚类划分方面，采用K-Means和BP神经网络结合的方法，能够在兼顾分类多样性和准确性下进行亚类划分。在鉴别玻璃类型方面，对小样本数据采用RUBoost机械学习方法，可以准确的鉴别玻璃类型。综上，本篇论文通过建立数学模型对鉴别古代玻璃类型，研究有无风化玻璃制品成分的变化规律，以及高钾、铅钡两种玻璃类型的化学成分统计规律和探索亚类方法方面做出了合理解答。