古代玻璃制品化学成分的分析与研究

摘要

古代玻璃极易受埋藏环境的影响而风化，并且在风化过程中，内部元素与环境元素进行着大量交换，导致其成分比例会发生变化，从而会影响对其类别的正确判断。玻璃在炼制的过程中需要添加助熔剂和稳定剂，我国古代有常用的助熔剂有草木灰、天然泡碱、硝石和铅矿石等，如果添加的助熔剂不同，其主要化学成分也会不同。研究古代玻璃化学成分，对研究古人当时时期生活面貌、审美观念、 生产能力、工艺技术水平、物产交流情况等有一定的意义[1]。为此，我们建立单因素方差分析，BP神经网络,多元回归和聚类分析等模型对古代玻璃制品成分进行分析和研究。

针对问题一，先将问题一分为三个小问。首先第一小问运用单因素方差分析模型，以表面风化作为因素，以玻璃类型、纹饰和颜色的情况作为检验指标，根据模型原理，得到玻璃类型的显著性p值为0.002,则有显著性影响，而玻璃颜色的显著性p值为0.070和玻璃纹饰的显著性p值为0.078，表现为无差异显著性影响；第二小问我们进行描述性统计分析，计算相关数值的平均数，中位数和标准差，进行折线图分析等，得到铅钡玻璃在风化后主要化学成分含量基本呈上升趋势，而高钾玻璃在风化后的主要化学成分含量基本上呈下降趋势的规律；第三小问利用BP神经网络，根据第二小问已得到的数据变化规律，总结出各个化学成分的变化情况，找到相关函数关系进而进行风化前的数据预测，具体预测值可在附录表19中查看。

针对问题二，同样先把问题分为三个小问。首先分析高钾玻璃和铅钡玻璃的分类规律，根据统计分析分为四类；接着针对铅钡玻璃和高钾玻璃进行亚类划分，先对这些化学成分的标志性数据进行分类汇总，以此作为亚类划分的依据，建立K-Means正在上传…重新上传取消 聚类分析模型，根据模型原理，将原来的四类划分为八个亚类，最终亚类划分结果已在图11呈现；最后进行数据的敏感性检验，也就是灵敏度分析，发现在添加30%的扰动时，铅钡未风化的数据准确率为92.3%，其他的数据准确度均为100%，这一定程度上体现了模型的准确性和合理性。

针对问题三，把该问题分为两个小问。需要先对未知玻璃文物的主要化学成分进行分析，鉴别这些文物属于什么类型，将表单3中的数据进行分类，初步得出这些玻璃的类别，接着对表格1和表格2中的数据进行回归分析，得到不同类型的玻璃各化学成分所占的权重，最终分类结果在正文中表13可知；然后将得到的结果进行灵敏度分析，即可以对这些化学成分中的某一个添加数据扰动，代入问题模型中，结果显示，在添加-1%至30%扰动时，分类结果均正确。

针对问题四，先进行不同玻璃类型的关联性分析，问题三中已得到不同玻璃化学成分所占的比重，我们将这些比重数据进行数据可视化，分析其中的关系，关联性主要体现为各化学成分的类别方面，具体在正文中可看；在解决差异性时，根据问题一所建立的单因素方差分析模型，以玻璃类型作为因素，相对应的14种化学成分作为检验指标，计算得到高钾玻璃在风化前后计算出来的p正在上传…重新上传取消 值大于0.05，说明高钾玻璃在风化后各化学成分之间的关联性变小，而铅钡玻璃风化后计算出来的p 值小于0.05，说明在风化后铅钡玻璃各化学成分之间仍然有较强的相互之间的影响。

关键词：单因素方差分析 BP神经网络 聚类分析 回归分析 灵敏度分析

一、问题的重述

1.1问题背景

