出土玻璃制品成分分析与鉴别研究

杜世强,乐靖雯,周海丽,杨松杰

(西北民族大学 数学与计算机科学学院,甘肃 兰州 730030)

0 引言

古代玻璃制品是考证丝绸之路历史发展变化的重要物品,而研究它的埋藏环境与风化之间的关系,对于充分挖掘中国古代对外商贸的历史渊源具有重要意义.文献[1]采用多元统计方法对四川、重庆、贵州、广西、广东地区出土的100余件古玻璃样品化学成分数据进行分析,以期能够提示我国汉代南方和西南地区的玻璃生产情况.文献[2]对广州出土的西汉早期至东汉的46件玻璃单色珠、耳珰等器物的化学成分进行了成分分析,探讨表面风化对化学成分定量分析的影响,并依据CaO、Al2O3含量,结合微量元素Rb和Sr的比例,对所占比例最高的钾硅酸盐玻璃进行了亚类划分.文献[3]对20多种较有代表性的中国古代铅玻璃样品进行实物考察和分析研究,根据化学检测结果,这批古玻璃可分为五个成分系统,并可归纳为三类:铅钡玻璃类有PbO-BaO-SiO2和Na2O-PbO-BaO-SiO2,铅硅玻璃类有PbO-SiO2,碱铅硅玻璃类有K2O-PbO-SiO2和Na2O-PbO-SiO2.文献[4]利用EDXRF探针无损分析技术,分析了河南禹县阳翟遗址出土的一批金元时期的玻璃样品,结果表明,这批古玻璃种类较为复杂,主要包括钾钙玻璃、钠钙玻璃和铅钾钙玻璃等.文献[1-4]均采用化学方法对古玻璃进行性能测定,其中文献[2-3]主要通过主观判断来划分古玻璃类别,文献[1,4]则采用聚类分析和因子分析验证了分类结果.如果对大批量样品的组成成分进行分类时,单凭主观判断易造成误差,且缺乏说服力,所以本文在参考文献[1-4]分类结果的基础上,利用多种统计分析方法进行数据分析,从而得出较为合理的对古玻璃分类的标准.

利用卡方检验和相关性分析玻璃文物表面风化和类型、纹饰、颜色的关系以及所含化学成分之间的相关性.结合文献[4],对化学成分的主要用途进行分类,再根据主要化学成分建立层次聚类和K-means聚类模型,得出玻璃文物分类、亚分类标准,最后构建随机森林分类器模型,进而对未知类型玻璃进行分类.

1 相关因素分析

1.1 数据预处理

分析2022年高教社杯全国大学生数学建模竞赛C题《古代玻璃制品的成分分析与鉴别》数据,发现玻璃制品颜色列存在缺失值,检测的化学成分存在较多空值.依据2022年高教社杯全国大学生数学建模竞赛C题所给背景条件,将化学成分比例累加和介于85%～105%之间的数据视为有效数据,通过计算,编号15和编号17为无效数据,把这两个编号从数据表中剔除,对于未检测到成分而缺失的空值,用“0”填充.对颜色列缺失,采用热卡填充方法对含有缺失值元组的其他属性与其他完整数据元组的属性求欧氏距离,将距离最短的元组中对应的值作为缺失值的估计量进行填充,计算公式如式(1).

其中,L表示欧氏距离,xi表示待填充的数据属性,yi表示完整的数据属性.经计算,“颜色”列中的空白处依次填充为“黑”“紫”“浅蓝”“浅绿”.

1.2 玻璃文物表面风化与其类型、纹饰和颜色的相关性分析

其中,A表示实际频数,T表示理论频数,且:

其中,R表示行,C表示列,TRC表示第R行和第C列的理论数,nR,nC表示第R行和第C列的合计数,N表示总的合计数.最后可根据χ2值查阅卡方界值表,得到P值［5］.

以风化情况作为变量X,以颜色、纹饰和类型作为变量Y,利用python编程实现上述模型,得卡方检验表1.

表1 表面风化与纹饰、玻璃类型、颜色相关性卡方检验表

从表1可看出,玻璃风化情况与纹饰、颜色不存在显著差异,但与类型存在显著差异.

1.3 玻璃文物中化学成分间的相关性分析模型建立与求解

由协方差方程Cor(X,Y)=E{[X-E(X)][Y-E(Y)]}可知其相关系数,14种化学成分的相关性分析模型如式(4):

其中,i=1,2,…,14,j=1,2,…,14.

对模型进行求解,得出不同类别玻璃文物化学成分之间的相关系数,如图1所示.

高钾玻璃中氧化锶与五氧化二磷间相关系数为0.78,呈显著正相关;氧化铁和五氧化二磷间相关系数为0.81,呈显著正相关;二氧化硅和氧化钾间相关系数为-0.86,呈显著负相关.而氧化锶和氧化钙,五氧化二磷和氧化铅、氧化锡,氧化钠和氧化铜、氧化铁,相关性最低,独立性最强.

