北京工商大学数学与统计学院，北京

收稿日期：2023年6月14日；录用日期：2023年7月9日；发布日期：2023年7月19日

摘要

基于已知玻璃类型的数据，研究玻璃文物表面风化的相关因素和由风化引起的主要化学成分的变化，由此基于主要化学成分建立玻璃的分类模型和亚类划分模型，并利用所得模型对未知类型的玻璃文物进行鉴别。研究发现：玻璃表面风化与其类型相关，与纹饰弱相关，而与颜色基本无关；与高钾玻璃相比，铅钡玻璃表面更易受风化影响；表面风化会引起高钾和铅钡玻璃一些成分比例的显著变化；精选4个重要成分PbO、BaO、SiO₂、K₂O，建立玻璃类型的随机森林(RF)和支持向量机(SVM)模型，分类准确率均为100%，使用这些模型对8个未知类型的玻璃文物进行鉴别，结果一致，均为：A1 A6 A7为高钾，其余为铅钡；对每个类型的玻璃各挑选8个重要成分，使用k-means聚类算法对每个类型玻璃可以很好地细分为两个亚类；最后，通过对成分数据增加微小随机扰动来证实所建模型是稳健的。

关键词

随机森林，支持向量机，K-Means聚类算法，成分数据

Classification Identification and Sub-Classification of Glass Relics

Meiping Xu^*, Xingxiu Wang^*, Mengmeng Jang, Siyu Qiao, Yunchuan You

School of Mathematics and Statistics, Beijing Technology and Business University, Beijing

Received: Jun. 14^th, 2023; accepted: Jul. 9^th, 2023; published: Jul. 19^th, 2023

ABSTRACT

Based on glass relics data with known glass types, factors and important chemical compositions related to the weathering of glass surface are studied. And classification models of glass types and sub-classification models are established based on their important chemical compositions. Then these classification models are applied to identify glass relics with unknown types. Results show that the weathering of glass surface is related to its type, weakly related to its decoration, and unrelated to its color. The surface of PbO-BaO glass is more susceptible to weathering than that of high K₂O glass. Surface weathering can cause significant changes in the proportion of some compositions of high K₂O and PbO-BaO glass. Four important compositions, PbO, BaO, SiO₂ and K₂O, are selected to make models of glass types via random forest (RF) and support vector machine (SVM), with classification accuracy of 100%. These models are employed to identify 8 unknown types of glass relics, same results are presented that A1, A6 and A7 are high K₂O glass, and the rest are PbO-BaO glass. Then 8 important compositions are selected for each glass type, both of which are divided into two subclasses via k-means clustering very well. Finally, the robustness of the proposed models is verified by adding a small random disturbance to the composition data.

Keywords:Random Forest, Support Vector Machine, K-Means Clustering Algorithm, Composition Data

Copyright © 2023 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

玻璃制品最初由丝绸之路传入中国，主要类别有高钾、钠钙玻璃等。高钾玻璃是以含钾量高的物质如草木灰作为助熔剂烧制而成的，主要流行于我国岭南以及东南亚和印度等区域。在中国本土发明制作的玻璃中，常用铅矿石作助熔剂，出土的玻璃文物往往铅钡含量较高。留存至今的古代玻璃制品往往长时间埋藏于墓穴等环境中，极易受环境的影响而风化。在风化过程中，内部元素与环境元素进行大量交换，导致其化学成分所占比例发生变化，从而影响对其类别的正确判断。另外，由于被埋藏环境和时间等条件的不同，不同类别玻璃文物风化的程度也有所差异，也给分类问题增加较大辨别难度。

关于风化，周良知 [1] 应用近代研究方法，研究影响硅酸盐玻璃风化的主要因素。王承遇等 [2] 对不同类型玻璃风化作用的机理进行了深入研究。王婕等 [3] 利用分析测试方法对战国时期八棱柱状铅钡玻璃器的风化产物、风化形貌和元素迁移进行了分析和探讨。这些工作对后续研究风化玻璃的主要成分提供了依据和借鉴。

