2022年国赛高教杯数学建模C题古代玻璃制品的成分分析与鉴别解题全过程文档及程序

2022年国赛高教杯数学建模

C题 古代玻璃制品的成分分析与鉴别

原题再现

丝绸之路是古代中西方文化交流的通道，其中玻璃是早期贸易往来的宝贵物证。早期的玻璃在西亚和埃及地区常被制作成珠形饰品传入我国，我国古代玻璃吸收其技术后在本土就地取材制作，因此与外来的玻璃制品外观相似，但化学成分却不相同。
玻璃的主要原料是石英砂，主要化学成分是二氧化硅（SiO2）。由于纯石英砂的熔点较高，为了降低熔化温度，在炼制时需要添加助熔剂。古代常用的助熔剂有草木灰、天然泡碱、硝石和铅矿石等，并添加石灰石作为稳定剂，石灰石煅烧以后转化为氧化钙（CaO）。添加的助熔剂不同，其主要化学成分也不同。例如，铅钡玻璃在烧制过程中加入铅矿石作为助熔剂，其氧化铅（PbO）、氧化钡（BaO）的含量较高，通常被认为是我国自己发明的玻璃品种，楚文化的玻璃就是以铅钡玻璃为主。钾玻璃是以含钾量高的物质如草木灰作为助熔剂烧制而成的，主要流行于我国岭南以及东南亚和印度等区域。
古代玻璃极易受埋藏环境的影响而风化。在风化过程中，内部元素与环境元素进行大量交换，导致其成分比例发生变化，从而影响对其类别的正确判断。如图 1 的文物标记为表面无风化，表面能明显看出文物的颜色、纹饰，但不排除局部有较浅的风化；图 2 的文物标记为表面风化，表面大面积灰黄色区域为风化层，是明显风化区域，紫色部分是一般风化表面。在部分风化的文物中，其表面也有未风化的区域。

现有一批我国古代玻璃制品的相关数据，考古工作者依据这些文物样品的化学成分和其他检测手段已将其分为高钾玻璃和铅钡玻璃两种类型。附件表单 1 给出了这些文物的分类信息，附件表单 2 给出了相应的主要成分所占比例（空白处表示未检测到该成分）。这些数据的特点是成分性，即各成分比例的累加和应为 100%，但因检测手段等原因可能导致其成分比例的累加和非 100%的情况。本题中将成分比例累加和介于 85%~105%之间的数据视为有效数据。
请你们团队依据附件中的相关数据进行分析建模，解决以下问题：
问题 1 对这些玻璃文物的表面风化与其玻璃类型、纹饰和颜色的关系进行分析；结合玻璃的类型，分析文物样品表面有无风化化学成分含量的统计规律，并根据风化点检测数据，预测其风化前的化学成分含量。
问题 2 依据附件数据分析高钾玻璃、铅钡玻璃的分类规律；对于每个类别选择合适的化学成分对其进行亚类划分，给出具体的划分方法及划分结果，并对分类结果的合理性和敏感性进行分析。
问题 3 对附件表单 3 中未知类别玻璃文物的化学成分进行分析，鉴别其所属类型，并对分类结果的敏感性进行分析。
问题 4 针对不同类别的玻璃文物样品，分析其化学成分之间的关联关系，并比较不同类别之间的化学成分关联关系的差异性。

整体求解过程概述(摘要)

