2022数学建模C题玻璃风化全流程解析：从数据预处理到灰色关联度-育师

简介：这份资源是2022年全国大学生数学建模竞赛C题的完整解题资料，面向备战数模竞赛的高校学生及指导教师，聚焦古代玻璃文物的成分分析与鉴别这一典型赛题。压缩包内共1个PDF文件，约3.55MB，内容涵盖赛题文档与配套代码，系统呈现了从数据预处理、数据探索到建模求解的全过程。资料围绕风化效应展开，涉及主成分分析降维、因子分析构建风化到未风化成分的转换矩阵、敏感性分析验证模型鲁棒性，以及灰色关联度分析比较不同类别玻璃的化学成分关联差异，并给出高钾玻璃与铅钡玻璃的亚类划分方法。读者可借此掌握完整的赛题求解思路、建模步骤与结果解释，理解分类预测与模型评估的实际应用。目前已有926人学习下载，适合需要系统复盘该赛题、提升建模与数据分析能力的学习者参考。

1. 2022全国大学生数学建模C题：一份能跑通全流程的玻璃风化建模资源

如果你正在准备数学建模竞赛，尤其是想找一套完整的、能从头跑到尾的真题代码和文档，那这份2022年全国大学生数学建模竞赛C题的资源包值得你花时间拆一遍。题目本身是古代玻璃制品的成分分析与鉴别，背景不复杂：玻璃埋在地下几千年，表面会风化，化学成分比例会变，你要根据风化后的检测数据反推它没风化时长什么样，还要把高钾玻璃和铅钡玻璃分开，再往下做亚类划分和成分关联分析。听起来像考古，实际上是一道标准的“数据预处理+降维+分类+关联度”综合题。资源包里包含完整论文文档和配套代码，覆盖了主成分分析、因子分析、灰色关联度、最短距离法聚类这几个核心方法。适合两类人：一是第一次打国赛、需要看完整解题链路的新手；二是做过几道题但总在“数据怎么洗、阈值怎么定、敏感性怎么证”上卡住的熟手。下面我按实际复现的顺序，把这份资源拆开讲清楚。

2. 数据预处理与风化系数：从表单1到表单2的关联逻辑

2.1 表单1的探索结论直接决定后续建模方向

表单1记录的是玻璃文物的基本信息：纹饰类型、颜色、玻璃类型、表面是否风化。资源里的论文先做了一轮频次统计和卡方检验，得出几个关键结论。纹饰B的样品全部是高钾玻璃、全部风化、颜色全是蓝绿色；纹饰A的高钾玻璃全部无风化，铅钡玻璃风化与无风化比例约11:5；纹饰C的高钾玻璃全部无风化，铅钡玻璃风化与无风化比例约17:7。这几个数字不是随便看看的，它们直接支撑了两个假设：同等饰纹条件下铅钡玻璃比高钾玻璃更容易风化；饰纹B代表的环境风化效应更强。后续所有建模都建立在这两个判断上。

卡方检验的结果也值得注意。纹饰与类型的卡方值16.394，P值0.000；颜色与类型的卡方值26.017，P值0.001；表面风化与类型的卡方值3.861，P值0.049。效应量化方面，颜色的Phi系数0.707，纹饰0.561，表面风化只有0.272。这说明颜色和纹饰与玻璃类型的关联很强，而表面风化与类型的关联相对弱一些。这个结论会影响你后面做分类时选哪些特征。

2.2 表单2的预处理步骤与风化系数定义

表单2是化学成分数据，每个样品有二氧化硅、氧化钠、氧化钾、氧化钙、氧化镁、氧化铝、氧化铁、氧化铜、氧化铅、氧化钡、五氧化二磷、氧化锶等十几个成分的百分比。原始数据不能直接拿来建模，资源里的做法是先在Excel里做几步关联：

' 在表单2右侧新增P列，计算B到O列的成分合计 =SUM(B2:O2) ' 从A列样品编号中提取编号放入Q列 =LEFT(A2, FIND(" ", A2)-1) ' 通过编号将表单1的玻璃类型、是否风化、颜色、饰纹关联到R、S、T、U列 =VLOOKUP(Q2, 表单1!A:E, 2, FALSE) =VLOOKUP(Q2, 表单1!A:E, 3, FALSE) =VLOOKUP(Q2, 表单1!A:E, 4, FALSE) =VLOOKUP(Q2, 表单1!A:E, 5, FALSE)

这几步的逻辑是：表单2只有化学成分，没有类型和风化标签，必须通过样品编号把表单1的信息挂过来，否则后面没法按类型分组建模。P列的合计用来检查数据是否有缺失或异常，正常样品的成分合计应该在95%到105%之间。

风化系数的定义是这份资源里比较有巧思的地方。论文针对铅钡玻璃和高钾玻璃分别定义了不同的公式：