关于分类，付强等 [4] 通过探讨风化对化学成分定量分析的影响，对所占比例较高的硅钾酸盐玻璃进行了亚类划分。赵凤燕等 [5] 通过pXRF方法对玻璃器的化学成分进行无损分析来进行分类。但他们仅基于成分数据做了一些简单的分析就依此给出类型识别，没有借助统计模型做更深入的探讨和分析，提供一些可以按照玻璃成分进行较为精准分类的泛化模型。近期，王祉皓等 [6] 基于机器学习的方法对风化硅酸盐原成分进行了预测及亚类划分，但其主要使用的神经网络算法对文中所采用的小样本数据很容易过拟合，容易导致模型的泛化能力和稳健性变差。

针对本文数据特点，我们借鉴何秋妍等 [7] 使用的有监督模型，选用随机森林(RF)和支持向量机(SVM)这两个既适用于小样本又有较好分类能力的模型对玻璃文物按高钾、铅钡两大类进行分类规律的讨论，使用k-means这一最简单且收敛速度较快的硬聚类算法对高钾/铅钡玻璃进行亚类划分，并对所建模型做稳健性分析以说明模型的可适用性，最后探讨表面风化对分类的影响。研究结果表明我们的模型在保证分类准确类的前提下，更简洁稳定，且具有很好的泛化能力。

2. 思路流程

如图1所示，本文首先对于所获取成分数据使用CLR (对数中心比)变换，然后基于转换后的数据进行独立性检验与单因素方差分析，选取最显著的8个变量，并使用随机森林中算法进行分类，同时使用具有线性和径向基核函数的SVM算法进行回判分析，佐证稳健性，最终得到高钾玻璃与铅钡两类玻璃的预测模型。然后使用K-means与多维尺度分析(cmd)对处理后的数据进行降维处理与亚类划分，将分类结果分别命名为铅钡(髙硅低铅)类玻璃、铅钡(髙铅低硅)类玻璃、高钾(髙硅)类玻璃与高钾(低硅)类玻璃。

图1. 研究思路图

3. 模型建立

3.1. 数据来源与预处理

通过查阅资料 [8] ，获得玻璃类型、纹饰、颜色与14个主要成分比例的玻璃文物相关数据。具体为：表单1中列出的玻璃文物已经知道其玻璃类型、颜色、纹饰和表面风化与否，表单2是表单1中文物在一些采样点的成分数据，表单3给出的是未分类玻璃文物的成分数据。表单2和3的数据具有成分性，但考虑到受限于当前的检测手段，可能存在检测得到成分比例和不等于100%的情况，项目方认为成分比例累加和位于85%至105%之间的数据有效，其余数据视作无效。

对表单2中编号18和20的两行无效数据进行删除，共获得67行有效观测数据；对表单2和3的未检测到的成分数据填补一个很小的正数1e−3，方便后续做几何平均；对由于采集误差等原因造成的14个主要成分数据占比和不为100%的问题进行归一化处理。由于这些成分数据取值范围是14维空间中的一个单纯形，我们通过式(1)的中心对数比变换(CLR)将其映射到14维实数空间。

$y_{ij}=\ln \frac{x_{ij}}{g\left(x_i\right)}$ , $g\left(x_i\right)={\left({\displaystyle \underset{j=1}{\overset{d}{\prod }}{x}_{ij}}\right)}^{1/d}$ (1)

其中， $x_{ij}\left(i=1,\cdots ,n;j=1,\cdots ,d\right)$ 表示第j列成分数据的第i个观测，n和d分别为有效观测个数和成分个数。由于这些主要成分之间具有完全多重共线性，为此在对训练数据建模前需要对其进行特征筛选，选出其中对区分不同类别玻璃的重要成分，方便后续分类模型的建立和亚类划分。本文所有计算及图形均使用R软件完成。

3.2. 风化与类型、纹饰和颜色的独立性检验