丝绸之路是我国古代连接中西方的交流通道，其中我国古代的玻璃成为了我国早期与各国贸易往来的宝贵物证。早期的古代玻璃在西亚和埃及地区常常被制作成珠形饰品传入我国内部，所以我国古代玻璃吸收其西方技术后在本土就地取材制作，但是不同的是我国古代玻璃与外来的玻璃制品外观相似，但是化学成分却大不相同。

玻璃的主要原料是石英砂，它的主要的化学成分是二氧化硅（SiO2正在上传…重新上传取消 ），在炼制的过程中需要添加助熔剂和稳定剂。在我国古代有常用的助熔剂有草木灰、天然泡碱、硝石和铅矿石等。其中以添加石灰石作为稳定剂，石灰石煅烧以后转化为氧化钙（CaO）可以作为稳定剂。如果添加的助熔剂不同，其主要化学成分也会不同。古代玻璃极易受埋藏环境的影响而风化。并且在风化过程中，内部元素与环境元素进行着大量交换，导致其成分比例会发生变化，从而会影响对其类别的正确判断。文物标记可表现为表面风化和表面无风化两种。其中的文物标记为表面无风化的特征是从表面上可以明显看出其文物的颜色和纹饰，但是不排除局部有较浅的风化，而文物标记表现为表面风化的特征为表面大面积灰黄色区域为风化层，其表面有明显的风化区域，紫色部分是一般的风化表面，但是在部分风化的文物中，其表面也有未风化的区域。

1.2问题要求

问题一：首先第一个小问是针对这些玻璃文物的表面风化与其玻璃类型、纹饰和颜色的关系进行分析，其次第二小问是对于结合玻璃的类型，分析文物样品表面有无风化化学成分含量的统计规律，最后第三小问是根据风化点检测数据，预测其风化前的化学成分含量。

问题二：依据附件数据分析高钾玻璃、铅钡玻璃的分类规律；对于每个类别选择合适的化学成分对其进行亚类划分，给出具体的划分方法及划分结果，并对分类结果的合理性和敏感性进行分析。

问题三：对附件表单3中未知类别玻璃文物的化学成分进行分析，鉴别其所属类型，并对分类结果的敏感性进行分。

问题四：针对不同类别的玻璃文物样品，分析其化学成分之间的关联关系，并比较不同类别之间的化学成分关联关系的差异性。

正在上传…重新上传取消

图 1 全文问题分析流程图

二、问题的分析

2.1问题的分析

2.1.1问题一的分析

在解决问题一时，根据题目要求，将题目分为三小问。首先需要对玻璃表面风化情况与玻璃类型、纹饰和颜色的关系进行分析，即进行差异显著性分析，运用单因素方差分析模型，以表面是否风化作为因素，以玻璃类型、纹饰和颜色的情况作为因素所处的水平，分析显著性p值是否小于0.05，进而分析其中的关系；然后结合玻璃的类型分析化学成分含量的变化规律，对此，我们进行描述性统计分析，计算相关数值的平均数，中位数和标准差，进行折线图分析等，总结其中的规律；最后利用BP神经网络，预测风化前的化学成分含量，根据前面已经得到的数据变化规律，总结出各个化学成分的变化情况，找到相关函数关系进而进行风花前的数据预测。

2.1.2问题二的分析

在解决问题二时，根据题目要求，需要我们针对铅钡玻璃和高钾玻璃进行亚类划分，并分析分类模型的合理性和敏感性。首先对两种玻璃的化学成分进行数值统计，对这些化学成分的标志性数据进行分类汇总，以此作为亚类划分的依据；然后在此基础上进行亚类划分，这里采用聚类分析模型解决该问题，观察玻璃的这些化学成分在风化前后的变化，以及玻璃的纹饰，颜色等的变化，通过聚类分析进行分类；最后进行数据的灵敏度检验，也就是灵敏度分析，得出检验结果，给出相应的合理性依据。

2.1.3问题三的分析