铅钡玻璃中氧化铅与氧化锶间相关系数为0.63,呈显著正相关;二氧化硫与氧化铅间相关系数为-0.85,呈显著负相关.氧化铜和氧化铁,氧化镁和二氧化硫,氧化钠和氧化锶相关性最低,独立性最强.

图1 高钾玻璃和铅钡玻璃化学成分相关系数图

2 玻璃文物中化学成分预测

文物在烧制过程中使用的原材料和助熔剂会直接与空气发生快速氧化还原反应,使得文物表面本身的钾元素或铅元素含量在最开始时就达到最大值.文物经过长时间埋藏,表面化学元素与周围环境所含化学元素进行大量交换,最终文物和埋藏环境的化学元素含量达到平均水平,即出土时间较长文物表面化学元素含量占比少于最开始出土时的化学元素含量,但仍具有可识别性.

铅钡玻璃前三类占比最高的化学成分为:氧化铅(PbO)、氧化钡(BaO)、氧化锶(SrO);高钾玻璃前两类占比最高的化学成分为:氧化钾(K2O)、氧化铁(Fe2O3).高钾玻璃预测数据部分选用这两个强相关性化学元素,铅钡玻璃预测数据部分选用这三个强相关性化学元素.

2.1 风化前化学成分预测模型的建立与求解

基于变量设定,构建一个二项逻辑回归模型,对影响文物类别划分的因素进行计量分析,规定因变量为0-1型二分变量,其中“1”表示为无风化,“0”表示为风化,解释变量为所给的化学成分含量数据类别.

Logit模型的具体形式如式(5)所示:

Logit(p)=b0+b0x0+…+bpxp

(5)

根据Logit变换定义,有:

将式(6)进行运算,可得:

其中b0为常数项,b1,b2,…,bp为偏回归系数.

由于因变量属于二分变量,因此运用极大似然法进行估计.文中使用了Logit模型变换,各自变量的偏回归系数bi(i=1,2,…,p)表示自变量xi每改变一个单位,类别划分为高钾玻璃和划分为铅钡玻璃的发生比自然对数值变化量.

利用所给数据采用极大似然函数进行结果溯因分析,估计出使得目前结果可能性最大的参数.根据该参数能够确定任意未知类别样品的划分概率和后验概率,最后得到极大似然估计函数［6］为:

根据玻璃类型分组,使用SPSS进行Logit回归,回归准确度如表2所示.

表2 Logit回归准确度表

由表2可知,回归总体准确度为84.6%,预测精度较高,因此该模型可用于预测风化点风化前的强相关化学成分含量.

通过该模型得到强相关化学成分含量预测结果,利用SPSS画折线图分析两类玻璃弱相关化学成分与强相关化学成分之间的关系,总结如下规律:

1)高钾玻璃和铅钡玻璃:氧化钙(CaO)与二氧化硅(SiO2)、氧化镁(MgO)与氧化铁(Fe2O3)存在某种函数关系,剩余弱相关化学成分与氧化钾(K2O)存在某种函数关系.

2) 两种玻璃:氧化钾(K2O)、氧化铜(CuO)、二氧化硫(SO2)与二氧化硅(SiO2)存在某种函数关系,氧化铝(Al2O3)与氧化铅(PbO)存在某种函数关系,氧化镁(MgO)、氧化钙(CaO)、氧化铁(Fe2O3)与氧化钡(BaO)存在某种函数关系.

以高钾玻璃中氧化钙(CaO)与二氧化硅(SiO2)为例,建立式(9)所示关系式:

y=-0.002901x3+0.5583x2-35.54x+754.3

(9)

利用python拟合后函数图像如图2所示,由预测所得风化前二氧化硅含量,进而求得氧化钙含量.铅钡玻璃文物类型风化前化学成分含量也可求得.

图2 高钾玻璃文物氧化钙(CaO)与二氧化硅(SiO2)拟合图像

3 玻璃文物的分类和亚分类规则

3.1 玻璃文物分类规则选择

对多次取样的文物,计算其化学成分的平均值,用于表征文物的主要化学成分.查阅相关资料［4］,将各化学成分的主要用途进行分类,如表3所示.

通过表3数据,将着色剂和含量过低的其他氧化物剔除,判断文物化学成分中二氧化硅、氧化铝、氧化钾、氧化钙、氧化铅和氧化钡6种氧化物为影响古玻璃文物分类的主要因素.

将所给数据分成风化和未风化两部分,去除其余8项特征,利用python软件通过层次聚类［7］得到如图3所示结果.表面风化的文物除编号2和编号48外,其他分类均正确,表面未风化文物除编号21外,其他分类均正确,该结果准确度较高.

图3 风化与未风化聚类结果图