通过对玻璃的表面风化状态与玻璃类型、纹饰和颜色绘制并列条形图(图2)进行可视化分析，观察其差异。我们发现：与高钾玻璃相比，铅钡玻璃表面更易受风化影响；与其它纹饰相比，好像纹饰B易受风化影响；而风化状态和颜色的相关性从图2看不出。考虑到样本量较小，采用χ²独立性检验分别对风化状态与类型、纹饰和颜色的3个二维列联表结果进行检验，给出相应的χ²值和相伴p值。考虑到风化与纹饰、颜色的列联表中出现观测频数小于5的情况，使用卡方检验容易出现偏差，还给出了Fisher精确检验的模拟p值(表1)。针对颜色分组太细的情况，把频数小于5的组按相近色合并为4个组后再次绘制其并列条形图并作上述检验，相应结果在图2(d)和表1的最后一行给出。

图2. 表面风化状态和玻璃类型、纹饰和颜色绘制并列条形图

表1. 卡方和Fisher精确检验结果

注：表1中^*和^#分别表示在0.05和0.1的水平下显著。

从表1的检验结果可得到结论：基于现有数据，认为玻璃表面风化与其类型不独立，与纹饰弱相关，而与颜色基本无关，与图2的可视化分析结果一致。

3.3. 成分含量与风化的统计规律

对高钾/铅钡玻璃经过CLR变换后的数据绘制其各成分数据在不同风化状态下的均值的并列条形图进行可视化分析，图3呈现的是其中有显著变化的成分的均值对比。

从图3可发现风化后高钾玻璃成分比例显著增加的有SiO₂、CuO、CaO、Fe₂O₃、Al₂O₃，显著下降的有K₂O、MgO、PbO、BaO；铅钡玻璃成分比例显著下降的有SiO₂、Na₂O、K₂O、Al₂O₃ ，显著增加的有P₂O₅、CaO、Fe₂O₃、SrO。为更好地评价上述有显著变化的成分对表面风化的显著性，我们也对它们做了箱型图及单因素方差分析进行比较，表2列出了相应的F值和相伴p值。需要说明的是，此处及下文SnO₂和SO₂的显著性我们均未予以考虑，原因是接近90%的采样点未检测出这两种成分。

图3. 高钾(上)/铅钡(下)玻璃在不同风化状态下部分成分的均值比较

表2. 高钾/铅钡玻璃部分成分对表面风化的方差分析

注：^***、^**、^*和^#分别表示在0.001、0.01、0.05和0.1水平下显著。

从表2可看出：风化会引起高钾玻璃SiO₂成分比例的显著变化，一定程度上引起CuO、Al₂O₃成分比例的变化；而铅钡玻璃由风化引起显著变化的成分比例较多，按显著性由强到弱排序，依次为SiO₂、P₂O₅、Al₂O₃、Na₂O、PbO、K₂O、CaO、SrO。这些结果与图3的可视化分析结论是相互印证的。

3.4. 分类模型的建立

为识别能够明显区分两类玻璃的重要成分，对14个CLR变换后的成分数据绘制了它们对类型的箱型图进行比较，并做各成分对类型的单因素方差分析，从中挑选出8个有显著差异的成分SiO₂、K₂O、Al₂O₃、Fe₂O₃、CuO、PbO、BaO、SrO (图4和表3)。

图4. 重要成分对玻璃类型的箱线图比较

表3. 重要成分对玻璃类型的方差分析

注：表3中显著性标注同表2。

使用表3中所列的8个重要成分，建立玻璃类型对它们的归一化数据的RF模型 [9] 。每次有放回地抽n个观测，建了500棵树，每个节点使用2个变量，算得2个特征重要性筛选的指标值：OOB数据作为测试集所做的交叉验证得到的平均减少精度(图5左)和平均减少Gini值(图5右)，从中精选出4个足以区分两类玻璃的重要成分：PbO、BaO、SiO₂、K₂O。使用此4个重要成分建立对玻璃分类的RF模型，误判率的OOB估计值为0%，分类准确率为100%。不考虑风化影响，直接运用该模型对表单3中8个未知类别的玻璃文物进行鉴别，给出其所属类型，结果见表4。

图5. 平均减少精度和平均减少Gini值下的重要特征排序

表4. 未知类型的玻璃文物进行鉴别结果