玻璃是人类最早批发明的人造材料之一，古代玻璃除少数天然玻璃外，其余均为人造玻璃。虽外观相似，但化学成分不同。研究玻璃制品的化学成分及鉴别对古代玻璃文物的保护具有重要意义。本文从成分数据分析视角出发，对玻璃的化学成分进行分析，并对不同类型的玻璃进行分类及预测。
首先对附件表单中数据进行预处理。表单 1 数据按照众数以及热卡填充进行缺失值插补。由于玻璃制品的化学成分属于成分数据，但因检测手段等原因可能导致其成分比例的累加和非 100%，另外很多成分检测不到或为 0 值。假设这种情况是由于仪器精度首先或四舍五入得到的。因此，选用乘法替换法对表单 2 和 3 中空自处及 0 值进行插补，最后得到的玻璃化学成分数据为成分数据
对于问题 1，首先基于 spearman 相关系数以及卡方检验分析玻璃文物表面风化与类型、纹饰、颜色的关系，结果均表明玻璃文物表面风化与玻璃类型有关，与纹饰、颜色无关。其次，根据成分数据的 Aitchison 几何结构，在单形空间分别计算每种玻璃风化与无风化的化学成分均值，结果表明高钾玻璃风化后二氧化硅占比增加较大，氧化钾显著减少，氧化钠、氧化铝、氧化铜、氧化钡占比均有所上升铅钡玻璃风化后二氧化硅占比降低，氧化铅、氧化钾升高，五氧化二磷降低。最后，以玻璃的化学成分为响应变量，纹饰、类型、颜色与表面风化为自变量，构建 Dirichle 回归模型，进而根据该模型预测风化点风化前的化学成分含量，具体结果见支撑材料。
对于问题 2，首先分析高钾玻璃与铅钡玻璃的分类规律，基于问题 1 结果，选择表面风化、化学成分为分类特征，以玻璃类型为类别。由于化学成分为成分数据，因此对化学成分做对称对数比率(𝑐𝑙𝑟)变换。采用决策树建立分类模型，以氧化铅为特征进行两类玻璃分类，若氧化铅值大于 1.5 则为铅钡玻璃，反之为高钾玻璃。同时建立偏最小二乘判别分析(PLS-DA)，通过不同特征的投影重要性，得到对分类有影响的特征分别为氧化铅、氧化钾、氧化钡和氧化锶。两种方法得到结果一致。其次，采用Kmeans 聚类分析分别对两种类型的玻璃进行聚类，确定最优聚类个数都为 3，然后基于 PLS-DA 对每种玻璃选择合适的化学成分，结果为高钾玻璃以氧化钾、二氧化硅、氧化钙、氧化锡和氧化钡为分类特征分为三类，铅钡玻璃以五氧化二磷、氧化铜、氧化钠和氧化铁为分类特征分为三类。
对于问题 3，依据问题 2 构建模型对未知类别的玻璃进行划分。方法一是基于决策树模型，选择氧化铅成分值对未知玻璃进行划分。方法二是对表单 3 数据做 cr 变换，选择氧化铅、氧化钾、氧化钡和氧化锶成分为自变量，玻璃类型为因变量，建立偏最小二乘回归(PLS)，通过交叉验证法选择主成分个数为 3，通过预测值确定玻璃类型。两种方法得到的结果一致，A1、A6、A7 为高钾玻璃，A2、A3、A4、A5、A8 为铅钡玻璃。进一步结合问题 2 中两类玻璃的亚类划分特征，基于 PLS 得到 A1、A6、A7 都是高钾玻璃的同一亚类，A2 和 A4 是铅钡玻璃的同一亚类，A3 和 A8 是铅钡玻璃的同一亚类，A5 是铅钡玻璃的另一亚类。
对于问题 4，首先分别对每种类型玻璃的化学成分计算相关系数，观察和分析相关性热力图。然后对两类玻璃化学成分之间的相关系数进行配对 Wilcoxon 检验，结果为两类玻璃化学成分之间的相关关系无显著差异。

模型假设：

1.假设化学成分指采样点处的化学成分。
2.假设未检测到的化学成分的原因为仪器精度受限,故暂时将化学成分缺失值记为 0，后续进行插补处理。
3.假设检测到的化学成分 0 值是由于四舍五入得到的，后续进行插补处理。

问题分析：

