简介:面向2022年全国大学生数学建模竞赛C题备赛者,这份PDF完整收录了古代玻璃文物成分分析与鉴别问题的解题方案与配套代码说明,聚焦数学建模在考古数据场景下的实际应用。内容围绕数据预处理、主成分分析、因子分析、敏感性分析及灰色关联度分析等核心方法展开,系统覆盖了从数据探索、风化系数定义、亚类划分、分类预测到不同玻璃类型化学成分关联性比较的全流程建模思路,并附有关键代码与结果解释,便于读者对照理解。资源共1个文件,大小3.55MB,PDF格式方便在线阅读与打印。目前已有927人学习下载,其问题拆解和模型选择逻辑具有较高参考价值,尤其适合备赛数学建模、需要参考完整赛题分析框架的高校学生与建模爱好者。读者可快速掌握风化数据场景下的降维、因子转换、鲁棒性检验与灰色关联度建模技巧,并将其迁移到其他成分分析类问题中。
1. 风化系数:为什么风化效应会让成分数据说反话
古代玻璃制品的成分分析,最反直觉的地方在于:风化并不是简单地让某些成分流失,而是让整组数据都在说谎。同样是埋藏环境下,高钾玻璃风化后二氧化硅占比显著升高、氧化钾显著降低;铅钡玻璃却完全反过来,二氧化硅占比显著降低、氧化铅和氧化钡显著升高。如果拿着风化后的检测数据直接做分类或回归,模型学到的很可能是"风化程度",而不是"玻璃类型"本身。这份 2022 年数学建模 C 题的标准解法,核心思路就是先定义风化系数把风化程度量化,再用主成分分析降维、因子分析反推未风化成分、灰色关联度分析化学成分关联——这一整套链路解决的不只是考古问题,更是"观测值被环境篡改后如何建模"的通用范式。
2. 数据探索先行:三个反直觉事实决定建模方向
建模比赛里最常见的翻车点不是模型不够高级,而是没看数据就套模型。这道题最值得借鉴的地方,恰恰是把数据探索放在了建模之前,而且探索结论直接决定了后续每个模型的选择。
2.1 两类玻璃的风化方向截然相反
玻璃的主要成分是二氧化硅,但为了降低熔融温度,古人会添加助熔剂。加入铅矿石、钡矿石得到的是铅钡玻璃,加入草木灰得到的是高钾玻璃。这些助熔剂残留物成为玻璃的次要成分:氧化铅(PbO)、氧化钡(BaO)、氧化钾(K2O)、氧化钙(CaO)、氧化镁(MgO)、氧化铝(Al2O3)等。
风化效应是玻璃与埋藏环境长期进行化学元素交换的结果。关键在于,两类玻璃的风化产物方向完全相反:
| 成分变化方向 | 铅钡玻璃 | 高钾玻璃 |
|---|---|---|
| 二氧化硅(SiO2) | 显著降低 | 显著升高 |
| 氧化铅 / 氧化钡 | 显著升高 | 变化不显著 |
| 氧化钾 | 变化不显著 | 显著降低 |
这个结论来自文献中对铅钡玻璃、钾玻璃风化机理的说明,也在附件表单中得到了验证。高钾玻璃风化后,表面富硅层形成,二氧化硅相对占比被抬高;铅钡玻璃风化后,铅、钡成分相对富集,二氧化硅被"稀释"。同一套风化机制,作用在两类玻璃上产生了方向相反的表观变化。
这带来一个直接后果:如果拿风化后的原始成分数据做聚类,高钾风化样本和铅钡未风化样本可能在某些维度上靠得很近,因为高钾风化后 SiO2 升高、K2O 降低,铅钡未风化时恰好也是 SiO2 高、铅钡低。数据探索的价值在这里就体现出来了——它告诉你,必须先给"风化程度"一个定量指标,否则后续所有分类都建立在流沙上。
2.2 风化系数的定义逻辑与阈值判定
论文里给出了两个风化系数定义,这是整个建模的地基。
对铅钡玻璃,风化系数定义为:
# 铅钡玻璃:风化导致 PbO、BaO 升高,SiO2 降低 # 所以比值越大,风化越严重 lead_barium_weathering = (PbO + BaO) / SiO2对高钾玻璃,风化系数定义为:
# 高钾玻璃:风化导致 SiO2 升高,K2O、CaO、MgO、Al2O3、Fe2O3、CuO 降低 # 所以 SiO2 在分子上,其他成分在分母上 potassium_weathering = SiO2 / (K2O + CaO + MgO + Al2O3 + Fe2O3 + CuO)这两个定义不是拍脑袋想出来的,而是严格对应风化方向:风化系数应该和风化程度单调相关。铅钡玻璃风化越重,PbO、BaO 越多,SiO2 越少,比值自然越大;高钾玻璃风化越重,SiO2 越多,其他成分越少,比值也越大。
有了风化系数之后,阈值判定就有了依据。观察数据分布可以发现一个典型现象:铅钡玻璃-A 纹类型的风化系数在 0.879 和 1.3 之间存在明显突变,中间几乎没有过渡值;高钾玻璃的风化系数在 7.23 和 12.41 之间也存在跳跃。这种"断层"说明风化样本和未风化样本在定量指标上是可分的,于是把阈值分别定在 1.0 和 10.0。
这里有个容易被忽略的细节:阈值不是从统计检验里算出来的,而是直接观察数据突变得到的。这个做法在数学建模里是合理的——你是在找一个能把样本分成两堆的分界点,而不是在估计一个总体参数。不过,论文里也承认存在少量标签不一致的情况,比如某样本的标注与风化系数判别的结果矛盾,最后以风化系数的判定为准做了手工调整。这种做法在比赛里无可厚非,但放到真实考古场景中,就意味着你需要第三方证据来验证"标注错了"而不是"风化系数定义错了"。
2.3 数据探索发现的其他规律
除了风化方向相反之外,数据探索还发现了三个具体现象:
第一,纹饰 B 的玻璃文物全部是高钾玻璃、全部风化、颜色全部是蓝绿色。这个组合的高一致性说明纹饰和产地、时代高度相关。第二,纹饰 A 的高钾玻璃全部无风化,而铅钡玻璃风化与未风化比例约为 11:5。第三,纹饰 C 的高钾玻璃全部无风化,铅钡玻璃风化与未风化比例为 17:73。
把这三个现象放在一起,可以得到两个推论:同等条件下铅钡玻璃比高钾玻璃更容易风化;纹饰 B 所代表的环境对玻璃有更强的风化效应,因为连相对耐风化的高钾玻璃都全部风化了。
这些结论直接在论文中作为问题一的回答基础,也为问题二的分类模型提供了先验知识:分类时必须把"是否风化"作为分层变量,先区分风化与未风化,再做类型判别。这就是为什么后面的主成分分析要分四类数据分别做——高钾风化、高钾未风化、铅钡风化、铅钡未风化——而不是把所有样本丢进一个 PCA 里。如果混在一起降维,主成分会被风化差异主导,类型差异反而被淹没。
3. 主成分与因子分析:把风化后的成分反推回未风化
题目最硬核的问题一是预测风化前的成分含量。这本质上是一个"反演"问题:你只能测到风化后的成分,但想还原风化前的成分。这个问题的解决思路非常巧妙——用因子分析建模风化过程,再对负载矩阵求逆反推。
3.1 为什么必须先做主成分分析降维
先看数据维度。表单 2 中每个样本有 14 个化学成分变量,而样本量只有几十个。直接用 14 维数据做因子分析,会遇到两个问题:一是运算量偏大,二是很多成分之间高度相关,存在多重共线性,直接做回归或因子分析,结果会被干扰成分带偏。
主成分分析的作用就是把 14 个原始变量压缩成少数几个互不相关的主成分。每个主成分是原始变量的线性组合,第一个主成分承载最大方差,第二个次之,依此类推。论文里的具体做法是,按玻璃类型和风化状态分成四类数据,对每一类分别做 PCA,选取累计贡献率达到 85% 以上的主成分。
以铅钡未风化类型为例,前 6 个主成分的累计贡献率达到 88.89%,对应的主要化学成分是 SiO2、Na2O、K2O、CaO、MgO、Al2O3。高钾未风化类型的前 5 个成分贡献率同样达到 88.89%。这意味着原始 14 维信息被压缩到 5-8 维,同时保留了绝大部分方差。
PCA 的具体实现并不复杂:
import numpy as np def pca_components(X, n_components=None, threshold=0.85): # X: 样本矩阵,每行一个样本,每列一个化学成分 # 第一步:标准化 X_mean = X.mean(axis=0) X_std = X.std(axis=0) X_norm = (X - X_mean) / X_std # 第二步:计算相关系数矩阵 # 注意这里用的是相关系数矩阵而非协方差矩阵, # 因为不同化学成分量纲相同但方差差异大,标准化后更公平 R = np.corrcoef(X_norm, rowvar=False) # 第三步:特征值分解 eigenvalues, eigenvectors = np.linalg.eigh(R) idx = np.argsort(eigenvalues)[::-1] eigenvalues = eigenvalues[idx] eigenvectors = eigenvectors[:, idx] # 第四步:按累计贡献率选主成分个数 total_var = eigenvalues.sum() cum_ratio = np.cumsum(eigenvalues) / total_var if n_components is None: n_components = np.searchsorted(cum_ratio, threshold) + 1 # 主成分得分:标准化后的数据投影到特征向量上 scores = X_norm @ eigenvectors[:, :n_components] return scores, eigenvectors[:, :n_components], cum_ratio这里有个关键细节:做 PCA 之前必须先标准化。如果直接对原始成分数据做特征值分解,SiO2 这种含量经常在 60% 以上的变量会主导第一主成分,而 PbO 这种含量波动大的变量却因为数值小被忽略。标准化之后,每个变量在降维决策中的权重才真正由其相关性决定。
3.2 因子分析模型的核心逻辑
因子分析假设可观测的化学成分向量 x 由公共因子 f 线性组合而成:
x = μ + A f + ε其中 μ 是成分期望向量,f 是公共因子向量(这里对应未风化时的成分特征),A 是负载矩阵,ε 是特殊因子即随机误差。
这个模型的物理含义是:风化的过程是连续的、由多个公共因子驱动的,某个时刻观测到的成分比例,是原始成分经过风化演化后在当前时刻的表现。公共因子 f 数量小于原始变量数量,这就是降维的本质。
关键的建模步骤在于如何获取负载矩阵 A。论文中的做法是:对标准化后数据的相关系数矩阵做特征值分解,取前 m 个特征值 λ 和对应的特征向量 u,构造:
A = [sqrt(λ1)·u1, sqrt(λ2)·u2, ..., sqrt(λm)·um]也就是把特征值开方后乘到特征向量上。这个 A 的列就对应各个公共因子在各个原始变量上的载荷。
获取公共因子的表达式需要对式 (3) 做变换。把 x = μ + A f + ε 稍作整理,得到 f 关于 x 的线性表示:
f = A^{-1} (x - μ)这就是"由风化后成分反推未风化成分"的转换公式。论文中给出的是铅钡玻璃-C 纹类型的转换系数矩阵——一个 8 行 10 列的系数表,每一行对应一个公共因子,每一列对应一个原始化学成分。把风化后的实测成分代入,就能得到未风化时的预测值。
用代码来演示这个反推过程会非常直观:
import numpy as np def factor_reverse(observed, mu, A): """ 由风化后成分反推未风化成分 observed: 风化后检测到的化学成分向量 mu: 训练集成分均值向量 A: 负载矩阵,shape (n_components, n_features) """ # 中心化:减去均值 centered = observed - mu # 求负载矩阵的伪逆 # 注意:A 不一定是方阵(主成分数 < 特征数), # 所以不能直接求逆,要用 pinv A_inv = np.linalg.pinv(A) # 反推公共因子值(即未风化时的主成分得分) factors = A_inv @ centered # 再由公共因子重构原始空间 recovered = mu + A.T @ factors return recovered这段代码里最值得注意的是np.linalg.pinv。比赛复盘时很多人会在这里踩坑:负载矩阵 A 的形状是 (主成分数, 特征数),不是方阵,直接调用np.linalg.inv会直接报错,必须用伪逆。伪逆的本质是最小二乘意义下的最优解,当 A 的列数多于行数时,它给出的解是范数最小的那一个。
3.3 主成分个数与累计贡献率的选择
论文中针对铅钡玻璃-C 纹选取了 8 个主成分,保留 97% 以上的成分构成信息。这个选择不是越高越好,而是有讲究的。
主成分个数太少,会丢失重要的成分变化信息;太多,则模型复杂度和过拟合风险上升。数学建模中通常以 85% 作为累计贡献率的下限,但对这道题来说,由于最终要做因子反演,主成分个数直接影响负载矩阵 A 的列数,进而影响反推结果的稳定性。
这里涉及一个数字上的权衡:如果只选前 3 个主成分,累计贡献率可能只有 70%,反推出来的未风化成分会严重失真;如果选到 8 个,贡献率到了 97%,但 A 的维度变大,伪逆计算对噪声的放大效应也会增强。论文里的选择是偏保守的,优先保证还原精度,牺牲了一部分运算效率。
从论证角度看,这个选择是站得住脚的。因子分析在这里不是一个"预测精确数值"的工具,而是一个"还原成分比例"的工具。对于这种目标,尽可能保留信息的优先级应该高于模型简洁性。这也解释了为什么问题三的分类只用区间判断就够了——分类关心的是样本落在哪个范围,不需要精确还原,而问题一需要精确的成分比例,必须用因子分析。
4. 分类、亚类划分与敏感性分析:建模落地与常见问题排查
问题二、三、四本质上是一套"分类 + 亚类划分 + 稳定性检验"的组合拳。这部分看起来没有因子分析那么惊艳,但恰恰是建模完整性的体现,也是数据挖掘类赛题最容易丢分的地方。
4.1 亚类划分:为什么用最短距离法聚类二氧化硅
问题二要求在每一类(高钾/铅钡、风化/未风化)内部再做亚类划分。论文的做法分两步:
第一步,用主成分分析找出贡献率最高的成分。四类数据各自的 PCA 结果中,SiO2 都是差异最大的一个成分。于是选择二氧化硅作为亚类划分的依据。
第二步,用最短距离法(single linkage)对二氧化硅数据聚类。最短距离法的核心定义是:两类之间的距离等于两类中最近的两个样本之间的距离。
D_pq = min { d(x_i, x_j) : x_i ∈ G_p, x_j ∈ G_q }合并后新类到其他类的距离递推公式为:
D_rk = min { D_pk, D_qk }这个方法的优点是计算简单,每次只要找距离矩阵中的最小非对角元素即可合并。缺点在算法层面也很明显:它对噪声和离群点非常敏感,只要有一个异常样本离得远,就可能让两个本不该合并的类被拉在一起。
论文实际是把二氧化硅数据聚类成两个聚点,然后取两个聚点的平均值作为分界,把数据分成两个区间。以铅钡未风化类型为例,两个聚点分别是 64.82 和 52.49,分界点在两者平均值附近,最终得到区间 [50.61, 58.70) 和 [58.70, 69.71]。
这种区间划分方法在数学建模中的合理性在于:它的结果可解释性极强。给定一个新样本,只需要检查它的 SiO2 落在哪个区间,就能直接判定亚类。不需要重新训练模型,不需要调参,一个 Excel 公式就能完成分类。
4.2 敏感性分析的"双面结论"说明了什么
敏感性分析是这道题的一个关键得分点,论文中给出了两个看似矛盾的结论:
问题二对四类已知数据的二氧化硅添加 3% 随机扰动后再分类,发现有三类数据的分类误差超过 20%,于是判定"分类结果的敏感性较强"。但问题三对表单 3 的未知数据做同样的 3% 扰动,分类结果与原来完全一致,判定"敏感性不强"。
同一个方法出现两种结论,这不是模型出错,而是两类数据本身的分布特性不同。训练集数据(表单 2)的二氧化硅值分布比较分散,很多样本靠近区间边界,3% 的扰动足以让它们越过边界跳到另一个区间;而表单 3 的未知样本则大多落在区间内部,距离边界足够远,同样幅度的扰动不会改变归属。
这里有一个值得深挖的细节:区间分类法天然带有"边界敏感"问题。只要样本落在边界附近,微小扰动就会导致分类翻转。解决方法不是改用更复杂的模型,而是检查样本在区间内的位置。如果所有样本都离边界较远,说明分类结果是稳定的;如果存在大量靠近边界的样本,则需要考虑把概率输出而不是硬分类作为结果。
4.3 灰色关联度分析的完整计算流程
问题四要求分析化学成分之间的关联关系,论文选用灰色关联度分析。这个方法特别适合考古数据,因为它不要求数据满足正态分布或线性关系假设,在样本量小、数据不完全的情况下也能给出稳定的相似性度量。
计算流程分为五步:
import numpy as np def grey_relational_analysis(reference, comparison, rho=0.5): """ reference: 参考数列,形状 (n,) comparison: 比较数列,形状 (m, n),每行是一个化学成分序列 rho: 分辨系数,一般取 0.5 """ # 第一步:无量纲化(均值化) ref = reference / reference.mean() comp = comparison / comparison.mean(axis=1, keepdims=True) # 第二步:计算差值序列 delta = np.abs(comp - ref) # 形状 (m, n) # 第三步:求两级最小差和最大差 min_delta = delta.min() max_delta = delta.max() # 第四步:计算关联系数 # 公式: ksi = (min + rho*max) / (delta + rho*max) ksi = (min_delta + rho * max_delta) / (delta + rho * max_delta) # 第五步:关联度 = 关联系数均值 degree = ksi.mean(axis=1) return degree关于分辨系数 ρ,论文特意提到 ρ ≤ 0.5463 时分辨力最好,一般取 0.5。ρ 越小,关联系数之间的差异越明显,分辨力越强;ρ 越大,所有系数向 1 靠拢,区分度降低。取 0.5 是权衡分辨力和稳定性的常见选择。
对四类玻璃分别计算关联度并排序,论文得到的结果是:
| 玻璃类型 | 化学成分关联度排序 |
|---|---|
| 铅钡未风化 | SiO2 > CaO > K2O > Al2O3 > MgO > Na2O |
| 高钾未风化 | SiO2 > CaO > K2O > MgO > Na2O |
| 铅钡风化 | SiO2 > Na2O > K2O > CaO > MgO |
| 高钾风化 | SiO2 > Na2O > K2O |
四类数据的共同点是 SiO2 都是关联度最高的成分,这符合玻璃以二氧化硅为主体结构的基本事实。差异在于次要成分的关联顺序不同:未风化样本中 CaO、K2O 与主体成分的关联更强,而风化样本中 Na2O 的关联排名上升。这说明风化过程确实改变了化学成分之间的耦合关系,也侧面印证了风化系数定义的合理性。
4.4 避坑指南:数学建模中最容易翻车的五个细节
拿这套题复盘,以下五个坑几乎每年都有参赛队踩进去。
坑一:直接对全部样本做聚类,不分风化状态。现象:聚类结果总是把风化样本单独分出来,类型反而混在一起。原因:风化差异的方差远大于类型差异,无监督算法优先捕捉方差最大的方向。解决:先按"类型 + 是否风化"分成四个子集,再分别建模。这套题里所有分析都是四类分开做的,这是最大的结构性决策。
坑二:风化系数定义方向错了。现象:系数大的样本反而显示"未风化"。原因:把分子分母放反了。铅钡玻璃风化后是 PbO、BaO 增加,所以要把它们放分子;高钾玻璃风化后是 SiO2 增加,所以要把 SiO2 放分子。解决:定义系数后先拿已知标签的样本做验证,如果系数和标签明显矛盾,优先查定义而不是查数据。
坑三:因子分析反推时用inv而不是pinv。现象:程序报错说矩阵不是方阵。原因:负载矩阵行数是主成分数,列数是特征数,两者不相等。解决:改用np.linalg.pinv,并在注释里写清楚是"最小二乘意义下的逆"。
坑四:敏感性分析只看"对不对",不看"差多远"。现象:随机扰动后分类结果变了,就断言模型不稳。原因:没有区分样本是落在区间中心还是区间边界。解决:除了报告分类准确率,还要报告每个样本到最近区间边界的距离分布。这是验证区间分类法可靠性的最直接指标。
坑五:灰色关联度分析没有做无量纲化。现象:关联度排序完全由数值大的成分主导,SiO2 永远是第一。原因:SiO2 的含量是其他成分的几十倍,绝对差被它垄断。解决:必须先做初值化或均值化处理,把每个序列拉到同一量级,否则整个分析没有意义。
5. 设计取舍:这套方法为什么能work,又会在哪里失效
把整套建模链路拆开看,最值得思考的设计取舍是"先定义风化系数,再走因子分析"这两步的配合。
风化系数本质上是一个领域知识驱动的特征工程,它把物理风化机理压缩成一个标量。它的作用是制造一个"分水岭",让后续所有模型在高信噪比的子空间里工作。这是整道题最关键的一步,因为风化方向相反的数据混在一起时,任何纯数据驱动的算法都会被误导。
因子分析的设计则是把"风化"建模成线性变换,然后靠矩阵求逆反推。这里面隐含了一个强假设:风化过程在严重风化之前是近似线性的。论文的问题假设部分明确写了这一点。这个假设在考古场景中是否成立?短期内规律是成立的,但长期风化会有非线性拐点,比如表面富硅层形成后,反应速率会明显下降。所以这套方法对中度风化样本的还原效果好,对严重风化的样本会系统性偏差。
再看区间分类法。它的优缺点在敏感性分析中暴露得非常充分:好处是透明,任何人都能看懂分类规则;坏处是边界处没有概率概念,样本落在 [58.69, 58.71] 和 [60.00, 60.01] 被同等对待,后者明明更稳定。如果换成高斯混合模型或核密度估计,可以给出更平滑的分类概率,但解释成本会上升。在数学建模竞赛的评分体系里,可解释性往往比精确度更值钱,这个取舍是正确的。
灰色关联度分析是整套方法里最"轻"的一环,它不需要分布假设,也不存在过拟合问题,计算量极小。它的局限在于:关联度是对"几何形状相似性"的度量,不一定是因果关系的度量。论文把关联度排序差异作为风化影响成分关联的证据,这个推断方向只能说合理,不能说严格。
从整体评价这套方法:它赢在流程的闭环。数据探索发现风化方向,定义系数量化风化程度,PCA 降维去噪,因子分析反演,区间分类加敏感性检验,灰色关联度做关系分析,每一步的输入输出都衔接得严丝合缝。它不依赖任何黑匣子模型,每一步都能被检验、被解释、被复现。正是这种"透明 + 可验证"的特质,让它成为数学建模题的一个典范解。
我自己的一个习惯是:拿到任何一份成分数据或观测数据,第一件事不是跑 PCA,而是先定义好"什么是一个样本的本质身份,什么又是被环境污染的噪声"。从那以后我每次做这类反演任务,都强制先走一遍"定义污染系数 → 验证方向 → 再看模型"这个流程,省下的返工时间远超想象。希望帮到你。
本文还有配套的精品资源,点击获取