基于多元线性回归的古代玻璃体系判别方法

钱鸿羽，岳思情，许 雷，郑心怡

（内江师范学院，四川 内江 641000）

玻璃作为人类最早期发明的人造材料之一[1]，在世界各国之间的贸易与交流活动中占据重要位置。由于其相对稳定的化学性质，迄今为止保存了大量古代玻璃制品，这是研究国家间早期贸易的宝贵物证。对古代玻璃制品的探究，对进一步研究各国间的交流活动具有重要意义。当前，学术界对古玻璃的研究主要集中在古代玻璃的成分体系研究、玻璃起源与产地研究、玻璃制作工艺研究等领域[2]，而对古代玻璃化学成分体系的研究为进一步研究玻璃生产技术及其起源奠定了基础，因而成为古代玻璃研究的重要一环。目前研究古代玻璃化学成分的方法主要有张欣睿[3]采用的X射线荧光分析、李玲等[4]采用的PIXE 技术、李颖等[5]运用拉曼光谱、红外光谱和ED-XRF 等技术，以及李曼[6]运用激光剥蚀电感耦合等离子体发射光谱进行化学成分分析的方法。由于古代玻璃易受外界影响而产生风化，从而导致内部化学成分发生变化，影响对古代玻璃化学成分体系的判别。在2022 年全国大学生数学建模竞赛C题的基础上，引入风化前玻璃各化学成分预测函数，通过风化后古代玻璃的化学组成含量预测相应的风化前化学成分，从而确定玻璃样品的成分。

1 理论基础

1.1 加性对数比变换（ALR）

加性对数比变换（ALR）是目前在成分分析中运用较为普遍的数据转换方法。为了消除定和约束的影响，通过计算各个组别所占成分比例与其均值的比值的对数，将成分数据从单形空间变换到欧式空间，从而能够应用经典统计方法解决成分数据问题。

其具体计算公式如下：

式中：i表示数据的组别，xi表示各组数据的成分比例，xp表示数据比例的平均数，yi则为新产生的数据变量。

1.2 独立样本Mann-Whitney U检验

Mann-Whitney U检验是用得最广泛的两独立样本秩和检验方法，其假设基础是：若两个样本有差异，则它们的中心位置将不同。具体步骤如图1 所示。

图1 独立样本M ann-W hitneyU 检验流程图

1.3 最小二乘法

定义2 个变量的n对样本数据为（x1,y1），（x2,y2），…，（xn,yn）设回归直线为：y=a+bx。

由yi=bxi+a+ei得代表样本点（xi,yi）到直线上点（xi,bxi+a）的距离。

当 |ei|=0 时，则表示点（xi,yi）在直线上，利用来刻画各样本数据与直线y=a+bx的“整体接近程度”；为了便于计算，通常利用来刻画“整体接近程度”。

2 风化前化学成分预测模型

由于古代玻璃易受风化作用影响而改变其各化学成分含量，因此根据样本数据建立以下风化前化学成分预测模型。

2.1 选取重要性指标

根据干福熹研究所得，我国早期古代玻璃体系可分为铅钡玻璃、高钾玻璃等[7]，故可按玻璃体系分别建立玻璃风化程度函数与风化前化学成分预测函数。

定义玻璃的风化程度函数为：

其中a为常量，Ci为重要性指标在玻璃制品中的含量比例，bi为重要性指标的对应系数，n为重要性指标的数量。

将样本数据进行预处理，将无法检测出化学成分的部分填充为0，去除成分比例累加总和处于85%～105%区间范围外的数据，以便于进行后续数据分析。

由于n元成分数据受到定和约束，故必然存在一组数据与其余数据呈现负相关情况，易与实际情况产生偏差，故针对普通数据的传统统计学分析方法难以适用于成分数据分析。为了进一步对数据进行分析，对数据进行加性对数比变换（ALR）处理，将原始数据比例关系转化为加性关系，从而可以更加充分的体现成分特性，使成分数据中的可解释性更强。

根据玻璃体系划分样品数据，提取处理后数据中玻璃体系属于高钾的样本，并对化学成分含量进行正态性检验。发现总体不符合正态分布。因此，选择非参数检验中的独立样本Mann-Whitney U 检验[8]，按照Cohen′s d 系数排序分析具有显著性差异的化学成分，从而得到文物样品表面是否存在风化与化学成分含量之间的统计规律，确定了判断玻璃风化程度的重要指标，检验结果见表1。

