基于成分数据的古代玻璃制品分析与分类.pdf

第1 2卷 第2期2023年6月数学建模及其应用M a t h e m a t i c a l M o d e l i n g a n d I t s A p p l i c a t i o n sV o l.1 2 N o.2J u n.2 0 2 3教学与竞赛基于成分数据的古代玻璃制品分析与分类马佩莹,韩雁来,李德兰,陈佳佳(山西财经大学 统计学院,山西 太原 0 3 0 0 0 6)摘 要:古代玻璃制品的化学成分属于成分数据,基于成分数据分析方法可以对玻璃制品的化学成分进行分析,研究其分类规律,并对未知玻璃文物鉴别其所属类型.首先,基于S p e a r m a n相关系数以及卡方检验分析玻璃文物表面风化与其类型、纹饰、颜色的关系;通过单形空间均值来分析玻璃表面有无风化化学成分含量的统计规律;构建D i r i c h l e t回归模型来预测风化点风化前的化学成分含量.其次,构建决策树、偏最小二乘判别分析两种模型对两类玻璃进行初分类特征选择;进一步,用K-m e a n s聚类对两类玻璃进行亚分类,并通过偏最小二乘判别分析对两类玻璃进行亚分类特征选择;进而,利用所得分类规律对未知类别玻璃鉴别其所属类型.最后,运用灰色关联分析分别探究两类玻璃化学成分之间的关联关系及其差异性.关键词:成分数据;决策树;偏最小二乘判别分析;K-m e a n s聚类;灰色关联中图分类号:O 2 9 文献标志码:A 文章编号:2 0 9 5-3 0 7 0(2 0 2 3)0 2-0 0 6 3-1 1 D O I:1 0.1 9 9 4 3/j.2 0 9 5-3 0 7 0.j mm i a.2 0 2 3.0 2.0 6收稿日期:2 0 2 3-0 1-1 9基金项目:山西省高等学校教学改革创新项目(J 2 0 2 2 0 5 7 0)通讯作者:陈佳佳,E-m a i l:c h e n j i a j i a 0 4 0 11 6 3.c o m引用格式:马佩莹,韩雁来,李德兰,等.基于成分数据的古代玻璃制品分析与分类J.数学建模及其应用,2 0 2 3,1 2(2):6 3-7 3.MA P Y,HAN Y L,L I D L,e t a l.A n a l y s i s a n d c l a s s i f i c a t i o n o f a n c i e n t g l a s s p r o d u c t s b a s e d o n c o m p o s i t i o n a l d a t a(i n C h i n e s e)J.M a t h e m a t i c a l M o d e l i n g a n d I t s A p p l i c a t i o n s,2 0 2 3,1 2(2):6 3-7 3.0 引言玻璃的发展历史悠久,是人类最早发明的人造材料之一,在中国经历了从舶来品到自主生产的过程.玻璃外观虽相似,但根据其添加的助熔剂不同,其主要化学成分也不同1-2.根据一批我国古代玻璃制品的相关数据研究玻璃制品的化学成分及鉴别,对古代玻璃文物的保护具有重要意义3.有学者根据中国古代玻璃的发展,将玻璃成分的演变分为5个阶段:从春秋到战国前期的K2O-C a O-S i O2系统,其中K2O/N a2O1;从战国到东汉时期的B a O-P b O-S i O2系统和K2O-S i O2系统;从东汉到唐代时期的P b O-S i O2系统;从唐代到元代时期的K2O-P b O-S i O2系统;从元代到清代时期的K2O-C a O-S i O2系统4.主成分分析方法已经被大量应用于古代文物的具体成分研究中,特别是古陶瓷、古玻璃器化学成分的分析处理5.此外,一些学者采用外束质子激发X荧光技术(P I X E)、电感耦合等离子体原子发射光谱分析(I C P A E S)方法,对新疆、湖北等地区出土的一批战国时期的玻璃珠(包含镶嵌玻璃珠)、玻璃璧样品进行检测6.史美光等7应用扫描电镜-能量色散X射线分析(S EM-E D X)、电感耦合等离子体发射光谱分析方法、密度测定和偏光显微镜观察等方法,对十多件有代表性的钾硅玻璃样品进行了研究.有学者采用光谱分析、同位素射线荧光分析和比重测试3种方法对西汉到北宋的5 2件玻璃样品进行了分析 8.对玻璃化学成分的研究大多基于物理、化学方法对其成分进行判别并研究其演化历程,仅有少数36教学与竞赛基于成分数据的古代玻璃制品分析与分类2 0 2 3年6月研究采用统计方法.本文从统计视角出发,基于成分数据9-1 0,运用统计方法建立成分预测、玻璃分类及进一步亚分类的模型,进而将这些模型推广,运用至新发现玻璃的成分研究与分类,有助于促进玻璃文物的保护和玻璃制品的生产制造.1 模型假设1)化学成分指采样点处的化学成分;2)未检测到化学成分的原因是仪器精度受限,故暂时将化学成分缺失值记为0,后续进行插补处理;3)检测到的化学成分0值是由于四舍五入得到的,后续进行插补处理.2 数据预处理2.1 附表1数据预处理附表1给出了玻璃文物编号、纹饰、类型、颜色和表面风化的基本信息.由于数据量较大且存在一定的缺失值,故对数据进行缺失值处理,防止对后续建模产生不利影响.其中缺失值均为玻璃文物颜色,按其他信息可分为“纹饰A、表面风化、铅钡”和“纹饰C、表面风化、铅钡”两类,按此分别进行填补.在所给文物中,符合“纹饰C、表面风化、铅钡”条件并已知颜色的1 5件文物颜色分布如表1所示.表1 颜色分布颜色蓝绿浅蓝浅绿深绿紫对应文物编号5 6,5 71 1,2 5,4 3,5 1,5 2,5 44 13 4,3 6,3 8,3 90 8,2 6 因数据为定性数据,故采用热卡来填充缺失值.对1 5件文物及待填充的4 0、5 8号文物的化学成分进行比较,用成分相似的文物颜色进行填充.方法是:一件文物在两个部位进行采样时取均值代表,在一个部位和严重分化点采样时按12的权重计算后作为代表,在一个部位和未分化点采样时将部位的化学成分作为代表.经比较可知,4 0号与3 9号文物成分相似,5 8号与5 1号成分相似,故将4 0和表2 颜色分布颜色黑蓝绿浅蓝对应文物编号4 9,5 02 30 2,2 8,2 9,4 2,4 4,5 35 8号文物颜色分别填充为深绿和浅蓝.符合“纹饰A、表面风化、铅钡”条件并已知颜色的9件文物颜色如表2所示.由附表2可知,9件文物中有6件文物采样点为未风化点,所得到的成分含量可参考性较低,故对此采用众数进行填补,即1 9和4 8号文物颜色填充为浅蓝.2.2 附表2、3数据预处理累加和介于8 5%1 0 5%之间的数据视为有效数据,故1 5号和1 7号文物为无效数据,将其剔除.假设表单2、3中出现的空白为仪器精度受限未检测到该成分,可看作近似零值;另外,假设表单2、3中的0值是由于四舍五入得到,也可看作近似零值.采用乘法替换对近似零值进行插补,具体过程为如下.考虑成分数据集X=xi j nD=x1 1x1 2x1Dx2 1x2 2x2Dxn1xn2xn D ,(1)其中:X中每一行为一个成分数据,每个成分数据含有D个部分,如第i行表示第i个玻璃文物的化学成分数据;xi,j表示第i个玻璃文物的第j个化学成分的含量.假定成分数据中有近似零值,且不同成分数据相同部分对应的探测范围是相同的.记探测范围向量为e=e1,e2,eD T,其中,ej为成分数据集X的第j个部分对应的探测范围.运用乘法简单替换法,对xi j替换后的数据为46第1 2卷 第2期数学建模及其应用V o l.1 2 N o.2 J u n.2 0 2 3xi j=i j,xi j=0,xi j1-k|xi k=0i k /c ,xi j0,(2)其中:i j为一个小于ej的数;c为成分数据的常数和约束,即Dj=1xi j=c.通过实验发现,当成分数据集中近似零值比例不高,i j等于探测范围的6 5%时插补效果最好,即i j=0.6 5ej.记附表2、3中空白处所在列的最小值为临界值,并将其乘以0.6 5,得到插补值,如表3和表4所示(仅列出部分插补值).表3 表单2成分数据插补缺失值表(部分)%文物采样点S i O2N a2OK2OS n O2S O20 16 9.3 3 0 00.5 2 7 59.9 9 0 00.1 5 1 70.3 9 0 00 3部位18 7.0 5 0 00.5 2 8 75.1 9 0 00.1 5 2 00.0 7 2 70 3部位26 1.7 1 0 00.5 2 3 91 2.3 7 0 00.1 5 0 60.0 7 2 05 72 5.4 2 0 00.5 2 6 40.0 7 2 40.1 5 1 30.0 7 2 45 83 0.3 9 0 00.5 2 3 90.3 4 0 00.1 5 0 60.0 7 2 0表4 表单3成分数据插补缺失值表(部分)%文物编号表面风化S i O2N a2OK2OS n O2S O2A 1无风化7 8.4 5 0 00.5 2 7 70.0 7 2 60.1 5 1 70.5 1 0 0A 2风化3 7.7 5 0 00.5 2 9 80.0 7 2 80.1 5 2 30.0 7 2 8A 3无风化3 1.9 5 0 00.5 2 3 91.3 6 0 00.1 5 0 60.0 7 2 0A 4无风化3 5.4 7 0 00.5 2 4 00.7 9 0 00.1 5 0 70.0 7 2 1A 5风化6 4.2 9 0 01.2 0 0 00.3 7 0 00.4 9 0 00.0 7 1 6A 6风化9 3.1 7 0 00.5 2 7 81.3 5 0 00.1 5 1 70.0 7 2 6A 7风化9 0.8 3 0 00.5 2 8 10.9 8 0 00.1 5 1 80.1 1 0 0A 8无风化5 1.1 2 0 00.5 2 4 80.2 3 0 00.1 5 0 92.2 6 0 0因成分数据相加需为1 0 0,故利用相对信息对数据进一步处理.其中相对信息是指成分数据仅有的信息反映在成分间的比率中,与每个成分的绝对数据是无关的.因为成分数据的每个成分乘以相同的正常数,成分间的比率是不变的,因此成分数据可以看作等价类,且类中数据含有相同的信息,可以通过适合的尺度因子表示为相同比例向量.闭合运算是对初始向量乘以合适的尺度因子,使得闭合后的成分和为常数k(这里为1 0 0),定义为C(x)=C(x1,x2,xD)T=(kx1/Di=1xi,kx2/Di=1xi,kxD/Di=1xi)T.(3)对于任意两个向量x,yRD+,如果C(x)=C(y),则x和 y是成分等价的.故将每个玻璃文物化学成表5 虚拟变量处理表变量处理结果变量处理结果纹饰A1浅蓝色1纹饰B2深蓝色2纹饰C3蓝绿色3铅钡1浅绿色4高钾2深绿色5无风化0绿色6风化1紫色7黑色8分乘以计算出的相应因子得到最终数据.3 模型建立与求解3.1 问题13.1.1 玻璃表面风化与纹饰、玻璃类型、颜色的相关性分析首先进行S p e a r m a n相关系数分析.在分析之前对定性数据进行虚拟变量处理,见表5.利用S P S S软件分别对纹饰、玻璃类型、颜色与表面风化进行S p e a r m a n相关性检验,结果见表6.表6结果表明,纹饰、颜色与表面是否风化无相关关系,而玻璃类型与其存在相关关系.56教学与竞赛基于成分数据的古代玻璃制品分析与分类2 0 2 3年6月表6 其他变量与表面风化相关性表纹饰表面风化玻璃类型表面风化颜色表面风化纹饰1.0 0 00.0 8 0玻璃类型1.0 0 0-0.3 0 1颜色1.0 0 00.0 8 8表面风化0.0 8 01.0 0 0表面风化-0.3 0 11.0 0 0表面风化0.0 8 81.0 0 0 其次,为保证结果的可信度,利用卡方检验对数据进行分析,在分析之前对数据进行频数统计,见图1-图3.图1 纹饰与表面风化的频数统计图2 玻璃类型与表面风化的频数统计图3 颜色与表面风化的频数统计由图1可知,在无风化与风化两种情况下,纹饰类型A、B、C变化差异较小,故可初步推断表面风化与纹饰类型无关;由图2可知,风化与无风化相比,铅钡玻璃类型占比增加,高钾玻璃类型占比减少,且变化幅度较大,故可初步推断表面风化与玻璃类型有关;由图3可知,风化与无风化相比,各表7 卡方检验结果表变量卡方检验p值纹饰与表面风化0.0 8 5玻璃类型与表面风化0.0 2 4颜色与表面风化0.3 2 5种颜色变化较小,故可初步推断表面风化与颜色无关.为进一步验证推断是否合理,利用S P S S软件进行卡方检验,结果见表7.因只有玻璃类型与表面风化的卡方检验p值小于0.0 5,故玻璃类型与表面风化有相关关系,纹饰、颜色与表面风化的卡方检验p值均大于0.0 5,故纹饰、颜色与表面风化无相关关系.3.1.2 不同类型玻璃表面有无风化时化学成分含量的统计规律分析运用附表2处理后的成分数据进行分析,分别计算出铅钡、高钾玻璃的无风化与风化检测点数据均值.具体计算方法如下.因成分数据为几何结构,故可对其进行扰动运算(类似实数空间上的加法运算).对于任意成分数据x=(x1,x2,xD)T,y=(y1,y2,yD)TSD,x与y的扰动运算定义为:xy=C(x1y1,x2y2,xDyD)TSD.(4)由扰动运算可进一步计算均值.运用R语言对成分数据进行均值计算,结果如表8所示(仅列出部分化学成分).表8 高钾玻璃、铅钡玻璃风化与无风化均值表(部分)%类型S i O2N a2OK2OS n O2S O2高钾未风化7 4.9 5 2 90.8 8 1 87.2 3 2 60.2 1 1 00.1 2 3 4高钾风化9 3.6 6 1 50.5 2 6 80.3 5 7 90.1 5 1 50.0 7 2 4铅钡未风化5 9.5 0 4 61.3 0 9 00.1 8 5 10.1 8 7 10.0 9 5 3铅钡风化2 7.0 3 1 70.7 4 3 20.1 5 0 20.2 1 0 50.1 7 4 5 为更直观地观察化学成分的变化,将结果可视化,如图4和图5所示.其中,外圈为风化,内圈为无风化.对于高钾玻璃,由图4可知风化后S i O2占比增加较大,K2O显著减少,A l2O3和C a O占比均有所下降;对于铅钡玻璃,由图5可知风化与无风化相比S i O2占比降低,P b O、B a O、P2O5和C a O占比升高.66第1 2卷 第2期数学建模及其应用V o l.1 2 N o.2 J u n.2 0 2 3图4 高钾玻璃风化前后元素对比图(彩图见封三)图5 铅钡玻璃风化前后元素对比图(彩图见封三)3.1.3 风化前化学成分预测首先对附表1中的数据进行取虚拟变量处理,见表5;再将附表1中的变量作为自变量,附表2中的化学成分作为因变量,作D i r i c h l e t回归.得到回归模型后,将自变量中的风化(1)替换为未风化(0),得到预测结果,如表9所示(因预测结果数据量较大,故此处只给出了部分预测结果).表9 预测风化前的化学成分含量表(部分)%文物编号S i O2N a2OK2OS n O2S O20 23 2.8 1 22.5 2 51.3 5 91.0 1 60.8 8 30 77 2.5 0 31.8 1 62.8 3 00.9 4 40.8 1 30 84 6.8 2 22.5 6 71.6 7 71.2 3 40.9 6 45 74 3.5 3 82.4 0 61.2 5 31.0 9 90.9 7 25 84 1.4 6 42.3 0 11.0 7 01.0 2 50.9 6 33.2 问题23.2.1 玻璃的分类规律分析决策树是根据已有样本的信息与现有分类,基于信息熵下降最快的原则,通过不断训练样本选择分类特征,最终得出包含节点与有向边的树状分类模型,可根据此模型对样本类别进行判断1 1.该方法易于理解且便于结果可视化,同时对于小样本具有很好的分类效果.考虑到本文样本量较小且需要选择不同类型玻璃的分类特征,故选用决策树模型.在进行建模之前首先要对成分数据进行中心对数比(c e n t e r e d l o g-r a t i o,C L R)变换.对于任意成分数据x=(x1,x2,xD)TSD,c l r变换将xSD变换为RD上的系数,c l r系数为c l r(x)=l o gx1gm(x),l o gx2gm(x),l o gxDgm(x)T.(5)记c l r变换后的数据为c l r(x)=(1,2,D)T,则c l r逆变换为x=c l r-1()=C(e x p1 ,e x p2 ,e x pD )T,(6)用转化后的成分数据建立分类模型.为使分类结果直观明了,本文使用决策树模型对其进行特征的选择,进一步探究主要分类特征,从而推断两种玻璃的分类规律.在此模型中,将风化、S i O2、N a2O、K2O、C a O、M g O、A l2O3、F e2O3、C u O、P b O、B a O、P2O5、S r O、S n O2和S O2作为分类特征,将玻璃的类型作为类别构建决策树模型,且将6 7个样本划分为训练集与测试集,其中4 7个样本为训练集,图6 决策树结果剩余则为测试集.最后采用R语言进行软件实现,结果如图6所 示.再 用 测 试 集 进 行 预 测,正 确 率 为1 0 0%.由图6可知,最终以P b O为特征进行两类玻璃分类,若P b O值大于1.3则为铅钡玻璃,反之为高钾玻璃.76教学与竞赛基于成分数据的古代玻璃制品分析与分类2 0 2 3年6月偏最小二乘判别分析是一种有监督的判别分析方法,可视为回归模型的扩展,主要用于多变量分析技术中的判别分析1 2.与偏最小二乘回归相比,偏最小二乘判别分析仅需1个数据集,但需要提前知道数据的标签,其余与偏最小二乘回归原理相同1 3.在此模型中,因变量为给定的类标签,自变量与传统回归模型相同,用此回归结果可推断出对因变量分类有较大影响的特征,从而可根据该模型对新给定的样本进行所属类别判断.通过计算不同变量的投影重要性(V I P)可反映每个变量的加载权重以及其对因变量的解释程度,从而通过该方法选择出对分类有较大影响的变量1 4.下面采用偏最小二乘判别分析对高钾玻璃和铅钡玻璃进行分类.由于c l r变换后数据求和为0,以玻璃类型为分类变量,风化、化学成分为特征变量,建立偏最小二乘判别分析,结果如图7所示,分析模型验证结果如图8所示.由图7可知,高钾玻璃和铅钡玻璃有明显的分离趋势.由图8可知,Q 2左侧交于Y轴负半轴,说明模型构建成功.图7 偏最小二乘判别分析结果图8 偏最小二乘判别分析模型验证表1 0 偏最小二乘判别分析的不同特征的V I P值变量V I P变量V I P风化0.2 0 9 8C u O0.5 6 9 4S i O20.8 9 0 7P b O2.4 6 1 4N a2O0.1 2 0 3B a O1.3 4 7 2K2O1.8 0 0 0P2O50.1 2 2 6C a O0.5 0 9 5S r O1.1 1 1 3M g O0.3 8 7 1S n O20.3 5 4 8A l2O30.6 1 5 7S O20.1 2 4 7F e2O30.7 2 8 4 计算不同变量的投影重要性V I P值,结果如表1 0所示.筛选V I P值大于1的特征,对于分类有影响的特征为P b O、K2O、B a O和S r O.上 面 图 表 结 果 是 在s i m c a软件操作完成.3.2.2 亚类划分方法结果及敏感性分析K-m e a n s是一种根据样本距离将样本分类的无监督学习方法,首先把样品随机分成k个初始类;再不断进行修改迭代最终将样品分到其最近均值的类中去.其中样本间的距离测算采用欧氏距离.欧式距离公式为:d2i j=pk=1xi k-xj k 2.(7)将两类玻璃分别进行更进一步的亚类划分.本次划分并不知亚类划分的结果,该批数据并未有明确的类型,因此该类问题为无监督学习问题,对其进行聚类分析.首先需要选择所分类别的个数,采用R语言中的f v i z_n b c l u s t函数进行判断,结果如图9和图1 0所示.根据其间断点可知两种玻璃的亚分类最终均选取3类.图9 高钾分类个数确定 图1 0 铅钡分类个数确定86第1 2卷 第2期数学建模及其应用V o l.1 2 N o.2 J u n.2 0 2 3 其次,在确定聚类个数的情况下用K-m e a n s进行聚类,结果如图1 1和图1 2所示.图1 1 高钾聚类情况 图1 2 铅钡聚类情况由图1 1、图1 2可知铅钡玻璃可被明显地分为3类,故其亚分类包含3类;高钾玻璃3类划分并未有非常明显的类别,仅可粗略将其分为3类,如表1 1所示.表1 1 亚分类玻璃类型亚分类文物编号高钾玻璃10 6部位1,1 822 1,0 7,0 9,1 0,1 2,2 2,2 730 1,0 3部位1,0 3部位2,0 4,0 5,0 6部位2,1 3,1 4,1 6铅钡玻璃12 0,3 7,5 0未风化点,0 8,0 8严重风化,1 1,1 9,2 6,2 6严重风化,3 9,4 0,4 3部位2,5 0,5 1部位1,5 2,5 4,5 4严重风化,5 6,5 822 8未风化点,2 9未风化点,3 0部位1,3 0部位2,3 1,3 2,3 5,4 9未风化点,0 2,4 1,4 8,4 9,5 1部位232 3未风化点,2 4,2 5未风化点,3 3,4 2未风化点1,4 2未风化点2,4 4未风化点,4 5,4 6,4 7,5 3未风化点,5 5,3 4,3 6,3 8,4 3部位1,5 7在确定分类个数时,不同的选择会出现不同的结果,故应对此模型的敏感性进行分析.因高钾玻璃并未有明显的分类效果,故尝试多种分类个数进行聚类,从而选择出较优的聚类情况.下面对每类玻璃的亚类划分选择合适的化学成分.对于高钾玻璃,亚类划分为3类.由于其中一类只有两种文物玻璃,因此删除这一类文物,以其余类为分类变量,化学成分为特征,构建偏最小二乘判别分析,如图1 3和图1 4所示.由图1 3可知,高钾玻璃的亚类有明显的分离趋势.由图1 4可知,Q 2左侧交于Y轴负半轴,说明模型构建成功.图1 3 高钾玻璃亚类偏最小二乘判别分析结果图1 4 高钾玻璃亚类偏最小二乘判别分析模型验证 计算不同变量的V I P值,结果见表1 2.筛选V I P值大于1的特征,表明对于高钾玻璃亚类分类有影响的特征分别为K2O、S i O2、C a O、S n O2和B a O.表1 2 高钾玻璃亚类偏最小二乘判别分析的不同特征的V I P值表变量V I P变量V I P变量V I P变量V I PS i O21.5 2 6 7C u O0.7 5 0 9M g O0.0 2 1 6S r O0.4 1 4 7N a2O0.3 8 1 5P b O0.2 0 6 8A l2O30.4 6 5 4S n O21.2 7 6 3K2O2.3 0 1 6B a O1.2 1 5 8F e2O30.5 2 3 5S O20.3 7 2 1C a O1.3 0 7 8P2O50.0 1 6 096教学与竞赛基于成分数据的古代玻璃制品分析与分类2 0 2 3年6月 对于铅钡玻璃,亚类划分为3类.以这3类为分类变量,化学成分为特征,构建偏最小二乘判别分析,结果如图1 5和1 6所示.由图1 5可知,铅钡玻璃的3类有明显的分离趋势.由图1 6可知,Q 2左侧交于Y轴负半轴,说明模型构建成功.图1 5 铅钡玻璃亚类偏最小二乘判别分析结果图1 6 铅钡玻璃亚类偏最小二乘判别分析模型验证表1 3 铅钡玻璃亚类偏最小二乘判别分析的不同特征的V I P值表变量V I P变量V I PS i O20.9 5 5 0C u O1.3 7 1 1N a2O1.3 0 4 6P b O0.2 8 1 3K2O0.5 8 2 5B a O0.8 6 0 5C a O0.6 6 9 6P2O52.0 5 1 2M g O0.7 7 1 6S r O0.4 4 3 5A l2O30.8 3 3 9S n O20.4 6 6 3F e2O31.2 1 1 9S O20.7 1 9 3 计算不同变量的V I P值,结果见表1 3.筛选V I P值大于1的特征,对于铅钡玻璃亚类分类有影响的特征分别为P2O5、C u O、N a2O和F e2O3.3.3 问题3方法一:表单3中未知玻璃类型的类别划分.基于决策树的结果,以P b O为特征进行分类,若P b O值大于1.5则为铅钡玻璃,反之为高钾玻璃.结果表明,A 1、A 6、A 7为高钾玻璃,A 2、A 3、A 4、A 5、A 8为铅钡玻璃.方法二:选取问题2中对于高钾玻璃和铅钡玻璃筛选的特征P b O、K2O、B a O和S r O,建立这4个特征与玻璃类型的偏最小二乘回归分析,其中因变量玻璃类型中高钾玻璃取值为1,铅钡玻璃取值为0.表1 4 不同主成分个数下偏最小二乘回归结果表评价准则1个2个3个4个交叉验证的预测均方根误差0.1 4 7 1 0.1 3 8 9 0.1 3 7 1 0.1 3 7 2调整的交叉验证的预测均方根误差0.1 4 7 0 0.1 3 8 1 0.1 3 6 5 0.1 3 6 5累计贡献率/%8 0.1 98 7.1 59 3.5 51 0 0基于交叉验证方法计算预测均方根误差,使用所有主成分进行回归.得到的结果如表1 4所示.从回归结果可以看出,主成分个数为3时,模型在经交叉验证后得到的预测均方根误差最小,同时3个主成分对各变量的累计贡献率已经达到9 3%,因此将偏最小二乘回归的主成分个数设定为3.主成分个数确定后,计算得到P b O、K2O、B a O和S r O的偏最小二乘回归系数分别为0.0 1 3 4、-0.1 5 8 2、-0.0 2 8 8和-0.0 2 0 6.将表单3中P b O、K2O、B a O和S r O的数据代入偏最小二乘回归模型,得到不同文物的预测值,若预测值接近1,为高钾玻璃;预测值接近0,为铅钡玻璃.预测结果见表1 5.表1 5 表单3未知玻璃文物的类型预测文物编号预测值预测玻璃类型A 11.0 4 0 8高钾A 20.0 1 7 7铅钡A 30.0 9 8 0铅钡A 40.1 8 9 8铅钡A 50.2 5 5 8铅钡A 60.9 6 1 5高钾A 70.9 6 3 1高钾A 80.1 1 2 2铅钡通过上述分析,可以看出两种方法预测结果一致.为了更进一步分析高钾玻璃与铅钡玻璃的亚类划分,基于上述亚类划分结果:1)对于高钾玻璃,以3个亚类为因变量,化学成分K2O、S i O2、C a O、S n O2和B a O为自变量,建立偏最小二乘回归,通过交叉验证确定主成分个数为2,对表单3中A 1、A 6和A 7文物进行预测,结果表明A 1、A 6和A 7都是高钾玻璃的同一亚类;2)对于铅钡玻璃,以3个亚类为因变量,化学成分P2O5、07第1 2卷 第2期数学建模及其应用V o l.1 2 N o.2 J u n.2 0 2 3C u O、N a2O和F e2O3为自变量,建立偏最小二乘回归,通过交叉验证确定主成分个数为2,对表单3中A 2、A 3、A 4、A 5和A 8文物进行预测,结果表明A 2和A 4是铅钡玻璃的同一亚类,A 3和A 8是铅钡玻璃的同一亚类,A 5是铅钡玻璃的另一亚类.3.4 问题43.4.1 相关性分析灰色系统着重内涵不明确的对象,将不确定的“灰”转化为可解读的“白”信息1 5.其中灰色关联通过序列曲线集合形状可直观地对关联度进行分析.将玻璃按照高钾玻璃与铅钡玻璃分类讨论.因为题中要求分析不同类别文物化学成分之间的关联关系,且变量较多,热力图更可直观地体现两两化学成分之间的相关关系,故分别制作高钾玻璃与铅钡玻璃的化学元素热力图,如图1 7所示.图1 7 高钾与铅钡玻璃热力图由图1 7(a)可以看出,按照热力图上方从左到右元素顺序,每个元素与左右相邻元素相关性较强,S i O2与B a O关联性较强,其余两两成分之间关联性相对较弱.由图1 7(b)可以看出,S n O2、S i O2、N a2O、F e2O3、K2O、M g O、A l2O3这7个成分两两之间均存在相关性(F e2O3和N a2O之间除外).C u O、B a O、S O2两两之间均存在较为强烈的相关性.表1 6 两类玻璃化学成分关联度对比表评价项(铅钡)关联度排名评价项(高钾)关联度排名 N a2O0.7 9 0 41 1N a2O0.6 2 9 11 1K2O0.9 6 3 95K2O0.8 9 9 18C a O0.8 2 2 71 0C a O0.9 1 3 27M g O0.9 2 0 59M g O0.8 5 7 19A l2O30.9 6 0 36A l2O30.9 7 6 52F e2O30.9 2 7 88F e2O30.5 6 8 51 2C u O0.6 7 4 21 2C u O0.9 2 0 46P b O0.9 7 8 51P b O0.9 5 7 74B a O0.9 6 9 44B a O0.5 6 5 51 3P2O50.5 4 0 01 3P2O50.8 3 4 11 0S r O0.9 5 8 77S r O0.9 7 9 81S n O20.9 7 2 92S n O20.9 4 9 25S O20.9 6 9 83S O20.9 6 7 33考 虑 到 玻 璃 的 主 要 化 学 成 分 为S i O2,故将S i O2作为母序列,其他化学成分作为子序列,运用M a t l a b进行灰色关联分析,最终得出关联度值.使用关联度值对4 9个评价对象进行排序,结果如表1 6所示.关联度值介于01之间,该值越大代表其与母序列S i O2之间的相关性越强.从表1 6可以看出:本次1 3个评价项中,铅钡玻璃中,P b O的综合评价最高(关联度为0.9 7 9 8),P2O5的综合评价最低(关联度为0.5 4 0 0);高钾玻璃中,S r O的 综 合 评 价 最 高(关 联 度 为0.9 7 8 5),B a O的综合评价最低(关联度为0.5 6 5 5).3.4.2 差异性分析比较不同类别之间的化学成分关联关系的差异性.由于热力图中相关系数为对称矩阵,因此仅选取高钾玻璃相关系数上三角矩阵与铅钡玻璃相关系数上三角矩阵.由于每个化学成分之间的相关关系17教学与竞赛基于成分数据的古代玻璃制品分析与分类2 0 2 3年6月是配对的,因此选取配对样本的非参数w i l c o x o n检验,所得p值大于(0.0 5),因此不拒绝原假设,即两种类型玻璃的化学成分之间的关联关系没有显著差异.通过R作出关于铅钡玻璃与高钾玻璃的关联度折线图和箱线图,如图1 8所示.由图1 8可知,高钾玻璃与铅钡玻璃相比分散程度较小,且通过表1 6的排名可知,两者的化学成分之间的关联关系没有显著差异,进一步验证了上述结论.图1 8 铅钡玻璃与高钾玻璃的关联度折线图和箱线图4 模型检验及可靠性分析4.1 针对问题1的检验在解答问题1之前,对表单中数据首先进行预处理.对于玻璃的化学成分,近似零值插补后转换为成分数据,化学成分比例和为1 0 0%,数据预处理合理.对于问题1,考虑到成分数据的特殊结构,在成分数据单形空间上计算均值,选择适用于成分数据的D i r i c h l e t回归模型,分析方法相比传统分析方法更加合理.4.2 针对问题2的检验对于问题2,选择两种方法分析高钾玻璃与铅钡玻璃的分类规律,两种方法结果一致,进一步验证了分类模型的合理性.对每类玻璃亚类划分时,基于K-m e a n s聚类分析确定了最优聚类个数,分析结果真实可靠.4.3 针对问题3的检验偏最小二乘回归方法基于交叉验证确定了最优主成分个数.基于两种方法对未知玻璃文物进行类别预测,结果一致,因此预测结果合理.4.4 针对问题4的检验采用P e a r s o n相关系数分析化学成分之间的关联关系,基于w i l c o x o n检验比较两种类型玻璃的化学成分之间关联关系的差异性.结果真实可靠.5 模型评价与展望5.1 模型的优点本文基于成分数据对玻璃化学成分进行分析,并对不同玻璃类型进行分类,主要优点如下:1)采用D i r i c h l e t回归模型对风化前的化学成分进行预测;2)对于玻璃类型的分类,考虑到成分数据求和为1 0 0%的约束,在模型构建前,首先对化学成分进行c l r变换,使得成分数据变换为欧式空间上的普通数据;3)采用不同分类模型对玻璃类型进行分类,不同模型结果一致.5.2 模型的缺点对于问题4,不同玻璃类型的化学成分之间的关联关系,采用P e a r s o n相关系数来说明.但P e a r s o n相关系数只能度量变量之间的线性相关关系,在使用之前未进行线性相关关系检验.5.3 模型的展望对于问题2,分析不同玻璃之间的分类规律时,由于两种类型玻璃的样本量不是很接近,因此后续可以考虑不平衡样本分类模型,通过对训练集样本重采样或方法修正来进行分类.对于问题4,可27第1 2卷 第2期数学建模及其应用V o l.1 2 N o.2 J u n.2 0 2 3考虑利用其他相关系数来度量化学成分之间的关联关系,例如最大距离相关系数、互信息等.参考文献1 赵志强.新疆巴里坤石人子沟遗址群出土玻璃珠的成分体系与制作工艺研究D.西安:西北大学,2 0 1 6.2 安家瑶.玻璃器史话M.北京:社会科学文献出版社,2 0 1 1:7-1 1.3 全国大学生数学建模组委会.2 0 2 2“高教社杯”全国大学生数学建模竞赛赛题E B/O L.2 0 2 2-0 9-1 5.h t t p:/w w w.m c m.e d u.c n/h t m l_c n/n o d e/5 2 6 7 f e 3 e 6 a 5 1 2 b e c 7 9 3 d 7 1 f 2 b 2 0 6 1 4 9 7.h t m l.4 干福熹.中国古代玻璃的起源和发展J.自然杂志,2 0 0 6,2 8(4):1 8 7-1 9 3.5B a x t e r M J.E x p l o r a t o r y m u l t i v a r i a t e a n a l y s i s i n a r c h a e o l o g yM.E d i n b u r g h:E d i n b u r g h U n i v e r s i t y P r e s s,1 9 9 4:4 8-5 2.6 李青会,黄教珍,李飞,等.中国出土的一批战国古玻璃样品化学成分的检测J.文物保护与考古科学,2 0 0 6,1 8(2):8-1 3.7 史美光,何欧里,周福征.一批中国汉墓出土钾玻璃的研究J.硅酸盐学报,1 9 8 6,1 4(3):3 0 7-3 1 3.8 史美光,曲长芝,张日清,等.中国早期玻璃器检验报告J.考古学报,1 9 8 4,4:4 4 9-4 5 7.9P a w l o w s k y-G l a h n V,B u c c i a n t i A.C o m p o s i t i o n a l d a t a a n a l y s i s:t h e o r y a n d a p p l i c a t i o n sM.C h i c h e s t e r:W i l e y,2 0 1 1.1 0P a w l o w s k y-G l a h n V,E g o z c u e J J,T o l o s a n a-D e l g a d o R.M o d e l i n g a n d a n a l y s i s o f c o m p o s i t i o n a l d a t aM.C h i c h e s