在解决问题三时，根据问题要求，我们需要对未知玻璃文物的主要化学成分进行分析，鉴别这些文物属于什么类型的，并进行敏感性分析，即灵敏度分析。首先将表单3中的数据进行分类，然后结合第二问的结论，对这些玻璃进行分类鉴别，初步得出这些玻璃的类别；接着对表格1和表格2中的数据进行回归分析，得到不同类型的玻璃各化学成分所占的权重，最后将得到的结果进行灵敏度分析，即可以对这些化学成分中的某一个添加数据扰动，代入原问题二模型中，来验证模型的稳定性和敏感性。

2.1.4问题四的分析

在解决问题四时，根据问题要求，要分析不同玻璃化学成分之间的关系，然后进行不同玻璃化学成分之间的差异性。在分析化学成分的关系时，运用第三问的回归分析模型，问题三中已得到不同玻璃化学成分所占的比重，我们将这些比重数据进行数据可视化，分析其中的关系；在解决不同玻璃化学成分之间的差异性时，根据问题一所建立的单因素方差分析模型，以玻璃类型作为因素，相对应的化学成分作为因素所处的水平，计算出不同玻璃类型对应的p值，进而分析差异性。

三、模型的假设

为了便于模型求解，现做如下假设：

1. 忽略各种玻璃类型在检测过程中可能受到的风化影响；
2. 检测玻璃化学成分含量时的环境一致；
3.  假定实验数据在进行处理的过程中没有人为的操作误差；
4. 假定在试验过程中除因素自身外其他影响指标的因素都保持不变；
5. 在进行亚类划分时不考虑无效数据组。
6. 

五、模型的建立与求解

5.1 问题一的求解 首先根据题目的要求将题目分为三个小问。第一小问是需要对玻璃表面风化情况 与玻璃类型、纹饰和颜色的关系进行分析，即进行差异显著性分析，经过问题的研究 分析可以运用单因素方差分析模型，以表面是否分化作为因素，以玻璃类型、纹饰和 颜色的情况作为因素所处的水平，利用工具分析计算出 p 值，观察分析显著性 p 值的 值是否小于 0.05，进而可以分析出其中玻璃表面风化情况与玻璃类型、纹饰和颜色之 间的关系，对于第二小问是结合玻璃的类型分析化学成分含量的变化规律，对于这一 问题，我们可以进行描述性统计分析，计算相关数值的平均数，中位数和标准差，进 行折线图分析，正态分布检验等的研究计算，总结其中的规律，最后一个小问是预测 风化前的化学成分含量，根据前面已经得到的数据变化规律，总结出各个化学成分的 变化情况，这是可以找出来相关的函数关系，从而进行风花前的数据的预测。

5.1.1 数据预处理 为了方便数据的分析与计算，首先对数据进行预处理，从而更加简便而且更加直 观的分析出玻璃文物的表面风化与其玻璃类型、纹饰和颜色等之间的关系，也更加便 于通过结合玻璃的类型，分析出文物样品表面有无风化化学成分含量之间的统计规律， 那么同样的也可以更加简洁的根据相关函数进行风化点检测数据，从而预测出其风化 前的化学成分含量，所以先对数据进行预处理，根据题目中的要求将成分将成分比例 累加和介于 85%~105%之间的数据视为有效数据，经过对数据的成分比例和进行计算， 计算结果的前 24 组数据如下表，其它数据在附录查看。

         Step 1:第一步首先对数据进行预处理操作，根据题目中要求：将成分将成分比例 累加和介于 85%~105%之间的数据视为有效数据，从制出的数据处理表格可以知道， 表单 2 中的数据，编号 15 和编号 17 的总成分小于 85%，所以很明显，数据为无效数 据，所以剔除掉这些数据。

        Step 2:通过观察表单 1 可知，颜色列中的数据有四个数据为空，此处我们为了处 理数据方便，在研究表面风化与颜色的关系时，将这四组数据视为无效数据，进行剔 除。

        Step 3:通过观察表单 2 可以知道，化学检测成分有很多都是空缺，根据题中的要 求可以知道这些是未检测到该成分，此处我们为了处理数据更加方便，将这些空白处 的数据进行补“0”操作，然后进行接下来的计算。