再分别应用具有线性和径向基核函数的SVM建立玻璃类型对上述4个重要成分的分类模型，分类准确率也是100%，对表单3中未知类别的玻璃文物进行鉴别，给出的结果与表4一致。这说明RF和SVM模型对本文数据的识别度很高，且只需通过4个重要成分(PbO、BaO、SiO₂、K₂O)的占比就可以对玻璃文物类别做出可靠的判别。

实际上，仔细观察表单3中被鉴别为高钾玻璃的编号为A1，A6，A7的化学成分，会发现它们不论风化与否，都有明显高于铅钡玻璃的SiO₂占比，且都未检测到PbO、BaO成分。事实上，较高的SiO₂占比也是高钾玻璃区别于铅钡玻璃的一个显著特点，而铅钡玻璃的PbO、BaO占比要明显高一些，这一点从2.5中图6容易看出。

另外，2.3的研究表明表面风化会对玻璃主要化学成分的占比有一定影响，但在筛选重要变量环节，我们发现表面风化对玻璃类型的影响排在上述4个重要成分之后，因而使用我们的模型鉴别玻璃类型无需特别关注表面风化因素，但这并不意味着它对鉴别玻璃类型不起作用。事实上，由于风化会引起SiO₂、Al₂O₃在高钾玻璃中明显富集，而在铅钡玻璃中显著贫乏；PbO在高钾玻璃中则显著减少、在铅钡玻璃中显著增加，从而使得玻璃类型的鉴别更加容易。

3.5. 亚类划分

亚类划分是一个无监督分类问题，同样由于14个化学成分的相关性，在使用聚类模型前，需要先对每个类型的玻璃选择合适的成分。通过比较高钾/铅钡玻璃各成分均值大小的方法给出比例由高到低的前8个重要成分的排序结果：高钾玻璃为SiO₂、K₂O、Al₂O₃、CaO、CuO、Fe₂O₃、P₂O₅、MgO，铅钡玻璃为SiO₂、PbO、BaO、Al₂O₃、P₂O₅、CaO、CuO、Na₂O (图6)。

图6. 高钾/铅钡玻璃前8个重要成分的均值排序

使用k-means算法对高钾/铅钡进行亚类划分 [9] 。使用类内总平方和(WSS)准则，高钾、铅钡各分两类时，WSS的下降基本平缓，其占总平方和的比例分别达到80.6%和61.9%。考虑到高钾、铅钡玻璃分别只有18和49个采样点，样本量偏小，我们只对各个类型玻璃再细分为两个亚类。仍以高钾玻璃为例阐述具体聚类过程：先对所选出的8个重要成分比例数据的每一列 $x_j$ 使用0-1正规化方法

$\left[x_j-\min \left(x_j\right)\right]/\left[\max \left(x_j\right)-\min \left(x_j\right)\right]$ 做同一尺度化处理；再利用高维空间中的欧氏距离算出数据的相似

性矩阵；然后使用cmd进行数据降维，图7呈现的是降到二维的相应分类散点图，其中同一形状的点表示在同一类。

从图7 (左)可看出，从cmd的第一个维度判断是否小于0就可以把高钾玻璃清晰地分为两类，通过对它们各自所含重要成分的均值比较，发现SiO₂的含量足以区分这两类，因而按照SiO₂含量高或低依次命名为高钾(髙硅)类和高钾(低硅)类，分别对应图7 (左)左侧的8个点和右侧的10个点，详见表5。

对铅钡玻璃，从图7 (右)可看出，从cmd的第一个维度判断是否大于0.1可以把它们清晰地分为两类，再通过对它们各自所含重要成分的均值比较，发现右侧21个点的SiO₂含量偏高，PbO含量偏低，而左侧28个点的PbO含量高，SiO₂含量偏低，故将它们分别命名为铅钡(髙硅低铅)类和铅钡(髙铅低硅)类，详见表5。

图7. 高钾/铅钡玻璃的k-means聚类结果

表5. 高钾/铅钡玻璃的亚类划分结果

3.6. 模型的稳健分析