从聚类分析树形图可知,如从阀值λ约为0.7处 (图3 b)划分,可将所有样品划分为两大类.

1)古玻璃样品属于高钾玻璃:PbO、BaO含量较低,K2O、CaO含量较高.K2O是主要助熔剂氧化物,CaO可能是石灰石稳定剂生成的产物或是其他助熔剂氧化物.

2)古玻璃样品属于铅钡玻璃:主要特点是PbO、BaO含量很高,K2O、CaO含量较低.PbO、BaO是主要助熔剂氧化物.

上述两类玻璃的Al2O3都较高,但烧制过程中形成的中间氧化物不能作为分类依据.

结合层次聚类结果,观察可知两种玻璃的分类规律.

1)高钾玻璃分类规律:当文物表面氧化钾成分占比最高时可分为高钾玻璃类.

2)铅钡玻璃分类规律:当文物表面氧化铅和氧化钡成分占比最高时可分为铅钡玻璃.

3.2 玻璃文物亚分类规则选择

采用K-means算法［8］分别对高钾玻璃和铅钡玻璃进行亚分类聚类,具体算法如下.

Step1:在样本数据集D中选择k个样本点,将k个样本点的值分别赋值给初始聚类中心:

Step5:计算数据集D中所有点的平方差Ei,并且与前一次误差Ei-1作比较:

判断|Ei-1-Ei|<δ是否成立,若成立则算法结束,否则转入Step2进行二次迭代,直到|Ei-1-Ei|<δ,成立.

高钾玻璃选取氧化钠(Na2O)、氧化钙(CaO)、氧化镁(MgO)、氧化铝(Al2O3)和五氧化二磷(P2O5)成分进行分析,二氧化硅(SiO2)、氧化钾(K2O)是高钾玻璃的主要成分,氧化铁(Fe2O3)、氧化铜(CuO)是着色剂,氧化铅(PbO)、氧化钡(BaO)、氧化锶(SrO)、氧化锡(SnO2)、二氧化硫(SO2)等含量极低,可忽略不计.

铅钡玻璃选取氧化钙(CaO)、氧化铝(Al2O3)成分进行分析(此处删去了另外三种化学成分,更有利于类别命名).

利用python求解模型,结合手肘法和轮廓系数法可知,高钾玻璃分成四类较合理,铅钡玻璃分成三类较合理.

高钾玻璃聚类结果如表4所示,根据高钾玻璃中CaO 和 Al2O3的浓度以及是否含钠元素［2］,分别将0、1、2、3命名为含钠中铝中钙高钾玻璃、低钙低铝高钾玻璃、低钙高铝高钾玻璃、不含钠中铝中钙高钾玻璃.

表4 高钾玻璃分类结果

铅钡玻璃聚类结果如表5所示.根据铅钡玻璃中 CaO 和 Al2O3浓度,将0、1、2分别命名为低钙低铝铅钡玻璃、高钙中铝铅钡玻璃、中钙高铝铅钡玻璃.

表5 铅钡玻璃分类结果

3.3 合理性分析

表6 高钾方差分析

表7 铅钡方差分析

从表6、表7可以看出,不同类别样本选取的化学成分全部呈现出显著性(P<0.05),方差分析对参与聚类分析的变量能很好地区分类别,类间差异足够大.

3.4 玻璃文物分类模型和亚分类模型的建立与求解

随机森林预测算法流程:基于bootsrap方法重采样,随机生成N个训练集,产生的训练集与原始训练集P样本总数相等.利用训练集构建相对应的决策树K1,K2,…,KN.在每一个内部节点选择分裂特征前,从训练集中的M个特征中随机选取m(m≤M)个作为当前内部节点的分裂特征集,选择对应基尼值最小的特征作为分裂特征.每棵树不进行剪枝操作,都生长到底.对于测试集数据X,每个决策树给出独立的预测结果K1(X),K2(X),…,KN(X),即投出自己的一票.统计每一个决策树的预测结果,将所得票数最多的预测值作为最终的预测结果［9］.

根据制定的分类规则,构建两个随机森林分类器,将待分类数据分成风化和未风化,并选取相应化学成分,分类结果如表8.

表8 分类结果(部分)

4 结语

不同文物因组成材料不同、埋藏地点不同,受到的侵蚀风化程度也不相同.为了较好地修复文物,对古代玻璃制品成分进行分析,可依靠数学分析工具分析出文物表面风化情况与客观因素之间的关系.文化风化后,主成分含量均有所降低.根据检测成分对未知玻璃制品进行鉴别,结合k-means聚类算法建立随机森林分类器,推导出文物本身具有的化学元素,对文物修复保护着指导作用,同时也对研究我国古代玻璃制品的风化条件与玻璃制品表面特征和埋藏环境之间的关系具有积极意义.

本文研究成果获全国大学生数学建模大赛甘肃赛区一等奖.