# 铅钡玻璃风化系数 def weathering_coeff_PbBa(sio2, pbo, bao): return (pbo + bao) / sio2 # 高钾玻璃风化系数 def weathering_coeff_K(sio2, k2o, cao, mgo, al2o3, fe2o3, cuo): return sio2 / (k2o + cao + mgo + al2o3 + fe2o3 + cuo)

铅钡玻璃风化后氧化铅和氧化钡含量升高、二氧化硅降低，所以用(PbO+BaO)/SiO2做比值，值越大风化越严重。高钾玻璃反过来，风化后二氧化硅升高、助熔剂成分降低，所以用SiO2除以其他成分之和。这个定义不是拍脑袋来的，是根据参考文献里两种玻璃风化机理相反这一事实推导的。

2.3 风化阈值的确定与标签修正

算出风化系数后，按玻璃类型和饰纹分成三组分别排序观察。铅钡玻璃-A纹的风化系数在0.879和1.3之间有一个明显跳跃，所以阈值定在1.0；铅钡玻璃-C纹在0.8和1.2之间跳跃，阈值也定在1.0；高钾玻璃在7.23和12.41之间跳跃，阈值定在10.0。这个“找跳跃点定阈值”的做法比拍脑袋定0.5或者1要靠谱得多。

注意：原始数据里有三个样品（编号20、48、30部位1）的“是否风化”标签与风化系数判断不一致，资源里以风化系数为准做了手工调整。这一步在实际操作中很容易被忽略，但不调整的话后面分类会出现标签噪声。

3. 问题一建模：主成分降维+因子分析反推未风化成分

3.1 为什么选因子分析而不是直接回归

问题一的核心任务是：根据风化后的成分数据，预测未风化时的成分。这本质上是一个“逆变换”问题。资源里选的是因子分析模型，思路是把成分向量X分解为公共因子和特殊因子的线性组合，然后通过负载矩阵A的逆矩阵，从风化后的观测值反推公共因子，再还原未风化成分。

为什么不用多元线性回归？因为回归需要有成对的“风化前-风化后”训练数据，但题目只给了风化后的检测值，没有对应的未风化真值。因子分析的好处是它不需要配对标签，只需要成分之间的协方差结构就能估计负载矩阵。当然这个假设比较强——它默认风化过程是线性的、连续的，论文里也明确写了这条假设。

3.2 主成分分析确定保留几个因子

在做因子分析之前，先用主成分分析降维。以铅钡玻璃-C纹为例，论文选取了8个主成分，累计贡献率超过97%。这个数字怎么来的？看特征值排序后的累计贡献率曲线，一般取到85%以上就够了，但这里取到97%是为了保证反推精度。

import numpy as np from sklearn.decomposition import PCA from sklearn.preprocessing import StandardScaler # 假设 data 是 n x m 的成分矩阵，每行一个样品 scaler = StandardScaler() data_std = scaler.fit_transform(data) pca = PCA() pca.fit(data_std) # 打印累计贡献率 cumvar = np.cumsum(pca.explained_variance_ratio_) for i, v in enumerate(cumvar): print(f"前{i+1}个主成分累计贡献率: {v:.4f}") # 选取累计贡献率超过97%的主成分个数 n_components = np.searchsorted(cumvar, 0.97) + 1 print(f"保留主成分个数: {n_components}")

参数说明：StandardScaler做标准化是必须的，因为各成分的量纲虽然都是百分比，但数值范围差异大（二氧化硅可以到90%以上，氧化锶可能只有0.1%），不标准化的话主成分会被大量纲变量主导。n_components取searchsorted找到的第一个超过97%的位置，比手动指定更稳妥。

3.3 负载矩阵求逆与未风化成分预测

拿到主成分后，计算负载矩阵A。A的每一列是一个主成分对应的特征向量乘以特征值的平方根。然后对A求逆，得到从风化成分到公共因子的转换系数矩阵。

# 获取负载矩阵 A eigenvalues = pca.explained_variance_[:n_components] eigenvectors = pca.components_[:n_components].T A = eigenvectors * np.sqrt(eigenvalues) # 对 A 求逆（注意 A 不是方阵，用伪逆） A_inv = np.linalg.pinv(A) # 从风化成分反推公共因子 # x_weathered 是标准化后的风化成分向量 factors = A_inv @ x_weathered # 还原未风化成分（逆标准化） x_unweathered_std = A @ factors x_unweathered = scaler.inverse_transform(x_unweathered_std.reshape(1, -1))

逻辑说明：A是m×k矩阵（m个成分，k个主成分），pinv求的是伪逆，得到k×m的转换矩阵。factors是k维公共因子向量，再乘回A得到m维的未风化成分估计。最后用inverse_transform还原到原始量纲。资源里给出了铅钡玻璃-C纹的完整转换系数矩阵和预测结果表，可以直接对照验证。