表1 M ann-W hitneyU 检验分析结果表

2.2 建立多元线性回归模型

对Cohen′s d系数进行排序，选取排名靠前的部分化学成分，认为它们有显著差异，并将其作为函数的自变量。高钾玻璃风化程度函数的自变量为玻璃中二氧化硅、氧化钾、氧化钙、氧化铝的对应含量比例。

根据玻璃样本大体的风化情况对风化程度进行赋值。根据多次拟合结果，确定未风化样品的值为0，一般风化样品赋值为0.8。使用最小二乘法进行函数拟合，以确定相应化学成分的系数与未知常数。对获得的残差进行检验，去除异常值，并进行多次线性拟合，直至总体残差的置信区间均包含零点，表明回归模型可以更好地符合原始数据，从而得到最优的多元线性回归函数。

拟合所得的高钾玻璃风化程度函数为：

其中YK代表高钾玻璃的风化比例，CSiO2代表SiO2在玻璃制品中的含量比例，CK2O，CCaO，CAl2O3同理。函数拟合程度如图2所示。

图2 高钾玻璃风化程度函数拟合结果图

同理，其余体系的玻璃风化程度函数也可按照上述方法得出。

通过对铅钡玻璃样本中的数据进行独立样本Mann-Whitney U 检验，按照Cohen′s d 系数进行排序，可得铅钡玻璃风化程度函数的自变量为玻璃中二氧化硅、氧化钙、氧化铝、氧化铅、五氧化二磷的对应含量比例。根据多次拟合结果，确定将未风化样本赋值为0，一般风化样品为0.75，严重风化样品赋值为1。经过多次拟合与剔除异常值，拟合得到的铅钡玻璃风化程度函数为：

其中YPb代表高钾玻璃的风化比例，代表SiO2在玻璃制品中的含量比例同理。函数拟合程度如图3所示。

图3 铅钡玻璃风化程度函数拟合结果图

计算函数预测结果与实际结果的差值，得到预测数据的偏差情况，所得结果如图4所示。

图4 样本预测数据偏差情况图

2.3 风化前化学成分预测模型

处理后的样品数据应根据玻璃系统进行分类，并根据风化情况进行二次分类。去除异常数据后，计算不同风化情况下各化学成分自身的差值。由于在化学反应中涉及的化学成分是成比例的，因此可以通过简化得到不同玻璃体系下风化过程中各化学成分反应比例。高钾玻璃与铅钡玻璃风化反应比例分别为：

分别探究玻璃在不同玻璃体系下的风化程度与变化程度最大的化学成分含量之间的线性关系，得到高钾玻璃和铅钡玻璃的风化程度与二氧化硅含量的函数分别为：

通过反解该线性函数，可以得到相应风化程度下参与反应的二氧化硅含量。将获得的数据乘以各化学成分反应比例，即可得对应样本在风化过程中参与反应的各化学成分含量。风化前的数据可以通过从风化后的数据中减去变化来获得。

高钾玻璃与铅钡玻璃各化学成分含量计算公式分别为：

其中Ci为风化后对应化学成分含量，为所求风化前对应化学成分含量，qi为对应反应比例。

将预测函数所得的数据代入拟合的风化程度函数，发现5.1%以内的计算偏差比例为94.12%，其中严重风化情况下的样本计算偏差相对较大，故认为预测函数较为有效。

3 玻璃化学成分体系判别模型

3.1 选取重要性指标

为避免风化引起的化学成分含量变化对判定玻璃体系的影响，本次从样本数据中选取无风化样本作为数据拟合的内容。筛选后，对类别进行赋值处理，将“-1”赋值给高钾玻璃，将“1”赋值给铅钡玻璃，使得定类变量“类别”能够转化为定量变量，便于后续判断和筛选对玻璃体系影响较大的化学成分。

通过对化学成分进行正态性检验，发现大部分数据并不满足正态分布。故同样采用独立样本Mann-Whitney U 检验，计算各个化学成分的Cohen′s d 系数并对其进行排名，所得结果见表2。

表2 M ann-W hitneyU 检验分析结果表