为考察2.4所建分类模型和2.5亚类划分模型的稳健性，我们通过对数据增加微小随机扰动对模型的稳健性进行分析。具体为：在预处理后的表2和表3的成分比例数据的每个数值上，加上一个均值为0的任意正态随机数，标准差依次取1e−4，1e−3，1e−2，1e−1，以表示对数据的干扰程度逐渐加强，然后对扰动后的数据重复2.4和2.5的工作共4次，结果均显示2.4所给几个模型的分类准确率、对表单3的分类鉴别结果和2.5的亚类划分结果没有改变，这说明我们的模型是稳健的，在后续玻璃文物的鉴别和亚类划分中均可放心使用。

4. 补充说明

我们也根据风化点检测成分数据，先预测其风化前的成分含量。以高钾玻璃为例，对由风化引起比例有显著变化的成分，先通过QQ图和KS检验检查其风化与未风化组各自中心尺度化后的分布匹配情况，此处选择中位数和MAD分别作尺度化的中心和尺度，主要是因为它们较均值和标准差更稳健，分布的匹配程度也更好。我们发现大多数成分数据的分布是匹配的，对风化成分数据利用式(2)的配对方法预测其无风化的数值。

${\hat{y}}_{ij}=m\left(Y_{nofen,j}\right)+\frac{s\left(Y_{nofen,j}\right)}{s\left(Y_{fen,j}\right)}\left[y_{ij}-m\left(Y_{fen,j}\right)\right]$ (2)

其中 $Y_{fen,j}={\left(y_{1j},\cdots ,y_{n_1,j}\right)}^{\text{T}}$ ， $n_1$ ， $m\left(Y_{fen,j}\right)$ ， $s\left(Y_{fen,j}\right)$ 分别表示风化玻璃的第j个成分向量及其样本个数、中位数和MAD，未风化玻璃的相应记号为 $Y_{nofen,j}$ ， $m\left(Y_{nofen,j}\right)$ ， $s\left(Y_{nofen,j}\right)$ ， $j=1,\cdots ,d$ 。铅钡玻璃的分析

和预测类似，故略去。

接着，对高钾/铅钡玻璃的风化成分数据使用其预测数据，通过CLR的逆变换将其还原为比例数据后与未风化成分比例数据合并，作为2.4和2.5的建模数据进行分析，发现对于那些由检测误差造成的成分比例之和明显偏低的观测值，由于样本观测太少及匹配变换产生的估计误差使得成分比例估计偏差较大，从而造成模型误判。同样采用4个成分PbO、BaO、CaO、Fe₂O₃，建立对玻璃分类的RF模型，误判率的OOB估计值为4.48%，有3个观测被误判，运用该模型对表单3的玻璃文物进行鉴别的结果同表4，SVM模型所给结论类似，但亚分类的结果与2.5有较大差异。这意味着在检测和分析的各环节尽量减少误差对玻璃文物类别的准确识别是非常必要的。

5. 结束语

本文基于玻璃文物的主要化学成分，应用RF和SVM这两种适用于小样本的分类模型对其进行类型鉴定，给出了仅使用4个重要成分、分类准确率100%的简约模型，并对未知类型的玻璃给出分类结果；使用k-means聚类算法对两类玻璃进行亚类划分，分别基于8个重要成分对所有样本给出了清晰的划分；通过对成分数据增加微小随机扰动证实了所建模型的稳健性；最后，通过补充说明对另一种可行的方法阐述了我们的观点。对于本文这种小样本数据，RF和SVM的泛化能力和稳健性一般优于神经网络等复杂机器学习模型，这是我们采用它们的重要原因。由于文物本身具有稀缺性的特点，考古学与其它学科的交叉研究明显偏少，本文把统计模型这一便捷工具应用于文物鉴定方面，旨在抛砖引玉，希望可以在一定程度上填补这方面的不足。本文的研究方法也适用于其他行业的样本分类问题。

基金项目

2023年北京市大学生创新创业训练计划(B008)。

文章引用

徐美萍,王行修,江萌萌,乔斯昱,由云川. 玻璃文物的类别鉴定与亚类划分
Classification Identification and Sub-Classification of Glass Relics[J]. 应用数学进展, 2023, 12(07): 3188-3199. https://doi.org/10.12677/AAM.2023.127319

参考文献