提示：因子分析反推的精度高度依赖主成分个数和线性假设。如果风化过程非线性严重，预测值会偏差很大。资源里没有做交叉验证，这是可以改进的地方。

4. 问题二与问题三：主成分降维+最短距离法聚类的分类链路

4.1 四类数据的主成分分析与区间构建

问题二要求对高钾玻璃和铅钡玻璃分别做亚类划分。资源里的做法是先把表单2的数据分成四类：铅钡无风化、铅钡风化、高钾无风化、高钾风化。然后对每一类分别做主成分分析，选出贡献率高的几个化学成分，去掉最大值和最小值后，用剩余数据的最大值和最小值构成一个区间。

铅钡无风化类型选了6个成分（SiO2、Na2O、K2O、CaO、MgO、Al2O3），累计贡献率88.89%，区间分别是[50.61, 69.71]、[0.8, 5.74]、[0.15, 0.32]、[0.38, 2.82]、[0.71, 1.49]、[1.44, 13.65]。高钾无风化选了5个成分，累计贡献率也是88.89%。铅钡风化选了5个成分，累计贡献率87.36%。高钾风化只选了3个成分（SiO2、Na2O、K2O），累计贡献率就到95.41%。

def build_interval(data, feature_cols, cumvar_threshold=0.85): """ 对指定特征列做主成分分析，返回筛选后的特征区间 data: DataFrame feature_cols: 参与分析的列名列表 cumvar_threshold: 累计贡献率阈值 """ X = data[feature_cols].values scaler = StandardScaler() X_std = scaler.fit_transform(X) pca = PCA() pca.fit(X_std) cumvar = np.cumsum(pca.explained_variance_ratio_) n = np.searchsorted(cumvar, cumvar_threshold) + 1 # 用前n个主成分的载荷绝对值之和排序，选出重要特征 loadings = np.abs(pca.components_[:n]).sum(axis=0) important_idx = np.argsort(loadings)[::-1] intervals = {} for idx in important_idx: col = feature_cols[idx] vals = data[col].values vals_sorted = np.sort(vals) # 去掉一个最大值和一个最小值 trimmed = vals_sorted[1:-1] intervals[col] = (trimmed.min(), trimmed.max()) return intervals, n

参数说明：cumvar_threshold取0.85是常规做法，资源里实际取到了0.87到0.95之间，说明作者希望保留更多信息。去掉最大值和最小值是为了避免极端值把区间拉得太宽，导致分类边界模糊。返回的intervals字典可以直接用于后续分类判断。

4.2 最短距离法聚类的亚类划分

选出重要特征后，资源里选了数据差异最大的二氧化硅做亚类划分。用最短距离法聚成2类，得到两个聚点。铅钡无风化的聚点是64.82和52.49，高钾无风化是71.22和62.77，铅钡风化是33.96和15.63，高钾风化是94.66和92.57。然后以两个聚点的平均值作为分界，把二氧化硅数据分成两个区间。

from scipy.cluster.hierarchy import linkage, fcluster from scipy.spatial.distance import pdist def subcluster_by_sio2(sio2_values, n_clusters=2): """ 对二氧化硅含量做最短距离法聚类，返回两个聚点和分界值 """ X = sio2_values.reshape(-1, 1) # 最短距离法对应 method='single' Z = linkage(X, method='single') labels = fcluster(Z, n_clusters, criterion='maxclust') centers = [] for lab in range(1, n_clusters+1): centers.append(sio2_values[labels == lab].mean()) centers.sort() boundary = (centers[0] + centers[1]) / 2 return centers, boundary

逻辑说明：linkage的method参数选'single'就是最短距离法，它定义类间距离为两类中最邻近样本的距离。fcluster按最大类别数切分。centers是两个类的均值，boundary取均值的中点作为分界。这个分界值就是亚类划分的阈值。

4.3 敏感性分析的3%扰动验证

问题三要求对分类结果做敏感性分析。资源里的做法是给表单3的数据添加3%的随机扰动，再代入问题二求得的区间做分类，看结果是否改变。

def sensitivity_test(data, intervals, perturbation=0.03, n_trials=100): """ 对数据进行随机扰动，检验分类结果的稳定性 """ original_labels = classify(data, intervals) change_count = 0 for _ in range(n_trials): perturbed = data.copy() for col in data.columns: noise = np.random.uniform(-perturbation, perturbation, size=len(data)) perturbed[col] = data[col] * (1 + noise) new_labels = classify(perturbed, intervals) if not np.array_equal(original_labels, new_labels): change_count += 1 sensitivity = change_count / n_trials return sensitivity