3.2 建立玻璃体系判别模型

选取Cohen′s d 系数大于1 的化学成分，这些化学成分被认为具有较大显著性差异。因此，可以得出结论，对玻璃体系的判别能力具有显著性差异从大到小分别为氧化钾、氧化铅、氧化钙、氧化钡，氧化硅和氧化锶，它们在玻璃中相应化学成分的比例即为判别函数的自变量。

将量化后的玻璃体系作为函数因变量，用最小二乘法求解相应自变量的系数。经过多次拟合与异常值处理，可获得玻璃体系判别函数为：

函数拟合程度如图5所示。

图5 玻璃体系判别函数拟合结果图

选取部分样本数据代入判别函数，经过多次测试，发现当判断偏差设置为±0.4时，该测试的准确度约为93.75%，因此认为多元线性回归所得的判别函数具有较高的准确性。

将预测结果与实际数据的偏差情况可视化，其结果如图6所示。

图6 判别函数预测结果偏差图

4 玻璃类型判别方法

古代玻璃易受到外界影响而风化，从而导致内部化学成分含量发生较大变化，故难以直接通过玻璃化学成分体系判别模型得出结果，所以首先对该玻璃样本风化前的化学成分进行预测。根据探究得到高钾玻璃和铅钡玻璃的风化程度与二氧化硅含量的关系分别如公式（5）、（6）所示，即铅钡玻璃与高钾玻璃在风化后其二氧化硅和其他化学成分含量的变化趋势相反，故若高钾玻璃被判别为铅钡玻璃代入其风化前化学成分预测模型时，会与其风化前数据产生更大的差异。此时代入玻璃体系判别后所得结果必然与假设相矛盾。综上可得，当且仅当对玻璃类型的假设正确时，代入玻璃体系判别模型后所得的结果才与假设相一致。

因此，可以先对玻璃样品类型进行假设。首先，假设所有风化玻璃均为高钾玻璃。将数据代入高钾玻璃风化前后化学元素含量的变化函数，计算风化前各化学成分的含量。将获得的数据代入玻璃体系判别函数，判断所得结果是否在可接受范围内。如果是，则认为该假设成立，文物属于高钾玻璃；如果不是，则假设与结果相矛盾，假设不成立，再次假设该文物属于铅钡玻璃，重复上述步骤。

若两者皆不成立，则将该玻璃样品进行单独讨论，计算结果与判定标准之间的差距，以差距较小的一方作为判别结果。

5 实例计算

基于2022 年全国大学生数学建模竞赛C 题所提供的数据，对表单中未知体系的风化玻璃文物的化学成分进行分析，鉴别其所属体系，具体数据见表3。

表3 玻璃样品及其化学成分含量表

将提供的玻璃样本均假设为高钾玻璃，并代入对应的风化前各化学成分预测函数，得到风化前各化学成分的含量。将所得数据代入玻璃体系判别函数，观察是否与假设相一致。若不一致，则假设样本的体例为铅钡玻璃并重复以上步骤。所得数据与判断结果见表4。

表4 风化前化学成分含量预测与分类结果表

从表4中观察所得，发现编号为08与26的两份样品在两种假设下的偏差均高于设定标准，通过计算并比较其在不同玻璃体系情况下的偏差程度，并取较小一方作为判断结果，可得编号08 与26 的玻璃体系为铅钡玻璃。

6 总结

研究玻璃制品的化学成分及鉴别对保护古代玻璃文物具有重要意义。本文从成分数据视角出发，采用加性对数比变换以消除定和限制对后续分析的影响。引入了风化前玻璃各化学成分预测函数，并通过风化后古玻璃的化学组成含量预测相应的风化前化学成分含量，从而确定玻璃样品的成分体系，得到古代玻璃体系判别模型。在应用本模型时，先对玻璃样品的类别进行假设，在假设的基础上代入相应的风化前化学成分预测模型，得到在该假设下风化前各化学成分比例。随后将此数据代入古代玻璃体系判别模型，根据所得结果与假设的一致性确定玻璃样品所属的玻璃类别。模型拟合过程中所运用的多元线性回归方法降低了计算的复杂性，同时确保了结果的准确性。最后，通过实例说明了该方法的有效性和适用性，为古代玻璃化学成分体系的分析与判别提供了新的思路。