简介面向全国大学生数学建模竞赛参赛者的2022年C题完整方案资源聚焦古代玻璃制品的成分分析与鉴别。资源基于赛题附件数据完整呈现了从数据预处理、探索性数据分析到建模求解的全流程先围绕玻璃类型、纹饰、颜色与风化关系展开统计挖掘定义风化系数定量刻画风化程度再通过主成分分析降维结合因子分析构建从风化到未风化成分的转换矩阵预测风化前含量随后按成分规律进行亚类划分并开展敏感性与灰色关联度分析揭示不同玻璃类别的化学成分关联差异。数据处理阶段还涵盖缺失值处理、异常值识别与标准化并选取二氧化硅、氧化钾等关键特征。资源为1个PDF文档3.55MB内含论文正文、解题思路与关键代码实现覆盖问题重述、数据探索、模型求解与结果分析等章节适合参赛参考也可用于学习PCA、因子分析、灰色关联度在考古数据中的应用。已有927人学习/下载对需要系统掌握建模全流程的本科生、研究生及数学建模爱好者具有实用价值。1. 玻璃成分与风化预测一次把“预测风化前成分”讲透拿到这题的第一反应很多人会去查风化机理、矿物学资料但我建议把重心放在一个反直觉的结论上风化不是让所有成分均匀流失而是铅钡玻璃和高钾玻璃呈现出完全相反的成分变化方向——铅钡玻璃风化后二氧化硅显著降低、氧化铅和氧化钡升高高钾玻璃恰恰相反二氧化硅显著升高、氧化钾降低。这意味着想预测风化前的成分不能套同一个修正公式必须先给样品分类、再分别建模。这份数学建模C题的完整文档加代码核心就是用主成分分析降维、因子分析反推、区间划分做亚类、灰色关联度分析做成分关联排序四步解决成分预测与鉴别问题。适合正在备赛、或想用PCA和因子分析处理成分数据的从业者。2. 数据预处理与风化系数的确定先把“风化/未风化”这个标签校准2.1 表单关联与总成分核算先建一张可追溯的工作表原题给了三张表表单1是文物的纹饰、颜色、类型与是否风化标签表单2是每个采样点的化学成分百分比表单3是未知类别文物的成分数据。做分析之前我最先做的一件事是把表单2的采样点编号和表单1的样品信息关联起来而不是直接在两张表里分别看。原文档的做法是在表单2右侧增加P列计算成分合计再从A列抽取编号放入Q列然后把表单1的“玻璃类型、是否风化、颜色、饰纹”关联到R到U列。import pandas as pd df1 pd.read_excel(表单1.xlsx) df2 pd.read_excel(表单2.xlsx) # 抽取采样点编号原表A列形如 A01-1 / A01-2取前3位作为样品编号 df2[样品编号] df2[采样点].str.extract(r([A-Z]\d)) # 关联表单1的类别信息 df2 df2.merge( df1[[样品编号, 玻璃类型, 是否风化, 颜色, 纹饰]], on样品编号, howleft ) # 计算成分合计检查数据完整性 component_cols [SiO2,Na2O,K2O,CaO,MgO,Al2O3, Fe2O3,CuO,PbO,BaO,P2O5,SrO] df2[成分合计] df2[component_cols].sum(axis1)这段代码的逻辑是先把采样点编号统一到样品维度再做一次左连接。原表里一个样品可能有多个采样点比如“30 部位1”和“30 部位2”如果不加样品编号直接合并会把同一个样品的两个部位当成独立样品导致后续统计里同一件文物被重复计数。合并后务必检查成分合计是否接近100实践中这些数据是四舍五入后的百分比合计通常在90到105之间浮动不等于100是正常的但低于85就需要警惕是否存在缺失成分项。2.2 风化系数怎么定分类型定义阈值从数据里读原文档直接把风化系数定义成了两个公式铅钡玻璃用氧化铅加氧化钡的和除以二氧化硅高钾玻璃用的是氧化硅与若干易流失氧化物的比值方向刚好相反。为什么要分开定义因为两类玻璃风化的主导化学过程不一样。铅钡玻璃风化时二氧化硅流失氧化铅、氧化钡相对富集所以比值大于1说明风化程度高高钾玻璃则是氧化钾等易溶成分流失二氧化硅相对富集比值同样越大越风化。def weathering_index(row, glass_type): if glass_type 铅钡玻璃: return (row[PbO] row[BaO]) / row[SiO2] if row[SiO2] 0 else None elif glass_type 高钾玻璃: # 高钾玻璃风化后SiO2占比升高K2O降低 return row[SiO2] / (row[K2O] row[Na2O] row[CaO] row[MgO] row[Al2O3]) if ( row[K2O] row[Na2O] row[CaO] row[MgO] row[Al2O3]) 0 else None return None df2[风化系数] df2.apply( lambda r: weathering_index(r, r[玻璃类型]), axis1 )实际操作中会遇到一个问题高钾玻璃里有些成分本身就是0分母可能很小风化系数会爆炸式增长。所以计算后要对系数做排序检查看是否存在异常大的孤立值。原文档的做法值得借鉴——对“铅钡玻璃-A纹”和“铅钡玻璃-C纹”观察风化系数分布发现0.8到1.3之间存在明显跳跃于是把阈值取在1.0对高钾玻璃发现7.23到12.41之间有跳跃阈值取在10。这个思路的本质是不依赖专家经验定阈值而是从数据分布里找断点。2.3 数据探索里藏着的三个结论完成关联和系数计算后原文档从表单1中挖出了三个直接用得上的结论纹饰B的玻璃全部是高钾、全部风化、全部蓝绿色纹饰A的高钾玻璃全部无风化纹饰C的铅钡玻璃风化与未风化比17比7。注意力放在第二个结论上——同等纹饰条件下铅钡玻璃比高钾玻璃更容易风化。这个结论在后文有个重要作用它解释了两类玻璃在风化方向上的差异不是样本噪声而是系统性的化学行为差异。3. PCA因子分析反推风化前成分原理与完整Python实现3.1 为什么是主成分分析因子分析而不是回归问题一的目标是根据风化点的检测数据预测风化前的成分。直觉方案是用回归以风化后成分为自变量、未风化成分为因变量拟合。但这个方案在本题里走不通——未风化的样本数量太少铅钡未风化只有15个点高钾未风化只有11个点而成分指标有12个直接做回归会严重过拟合系数完全不可信。原文档的思路是先主成分分析再因子分析。主成分分析负责降维和去相关因子分析的定位更巧妙它假设当前观测到的风化成分是某个潜在“原始成分向量”经过线性变换后的结果。写成矩阵形式就是x μ Af ε其中f是未风化时的公共因子向量A是因子载荷矩阵。如果能把A求出来f A⁻¹(x - μ)就把风化后成分反推成了未风化成分。3.2 从相关系数矩阵到因子载荷矩阵因子载荷矩阵A不是直接估计出来的而是通过对成分数据的相关系数矩阵做特征值分解得到。原文档公式(5)写的是A [√λ1·u1, √λ2·u2, ..., √λm·um]其中λ是特征值u是对应特征向量。这是因子分析里的标准做法。需要特别注意一个细节如果只保留前k个主成分A不是方阵直接求逆不成立应该用广义逆伪逆。3.3 Python实现与参数说明import numpy as np import pandas as pd from sklearn.preprocessing import StandardScaler from numpy.linalg import pinv def predict_unweathered(pca_data, n_components8): pca_data: DataFrame风化样品的成分数据 n_components: 保留的主成分个数 返回: 未风化成分预测值 # 1. 标准化因子分析假设变量同量纲可比 scaler StandardScaler() X_scaled scaler.fit_transform(pca_data) # 2. 计算相关系数矩阵并做特征值分解 cov_matrix np.cov(X_scaled, rowvarFalse) eigvals, eigvecs np.linalg.eigh(cov_matrix) # 按特征值降序排列 idx np.argsort(eigvals)[::-1] eigvals eigvals[idx] eigvecs eigvecs[:, idx] # 3. 截取前k个主成分构造载荷矩阵A k n_components A eigvecs[:, :k] * np.sqrt(eigvals[:k]) # 4. 对载荷矩阵求伪逆得到转换矩阵 A_inv pinv(A) # 5. 反推公共因子再映射回未风化成分 F X_scaled A_inv.T X_unweathered_scaled F A.T X_unweathered scaler.inverse_transform(X_unweathered_scaled) return pd.DataFrame(X_unweathered, columnspca_data.columns) # 示例对铅钡风化样本预测未风化成分 weathered_pb df2[(df2[玻璃类型] 铅钡玻璃) (df2[是否风化] 风化)] pred_pb predict_unweathered(weathered_pb[component_cols], n_components8)参数说明n_components的选取参考累计贡献率原文档对“铅钡玻璃-C纹”取8个主成分保留97%以上的成分构成信息这个阈值是合理的一般取到85%以上即可。pinv(A)是关键载荷矩阵是10×8的矩形矩阵时只有伪逆可用。注意标准化用的是z-score逆变换时inverse_transform必须使用同一scaler对象否则预测值量纲会错。这里有个小细节值得解释反推得到的X_unweathered_scaled其实是F投影回原始变量空间的结果它不是真实的未风化成分而是“最接近原始成分空间的估计”。所以预测值要和实际考古知识做交叉验证不要盲目当精确值用。3.4 预测值的合理性判断原文档给出了铅钡玻璃-C纹的预测结果表观察表格可以看到一个规律预测的未风化二氧化硅含量大多在50到70之间氧化铅和氧化钡大多降到个位数这与铅钡玻璃未风化样本的统计区间高度重合。这就是一个很好的合理性检验——预测出的分布不应该脱离训练数据的现实范围。如果某条预测出现负值或者超过100那不是化学上的可能而是因子分析模型在边界样本上的数值振荡常见做法是直接截断到0到100之间。4. 亚类划分与未知样本鉴别区间是你唯一需要维护的产出物4.1 主成分选取与主成分区间问题二要求对高钾玻璃、铅钡玻璃做亚类划分。原文档的处理流程是先按“铅钡无风化、高钾无风化、铅钡风化、高钾风化”四类分别做主成分分析选出贡献率高的化学成分然后对选中成分去掉最大值和最小值后形成一个区间数值落在区间里或靠近区间就归入这个类型。from sklearn.decomposition import PCA def pca_interval(group_data, threshold0.85, exclude_outliersTrue): group_data: 某一类样品的成分数据 threshold: 累计贡献率阈值 exclude_outliers: 是否去掉最大值和最小值 pca PCA() pca.fit(group_data) # 累计贡献率达到阈值的主成分数和对应成分 cumvar np.cumsum(pca.explained_variance_ratio_) n_keep np.argmax(cumvar threshold) 1 # 取载荷最大的前k个化学成分 loadings np.abs(pca.components_[:n_keep]).mean(axis0) top_idx np.argsort(loadings)[::-1][:min(n_keep 2, len(loadings))] top_cols group_data.columns[top_idx].tolist() # 去掉最大值和最小值形成区间 intervals {} for col in top_cols: if exclude_outliers: trimmed group_data[col].sort_values()[1:-1] else: trimmed group_data[col] intervals[col] (trimmed.min(), trimmed.max()) return n_keep, top_cols, intervals n_keep, top_cols, intervals pca_interval(df2[df2[玻璃类型] 铅钡玻璃])这段代码里top_idx的选取逻辑值得注意主成分是成分指标的线性组合载荷绝对值越大说明该成分在主成分中贡献越大用前几个主成分的平均载荷来筛“重要化学成分”比单纯看成分方差更稳。去掉最大值和最小值再去求区间是为了防止离群点把区间边界拉得过宽。原文档里铅钡无风化的二氧化硅区间是[50.61, 69.71]如果不做截尾处理这个区间会被个别异常点撑到更宽后续判别时就会把明显不属于该类的样本误收进来。4.2 为什么选二氧化硅做亚类划分变量原文档在主成分分析后选择了二氧化硅作为亚类划分变量。原因是观察四类数据后发现二氧化硅的数据差异最大区分度最好。这里有一个可以直接借鉴的操作思路不是所有主成分都要保留而是从“差异最大的单一指标”入手划分再验证划分结果是否合理。对铅钡无风化样本二氧化硅区间是[50.61, 69.71]用最短距离法聚类得到两个聚点64.82和52.49取平均值58.70作为分界点于是得到两个亚类区间[50.61, 58.70)和[58.70, 69.71]。4.3 Python实现亚类划分from scipy.cluster.hierarchy import linkage, fcluster import numpy as np def subclass_split(data, col, methodsingle): 对指定列做最短距离法聚类返回聚点、分界点和两类区间 data: DataFrame col: 用于划分的化学成分列名 method: single表示最短距离法 X data[[col]].values Z linkage(X, methodmethod) # 聚成两类 labels fcluster(Z, 2, criterionmaxclust) cluster1 X[labels 1] cluster2 X[labels 2] c1_mean, c2_mean cluster1.mean(), cluster2.mean() # 分界点取两个聚点的平均值 split_point (c1_mean c2_mean) / 2 low, high sorted([c1_mean[0], c2_mean[0]]) return { 聚点1: low, 聚点2: high, 分界点: split_point, 区间1: [data[col].min(), split_point], 区间2: [split_point, data[col].max()] } # 对铅钡无风化样本按二氧化硅划分 pb_unweathered df2[(df2[玻璃类型] 铅钡玻璃) (df2[是否风化] 无风化)] result subclass_split(pb_unweathered, SiO2) print(result)linkage(X, methodsingle)里的methodsingle对应最短距离法它定义的类间距离是两类中最邻近样本之间的距离。这个方法的好处是对离群点不敏感——单个异常点不会像重心法那样拉动整个类的中心位置。分界点取两个聚点平均值是一种简化处理也可以用两个类的最小值最大值做区间但用平均值分界能保证两个区间在数轴上无缝衔接避免出现样本落在中间空隙的情况。4.4 未知样本分类与敏感性分析含3%扰动问题三直接复用问题二的区间。把表单3的数据代入区间看在哪个区间就属于哪类。原文档做了敏感性分析给数据加3%随机扰动后再分类看结果是否变化。结果让很多人意外——训练集的四类数据加了3%扰动后有三类误差超过20%说明训练集内的临界样本对扰动很敏感但表单3的数据加3%扰动后分类结果完全不变。这个对比在第5章会详细展开。rng np.random.default_rng(42) def add_noise(data, noise_level0.03): return data * (1 rng.uniform(-noise_level, noise_level, data.shape)) def classify_by_interval(sample, intervals): 按区间判别所属类别 for category, interval_map in intervals.items(): if all(interval_map[col][0] sample[col] interval_map[col][1] for col in interval_map): return category return 未命中 # 对表单3做扰动前后分类对比 df3_noisy add_noise(df3[component_cols].values) df3_noisy pd.DataFrame(df3_noisy, columnscomponent_cols) df3[扰动后分类] df3_noisy.apply(lambda r: classify_by_interval(r, intervals), axis1)扰动幅度是3%这是原文档固定的参数不要随意加大——扰动加大会让临界样本大量跨区间你无法区分是模型不稳定还是扰动本身过大。生成噪声时务必固定随机种子否则每次跑结果都不一样别人复现不了你的敏感性结论。5. 避坑稳定性检验看的是训练集不是验证集5.1 现象样本标签与风化系数标记不一致原文档在数据探索阶段发现样本20、48、30部位1的标签与风化系数标记不一致。按风化系数为准手工调整了标签。原因表单1的“是否风化”是人工标注的定性判断风化系数是定量指标两者在临界样本上会冲突。人工标注依赖肉眼判断表面状态而风化系数反映的是化学成分偏移程度存在主观与客观的系统性差异。解决定义风化系数后以系数阈值为准统一修正标签。但要注意修正逻辑要透明——修正了哪几个样本、为什么修正都要留在文档里。数学建模评审和项目复盘最忌讳的就是“悄悄改了数据”。5.2 现象因子分析反推的成分预测值出现负值或超100%原因主成分截断保留了k个主成分后载荷矩阵A是矩形矩阵伪逆计算的投影会引入数值误差再加上原始数据本身是四舍五入的百分比误差会被放大。解决对预测值做截断处理小于0的置为0大于100的置为100。更稳妥的做法是只对化学上必然存在的成分做反推对波动极大、含量极低的微量成分不做预测直接参考同类样本的中位数填充。5.3 现象训练集扰动误差超20%表单3扰动结果却全部一致原因训练集的亚类区间边界是紧贴样本值生成的临界样本就在边界附近加3%扰动后很容易跨过边界表单3的样本成分大多落在区间内部远离边界的位置小扰动无法推动它们跨区。解决这个现象不是bug而是敏感性分析的正确解读方式。报告里要分开写两个结论模型对临界样本敏感说明边界附近判别能力弱对远离边界的样本稳定说明主类判别可靠。不要试图掩盖训练集敏感性高的事实这恰好说明亚类边界需要留出缓冲带。5.4 现象表单3样本的成分组合与训练集四类都不完全匹配某些未知样本会同时接近两个区间或者所有区间都命中不了。原因真实文物的成分是连续分布的训练集的四类只是四个粗略的簇边界之外存在过渡样本。解决原文档的区间判别是“硬分类”命中哪个算哪个。实际落地时建议改成“软分类”对每个区间计算距离输出最近区间和距离值距离超过阈值就标记为“待人工复核”。不过原文档没有做这步直接用了硬分类这是一个可以留着说的改进点。5.5 现象不同的人跑主成分分析选出的“重要成分”不一样原因主成分分析是纯方差驱动的方法对标准化方式敏感。用z-score标准化和用min-max标准化得到的载荷分布不同选出的重要特征也可能不同。解决固定数据处理流程从标准化方法到主成分保留数都要写清楚。工程上更常见的做法是直接比较标准化前后的稳定性——如果标准化方式换了之后重要成分排序大变说明数据里存在主导方差的高幅值成分比如铅钡玻璃里氧化铅的方差远大于其他成分这时要结合业务知识判断而不是盲目信统计结果。6. 灰色关联度分析量化“成分之间谁和谁更近”6.1 参考数列与比较数列的构造问题四要求分析不同类别玻璃文物化学成分之间的关联关系。原文档选择灰色关联度分析法而不是皮尔逊相关原因是灰色关联度不要求数据线性、不要求正态分布对样本量的要求也低——高钾风化样本只有7个点12个成分做皮尔逊相关矩阵基本是噪声灰色关联度在这种小样本场景下更稳。做法是每个类别内把样本编号作为参考数列Y各化学成分作为比较数列Xi。先做均值化处理消除量纲差异然后计算每个化学成分与参考数列的关联系数再取平均得到关联度。6.2 分辨系数ρ的取值逻辑灰色关联度公式里的ρ是分辨系数一般取0.5。原文档直接用了这个默认值。ρ越小分辨力越强但同时会放大微小差异对排名的影响ρ越大排名越稳但可能把差异抹平。0.5是工程实践里最常用的值如果关联度排序结果在ρ取0.4和0.6时差异很大说明你选的比较序列本身区分度不够。6.3 Python实现与排序结果def grey_relation_degree(reference, compare_series, rho0.5): reference: 参考数列一维数组 compare_series: 比较数列二维数组每行一个化学成分 rho: 分辨系数默认0.5 返回: 每个比较数列的关联度 # 均值化处理 ref reference / reference.mean() comp compare_series / compare_series.mean(axis1, keepdimsTrue) # 计算差序列 diff np.abs(ref - comp) # 两极最小差和最大差 min_diff np.min(diff) max_diff np.max(diff) # 关联系数 coef (min_diff rho * max_diff) / (diff rho * max_diff) # 关联度 关联系数均值 degree coef.mean(axis1) return degree # 示例铅钡无风化类别的关联度排序 pb_unweathered_data df2[(df2[玻璃类型] 铅钡玻璃) (df2[是否风化] 无风化)] ref_series np.arange(1, len(pb_unweathered_data) 1, dtypefloat) compare_data pb_unweathered_data[[SiO2,CaO,K2O,Al2O3,MgO,Na2O]].T.values deg grey_relation_degree(ref_series, compare_data) comp_names [SiO2,CaO,K2O,Al2O3,MgO,Na2O] for name, d in sorted(zip(comp_names, deg), keylambda x: -x[1]): print(f{name}: {d:.4f})参数说明rho的默认值0.5在大多数场景够用如果你发现排序结果对rho过于敏感可以分别用0.4、0.5、0.6跑一遍把三种rho下的排序都列出来结论一致的条数越多越好。mean(axis1)计算的是每个化学成分在所有样本点上的平均关联系数这是灰色关联度的标准做法。6.4 一个收尾技巧用关联度排序反推分类区间合理性做完问题四后把四类的关联度排序放在一起对比铅钡玻璃排序是二氧化硅氧化钙氧化钾氧化铝氧化镁氧化钠高钾玻璃是二氧化硅氧化钠氧化钾氧化钙氧化镁。可以发现二氧化硅关联性在两个类别里都是最强的。从那以后我每次做这类成分数据题都会强制走一遍这个自检流程先看两类玻璃的风化方向是否相反再看亚类区间是否有重叠最后用灰色关联度排序验证分类区间是否和化学直觉一致。三个检查过一遍模型有没有硬伤基本心里有数。这个资源里的文档把每一步的处理原因都写得比较完整代码也能直接跑通省去了自己从零搭流程的时间。希望帮到你。本文还有配套的精品资源点击获取