5.1.2 问题一第(1)问模型的建立与求解 我们在研究玻璃表面风化与其玻璃类型、纹饰和颜色的关系时，显然需要将表面 分化作为因变量，将其他变化量作为自变量进行分析，这个地方采用了单因素方差析模型，以表面分化作为因素，将表面分化进行一定数据转化作为因素所处的水平， 颜色、服饰和类型作为检验指标，以此来检验表面分化与玻璃类型、纹饰和颜色的关 系。 我们建立了单因素方差分析模型，首先方差分析是研究一个（或多个）自变量对 一个（或多个）因变量影响的方法。若控制变量对观察变量产生了显著影响，可进一 步分析出该控制变量的某一水平对观察变量产生了显著影响[2]

想要全部的代码以及论文可以看我的github  。GitHub – Ggy-king/2022-Mathematical-Modeling-National-Code: 2022年数学建模国赛代码–matlab 其实也用了部分python 主要是用它做线性回归预测 再没其他的了 剩下还是偏向于matlab的的实用性

或者评论留下邮箱 私发给你。

代码一：

代码1：
# 使用bp用神经网络算法预测风化前后变化值

import numpy as gg
from numpy import dot, random, ones_like, exp, ones, zeros, multiply, zeros_like
import matplotlib.pyplot as plt
import math

# 若是10个作为数据 先以铅钡二氧化硅为例 使用风化后预测风化前
# 三类数据输入，每类取7个作为训练数据，3个作为测试数据 先是测试模拟
kind1 = gg.array([[51.26, 3.16, 2.87], [19.79, 2.32, -5.8],
                  [39.57, 2.18, -3.39], [35.78, 1.21, -4.73],
                  [20.14, 1.58, -4.78], [33.59, 1.01, -3.63],
                  [25.42, 1.40, -1.89], [29.15, 1.44, -3.22],
                  [17.98, 1.33, -4.38], [12.38, 1.33, -4.38]], dtype=float).reshape(-1, 3)

label1 = gg.zeros_like(kind1)
label1[:, 0] = ones([len(label1)], dtype=float)
kind1 = gg.hstack((kind1, label1))
ext = gg.ones(len(kind1))
ext = ext.reshape(10, -1)
kind1 = gg.hstack((ext, kind1))
kind2 = gg.array([[-0.24, 0.93, -1.01], [-1.18, 0.39, -0.39], [0.74, 0.96, -1.16],
                  [-0.38, 1.94, -0.48], [0.02, 0.72, -0.17], [0.44, 1.31, -0.14],
                  [0.21, 0.03, -2.21], [0.37, 0.28, -1.8], [0.18, 1.22, 0.16],
                  [0.46, 1.49, 0.68]]).reshape(-1, 3)

label2 = zeros_like(kind2)
label2[:, 1] = ones([len(label2)], dtype=float)
kind2 = gg.hstack((ext, kind2, label2))
kind3 = gg.array([[31.94, 1.17, 0.64], [36.93, 3.45, -1.33], [55.21, 0.99, 2.69],
                  [34.34, 3.19, 1.51], [44.26, 1.79, -0.87], [22.65, -0.22, -1.39],
                  [18.96, -0.44, -0.92], [32.65, 0.83, 1.97], [33.44, 0.68, -0.99],
                  [48.12, -0.45, 0.08]]).reshape(-1, 3)

label3 = zeros_like(kind3)
label3[:, 2] = ones([len(label3)], dtype=float)
kind3 = gg.hstack((ext, kind3, label3))
#all_kind = gg.vstack((kind1, kind2, kind3))
train_mate = gg.vstack((kind1[:7], kind2[:7], kind3[:7]))
test_mate = gg.vstack((kind1[7:], kind2[7:], kind3[7:]))

# 函数及导数预测
def ta_h(x):
    return math.tah(x)
def diff_tag_h(x):
    return 1.0 / (1 + pow(x, 2))