在回答问题之前首先对附件表单中数据进行预处理。表单 1 数据按照众数以及热卡填充进行缺失值插补。表单 2 和 3 中空白处及 0 值，选用乘法替换法进行插补，然后将化学成分数据转化为成分数据 1。
问题 1
首先基于 spearman 相关系数以及卡方检验分析玻璃文物表面风化与类型、纹饰、颜色的关系。其次，根据成分数据的 Aitchison 几何结构,在单形空间分别计算每种玻璃风化与无风化的化学成分均值，再通过环形图进行比较分析。最后，以文物的化学成分为响应变量，纹饰、类型、颜色与表面风化为自变量，构建 Dirichlet 回归模型，进而根据该模型预测风化点风化前的化学成分含量。
问题 2
为研究高钾玻璃与铅钡玻璃的分类规律，基于问题 1 的结果，选择玻璃基本信息、化学成分为分类特征，以玻璃类型为类别。由于化学成分为成分数据，因此对化学成分做对称对数比率(𝑐𝑙𝑟)变换。分别采用决策树和偏最小二乘判别分析进行分类，选择对于分类重要的特征。对每种类型的玻璃的亚类划分，采用 kmeans 聚类分析分别对两种类型的玻璃进行聚类分析，然后基于偏最小二乘判别分析对每种玻璃选择合适的化学成分。
问题 3
依据问题 2 构建的模型对未知类别的玻璃文物进行类别划分。方法一是基于决策树模型，选择筛选的特征对未知玻璃进行划分。方法二是对表单 3 数据做𝑐𝑙𝑟变换，选择偏最小二乘判别分析筛选的特征为自变量，玻璃类型为因变量，建立偏最小二乘回归，通过交叉验证法确定主成分个数,通过预测值确定玻璃类型。
问题 4
对于每种玻璃类型,我们首先采用 Pearson 相关系数分析化学成分之间的相关关系，并将此相关性以热力图的形式进行数据可视化展示，进而通过配对 Wilcoxon 检验确定两类玻璃化学成分之间的相关关系有无显著差异。

模型的建立与求解整体论文缩略图

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

程序代码：(代码和文档not free)

The actual procedure is shown in the screenshot

import openpyxl
from collections import Counter

defgetOriginData(min_row, max_row, min_col, max_col):
    workbook = openpyxl.load_workbook('data.xlsx').worksheets[0]

    data =[]for row in workbook.iter_rows(min_row=min_row, max_row=max_row, min_col=min_col, max_col=max_col):
        line =[]for cell in row:
            line.append(cell.value)
        data.append(line[0])return data

line1, line2, line3, line4, line5 =[],[],[],[],[]for i inrange(1,6):if i ==1:
        line1 = getOriginData(2,59, i, i)elif i ==2:
        line2 = getOriginData(2,59, i, i)elif i ==3:
        line3 = getOriginData(2,59, i, i)elif i ==4:
        line4 = getOriginData(2,59, i, i)elif i ==5:
        line5 = getOriginData(2,59, i, i)

a_num = line2.count("A")
b_num = line2.count("B")
c_num = line2.count("C")print(a_num, b_num, c_num)

a_fh, a_nfh =0,0
b_fh, b_nfh =0,0
c_fh, c_nfh =0,0for i inrange(len(line2)):if line2[i]=='A':if line5[i]=='风化':
            a_fh +=1else:
            a_nfh +=1elif line2[i]=='B':if line5[i]=='风化':
            b_fh +=1else:
            b_nfh +=1else:if line5[i]=='风化':
            c_fh +=1else:
            c_nfh +=1print(a_fh, a_nfh)print(b_fh, b_nfh)print(c_fh, c_nfh)

type1_num = line3.count("高钾")
type2_num = line3.count("铅钡")print(
    type1_num, type2_num
)

a1_fh, a1_nfh =0,0
b1_fh, b1_nfh =0,0for i inrange(len(line3)):if line3[i]=='高钾':if line5[i]=='风化':
            a1_fh +=1else:
            a1_nfh +=1else:if line5[i]=='风化':
            b1_fh +=1else:
            b1_nfh +=1print(a1_fh, a1_nfh)print(b1_fh, b1_nfh)

fh_color, nfh_color =[],[]for i inrange(len(line5)):if line5[i]=='风化':
        fh_color.append(line4[i])else:
        nfh_color.append(line4[i])print(dict(Counter(fh_color)),dict(Counter(nfh_color)))