参数说明：perturbation=0.03对应3%的扰动幅度，n_trials=100是重复次数。sensitivity越接近0说明模型越鲁棒。资源里问题三的结论是扰动后分类结果一致，敏感性不强；但问题二中亚类划分的敏感性较强，有三类数据误差超过20%。这个差异值得注意——分类和亚类划分的稳定性不是一回事。

5. 避坑与排查：这份资源里最容易翻车的五个地方

5.1 风化系数阈值在不同饰纹间不能直接套用

现象：直接把铅钡玻璃-A纹的阈值1.0用到C纹上，发现部分样品分类结果与标签矛盾。原因：不同饰纹代表不同埋藏环境，风化程度分布不同，A纹和C纹的风化系数跳跃点虽然都在1.0附近，但具体分布有差异。解决：按饰纹分组后分别观察跳跃点，A纹和C纹可以共用1.0，但高钾玻璃必须单独定阈值10.0，不能混用。

5.2 主成分分析前忘记标准化导致降维失效

现象：不做标准化直接跑PCA，第一主成分贡献率超过99%，但后续分类效果很差。原因：二氧化硅含量在60%到95%之间，氧化锶可能只有0.1%，不标准化的话PCA会被大量纲变量完全主导，降维等于没降。解决：用StandardScaler做Z-score标准化，让每个成分的均值为0、方差为1，再跑PCA。

5.3 因子分析求逆时用错矩阵维度

现象：对负载矩阵A直接调用np.linalg.inv报错，提示矩阵不是方阵。原因：A是m×k矩阵（m个成分，k个主成分），m不等于k，不能求常规逆。解决：用np.linalg.pinv求伪逆，或者先做A.T @ A再求逆。资源里用的是伪逆，结果一致。

5.4 敏感性分析的扰动方式影响结论

现象：对数据加3%绝对扰动和对数据乘(1+3%)相对扰动，得到的敏感性结论不同。原因：不同成分的数值范围差异大，绝对扰动对小数值成分影响更大。解决：统一用相对扰动，即乘以(1+均匀噪声)，噪声范围[-0.03, 0.03]。资源里用的是相对扰动。

5.5 灰色关联度分析的分辨系数取值

现象：分辨系数ρ取不同值时，关联度排序发生变化。原因：ρ越小分辨力越大，但太小会导致关联度对数据波动过于敏感。解决：按论文里的做法取ρ=0.5，这是灰色关联度分析的标准取值。如果要做对比实验，可以在0.3到0.7之间取几个值看排序是否稳定。

6. 问题四的灰色关联度分析与全流程验证技巧

问题四要求分析不同类别玻璃中化学成分之间的关联关系。资源里用的是灰色关联度分析，具体做法是把每个类别的样品数据作为参考数列，各化学成分作为比较数列，计算关联系数和关联度，然后排序。

灰色关联度的核心公式是关联系数：

% 灰色关联度分析核心代码 function [grey_degree] = grey_relation(X, Y, rho) % X: 参考数列 (n x 1) % Y: 比较数列矩阵 (n x m) % rho: 分辨系数，通常取 0.5 [n, m] = size(Y); % 初值化处理 X_norm = X / X(1); Y_norm = Y ./ Y(1, :); % 计算绝对差序列 delta = abs(X_norm - Y_norm); % 计算关联系数 delta_min = min(delta(:)); delta_max = max(delta(:)); xi = (delta_min + rho * delta_max) ./ (delta + rho * delta_max); % 计算关联度 grey_degree = mean(xi, 1); end

逻辑说明：初值化处理是为了消除量纲影响，每个数列除以自己的第一个元素。delta是参考数列和比较数列的绝对差。关联系数xi的公式里，rho取0.5是标准做法。最后对每个比较数列的关联系数求平均，得到关联度。关联度越大，说明该化学成分与参考数列的关联越强。

资源里的结论是：铅钡玻璃的关联度排序为二氧化硅>氧化钙>氧化钾>氧化铝>氧化镁>氧化钠；高钾玻璃为二氧化硅>氧化钠>氧化钾>氧化钙>氧化镁。两类玻璃的关联度排序差异明显，说明它们的成分关联结构不同。但四个类别之间比较时，差异其实不大，都是二氧化硅关联性最强。

这里有一个我自己的教训。第一次做这道题的时候，我直接把所有类别的数据混在一起跑灰色关联度，结果排序出来跟论文对不上。后来才发现，灰色关联度分析必须按类别分别做，因为不同类别的参考数列不同，混在一起算没有意义。从那以后我每次做关联度分析都强制走一遍“先分组、再分别计算、最后对比”的流程，再也没翻过车。

验证灰色关联度结果是否合理，可以做一个简单检查：把关联度排序与主成分分析的载荷排序对比。如果某个成分在主成分分析中载荷很高，在灰色关联度中排序也应该靠前。资源里二氧化硅在两个分析中都是第一，说明结果自洽。如果出现矛盾，优先检查数据标准化和分辨系数取值。