# sigmoid
def sigmoid(x):
    return 1.0 / (1 + exp(-x))

# sigmoid求导
def diff_sigmoid(x):
    out = sigmoid(x)
    return out * (1 - out)

代码二： 

# 线性函数
def linear(x):
    return x

# 线性函数求导数
def diff_linear(x):
    return ones_like(x)  # 对线性函数求导

# 定义神经网络 模型
class Gao:
    # 输入层、隐藏层、输出层的节点（数）
    def __init__(self, n_i, n_h, n_o):
        # 获取各层节点数量
        self.n_i = n_i
        self.n_h = n_h
        self.n_o = n_o
        
# 激活神经网络上全部节点 存储加权求和之后 对应 net
        self.mate_i = ones(self.n_i)
        self.mate_net_h = ones(self.n_h)
        self.mate_net_o = ones(self.n_o)
       
 # 对应风化前后的数据操作
        self.mate_y = ones(self.n_h)
        self.mate_z = ones(self.n_o)
        self.f0_net_k = ones(self.n_o)
        self.delta_k = ones(self.n_o)
      
  # 初始化权重矩阵
        self.wi = random.random((self.n_h, self.n_i))
        self.wo = random.random((self.n_h, self.n_o))
       
 # 待更新缓存
        self.delta_wi_temp = self.wi
        self.delta_wo_temp = self.wo
    def calculate_output(self, iggut):
        # iggut layer
        self.mate_i = iggut
        # in - hidden
        self.mate_net_h = dot(self.wi, self.mate_i)
        self.mate_y = gg.array(list(map(ta_h, self.mate_net_h)))
        # self.mate_y = self.mate_y.reshape(-1, 1)
        # hidden-output
        self.mate_net_o = dot(self.mate_y, self.wo)
        self.mate_z = list(map(sigmoid, self.mate_net_o))
        return self.mate_z  # 输出
    def BPBP(self, target, upmate_flag, rate_1, rate_2):
        # 得到误差
        error_t_k = target - self.mate_z
        for i in range(self.n_o):
        self.f0_net_k[i] = diff_sigmoid(self.mate_net_o[i])
        self.delta_k = gg.multiply(self.f0_net_k, error_t_k)
        mate_y_temp = self.mate_y.reshape(-1, 1)
        delta_wo = dot(mate_y_temp, self.delta_k.reshape(1, 3))
        epsilon = zeros(self.n_h).reshape(-1, 1)
        for i in range(self.n_h):
            epsilon[i] = multiply(self.delta_k, self.wo[i:i + 1][0]).sum()
        delta_wi = rate_2 * dot(epsilon, self.mate_i.reshape(1, -1))
        self.delta_wo_temp = self.delta_wo_temp + delta_wo
        self.delta_wi_temp = self.delta_wi_temp + delta_wi
        if upmate_flag == 1:
            # 测试即
            self.wo = self.wo + rate_2 * delta_wo
            # 测试即
            self.wi = self.wi + rate_1 * delta_wi
        error = 0.5 * dot((target - self.mate_z),
                          (target - self.mate_z).reshape(-1, 1))
        return error
    def train(self, patterns, iggut_mate, rate_1, rate_2):  # 全部样本
        # stop_flag = 0
        error_set = []
        acc_set = []
        step = 0
        sample_len = len(patterns)
        sample_num = 0
        rate_temp = 0
        # while stop_flag == 0:
        for m in range(5000):
            step += 1
            upmate_flag = 1
            for p in patterns:
                sample_num += 1
                igguts = p[1:4].reshape(-1, 1)
                targets = p[4:]
                if sample_num == sample_len:
                    upmate_flag = 1
                self.calculate_output(igguts)
                error = self.BPBP(targets, upmate_flag, rate_1, rate_2)
            rate = self.test(iggut_mate)
            rate_temp = rate_temp + rate
            if step % 100 == 0:
                error_set.append(error)
                print("error", error, "acc:", rate)
            if step % 10 == 0:
                rate_temp = rate_temp / 10
                acc_set.append(rate_temp)
                rate_temp = 0
        return error_set, acc_set