import openpyxl

# sheet1
wb1 = openpyxl.load_workbook('data.xlsx').worksheets[0]# sheet2
wb2 = openpyxl.load_workbook('data.xlsx').worksheets[1]defgetOriginData(min_row, max_row, min_col, max_col, who):
    data =[]for row in who.iter_rows(min_row=min_row, max_row=max_row, min_col=min_col, max_col=max_col):
        line =[]for cell in row:
            line.append(cell.value)
        data.append(line)return data

for i inrange(1,4):if i ==1:
        line1 = getOriginData(2,59, i, i, wb1)
        line1 =[item[0]for item in line1]
        line12 = getOriginData(2,70, i, i, wb2)
        line12 =[int(item[0])for item in line12]elif i ==2:
        line22 = getOriginData(2,70, i, i, wb2)
        line22 =[(item[0])for item in line22]elif i ==3:
        line3 = getOriginData(2,59, i, i, wb1)
        line3 =[item[0]for item in line3]
dic ={int(line1[i]):line3[i]for i inrange(len(line1))}

sheet_data = getOriginData(2,70,3,16, wb2)# remove none
new_sheet_data =[]for item in sheet_data:
    temp =[]for w in item:if w ==None:temp.append(0)else:temp.append(w)
    new_sheet_data.append(temp)# prob
pro12_data =[]for i inrange(len(new_sheet_data)):
    temp =[]
    total =sum(new_sheet_data[i])for j in new_sheet_data[i]:
        temp.append(round(j/total,6))
    pro12_data.append(temp)

type1, type2 =[],[]for i inrange(len(line12)):if dic[line12[i]]=='高钾':
        pro12_data[i].insert(0, line22[i])
        type1.append(pro12_data[i])elif dic[line12[i]]=='铅钡':
        pro12_data[i].insert(0, line22[i])
        type2.append(pro12_data[i])# save total prob
wb_type1 = openpyxl.Workbook()
ws_type1 = wb_type1['Sheet']
ws_type1.append(['文物采样点','二氧化硅','氧化钠','氧化钾','氧化钙','氧化镁','氧化铝','氧化铁','氧化铜','氧化铅','氧化钡','五氧化二磷','氧化锶','氧化锡','二氧化硫'])for i inrange(len(type1)):
    d =tuple(type1[i])
    ws_type1.append(d)
wb_type1.save('高钾.xlsx')

wb_type2 = openpyxl.Workbook()
ws_type2 = wb_type2['Sheet']
ws_type2.append(['文物采样点','二氧化硅','氧化钠','氧化钾','氧化钙','氧化镁','氧化铝','氧化铁','氧化铜','氧化铅','氧化钡','五氧化二磷','氧化锶','氧化锡','二氧化硫'])for i inrange(len(type2)):
    d =tuple(type2[i])
    ws_type2.append(d)
wb_type2.save('铅钡.xlsx')# node number
num =[]
point_num =[]for i inrange(len(line12)):if'严重风化点'in line22[i]:
        num.append(line12[i])
        point_num.append(line22[i])

final =[]
s =0for item in[[31,47],[24,46],[45,46]]:
    temp =[]for t in item:for ii inrange(len(line12)):if line12[ii]== t:
                temp.append(new_sheet_data[ii])
    temp1 =[]for m, n inzip(temp[0], temp[1]):
        temp1.append(m + n)
    temp1 =[round(item/2,6)for item in temp1]
    temp1.insert(0, point_num[s])
    final.append(temp1)
    s +=1

wb_final = openpyxl.Workbook()
ws_final = wb_final['Sheet']
ws_final.append(['文物采样点','二氧化硅','氧化钠','氧化钾','氧化钙','氧化镁','氧化铝','氧化铁','氧化铜','氧化铅','氧化钡','五氧化二磷','氧化锶','氧化锡','二氧化硫'])for i inrange(len(final)):
    d =tuple(final[i])
    ws_final.append(d)
wb_final.save('final_predict.xlsx')

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可