#     def test(self, iggut_mate):
#       
  #  此处为测试
#         ok = 1
#         for p in iggut_mate:
#             igguts = p[1:4].reshape(-1, 1)
#             targets = p[4:]
#             output = self.calculate_output(igguts)
#             out_kind = gg.where(output == gg.max(output))

#             if targets[out_kind] == 1:
#                 ok = ok + 1
#         rate = ok / len(iggut_mate)
#         return rate
#     def plot_plot(self, error_set0, error_set1, error_set2):
#         set_len = len(error_set1)
#         plt.plot(range(set_len), error_set0, range(set_len),
#                  error_set1, '-', range(set_len), error_set2, '--')
#         plt.legend(['error_set0', 'error_set1', 'error_set2'], loc='best')
#         plt.title("ErrorCurve")
#         plt.show()
# def Run(test_mate=test_mate):
#     pat = train_mate
#     test_mate = test_mate
#     rate_1 = 0.001
#     rate_2 = 0.004
#     rate_3 = 0.006
#     # 创建一个神经网络：输入层  隐藏层  输出层
#     n0 = Gao(3, 3, 3)
#     error_set0, acc0 = n0.train(pat, test_mate, rate_2, rate_2)
#     n1 = Gao(3, 6, 3)
#     error_set1, acc1 = n1.train(pat, test_mate, rate_2, rate_2)
#     n2 = Gao(3, 8, 3)
#     error_set2, acc2 = n1.train(pat, test_mate, rate_2, rate_2)
#     n2.plot_plot(error_set0, error_set1, error_set2)
# if __name__ == '__main__':
#     Run()

代码三： 

代码2：图10和图11的MATLAB代码
%%定义函数 数据另导
function [fitresult,gof] = createFit6(x,y,z)

%自拟合
[xDate,yDate,zDate] = prepareSurfaceDate(x,y,z);
ft = fittype('poly24')
[fitresult,gof] = fit([xDate,yDate],zDate,ft);
figure('Name','untitled fit 1');
%设量
h = plot(firesult,[xDate,yDate],zDate);

%备注暂且不标
% legend()
xlable('x1','Interpreter','none');
ylable('y1','Interpreter','none');
zlable('z1','Interpreter','none');
grid on
%导出
view(-28.7,39.5);

代码3：第二问第一小问的MATLAB代码
%%  读取处理后的数据集
data_result = xlsread('高钾铅钡分类.xlsx');
M = size(data_result, 1);

%%  添加截距项
L_train = data_result(:, 1: end - 1);
R_train = data_result(:, end);

%%  建立模型
ctree = fitctree(L_train, R_train, 'MinLeafSize', 8);

%%  查看决策树视图
view(ctree, 'mode', 'graph');

%%  仿真测试
T_sim1 = predict(ctree, L_train);
%%  计算准确率
error1 = sum((T_sim1 == R_train)) / M * 100 ;

%%  绘图
figure
plot(1: M, R_train, 'r-*', 1: M, T_sim1, 'b-o', 'LineWidth', 1)
legend('真实值', '预测值')
xlabel('预测样本')
ylabel('预测结果')
string = {'预测结果对比'; ['准确率=' num2str(error1) '%']};
title(string)
xlim([1, M])
grid

%%  混淆矩阵
figure
cm = confusionchart(R_train, T_sim1);
cm.Title = 'Confusion Matrix for Data';
cm.ColumnSummary = 'column-normalized';
cm.RowSummary = 'row-normalized';

想要全部的代码以及论文可以看我的github-网址如下(可按需下载):GitHub – Ggy-king/2022-Mathematical-Modeling-National-Code: 2022年数学建模国赛代码–matlab 其实也用了部分python 主要是用它做线性回归预测 再没其他的了 剩下还是偏向于matlab的的实用性

或者评论留下邮箱 私发给你。