简介一份围绕2022年全国大学生数学建模C题的完整方案包面向数学建模参赛者、数据分析与考古方向研究者。内容以古代玻璃文物成分分析为切入点完整展示从数据预处理、数据探索、特征工程到主成分分析、因子分析、敏感性分析及灰色关联度分析的全流程建模思路。方案基于真实赛题撰写针对四个问题逐一给出解决步骤、模型构建与结果解释帮助读者理解风化效应量化、亚类划分、未知样本鉴别及化学成分关联比较。资源共1个PDF文档压缩包大小3.55MB文件为完整论文全文章节结构清晰公式与图表完整。已有927人学习特别适合备赛参考、建模方法复盘、论文写作模板以及考古数据研究入门。1. 2022 全国大学生数学建模 C 题在古代玻璃风化数据里还原真相这道题是典型的“从被污染的数据里还原真相”的数学建模题。文物里出土的古代玻璃埋了几百年后表面成分早就和烧制时不一样了有人拿着风化严重的铅钡玻璃标本测出来氧化铅高达 47%就误以为它本身是高铅配方其实原始配方里二氧化硅占了 80% 以上。2022 年的这道 C 题核心就是三件事一是预测风化前的成分二是给玻璃分类型和亚类三是看化学成分之间的关联关系。适合正在备赛的学生也适合想练数据分析、降维、聚类、敏感性分析这套流程的开发者。整份资源包含完整论文和可运行的代码照着跑一遍基本能把主成分分析、因子分析、灰色关联度这几个模型串成一条完整链路。2. 数据准备与探索两张表单合并出可用的宽表拿到建模题先别急着上模型先把附件里的两张表搞清楚。2.1 表单结构先搞清楚数据和字段这道题给了三张表单。表单 1 记录每件文物的属性信息文物编号、玻璃类型高钾玻璃或铅钡玻璃、表面是否风化、颜色、纹饰。表单 2 是成分检测数据每一行是一个采样点列是各类氧化物的百分比含量比如二氧化硅、氧化钠、氧化钾、氧化钙、氧化铅、氧化钡等。表单 3 是待鉴别的未知文物成分数据。关键问题是表单 1 和表单 2 不是一一对应的。一个文物可能有多个采样点比如“30 部位 1”“30 部位 2”这种标注方式所以不能直接横向拼接必须通过文物编号做关联。我拿到资源后第一件事就是把表单 1 的属性字段合并到表单 2 的每一行上形成一张宽表后面的分析才方便import pandas as pd df_info pd.read_excel(附件.xlsx, sheet_name表单1) df_comp pd.read_excel(附件.xlsx, sheet_name表单2) # 从编号列中提取文物编号再关联表单1的属性 df_comp[文物编号] df_comp[采样点].str.extract(r(\d)) data df_comp.merge( df_info[[文物编号, 玻璃类型, 是否风化, 颜色, 纹饰]], on文物编号, howleft ) print(data.shape) data.head()merge 的关键是on文物编号howleft表示以成分表为主表保留所有行str.extract(r(\d))是为了把“30 部位 1”这种带中文的编号清洗成纯数字。如果直接拿原始字符串去关联大概率匹配不上。这一步做完每一行成分数据就带上了“玻璃类型、是否风化、颜色、纹饰”四个属性字段后续做分组统计就方便了。2.2 风化系数计算为每种玻璃单独定义预处理里最有价值的一步是定义了一个叫“风化系数”的一维指标。论文的思路是不同玻璃类型风化方向相反所以不能用同一个公式计算风化程度必须分别定义。铅钡玻璃风化后二氧化硅显著降低、氧化铅和氧化钡显著升高所以风化系数定义为氧化铅加氧化钡的和除以二氧化硅高钾玻璃相反风化后二氧化硅显著升高、碱性氧化物降低所以风化系数定义为二氧化硅除以氧化钾、氧化钙、氧化镁、氧化铝、氧化铁、氧化铜六项之和import numpy as np lead_barium_cols [氧化铅(PbO), 氧化钡(BaO)] high_K_cols [ 氧化钾(K2O), 氧化钙(CaO), 氧化镁(MgO), 氧化铝(Al2O3), 氧化铁(Fe2O3), 氧化铜(CuO) ] data[风化系数] np.nan mask_lead data[玻璃类型] 铅钡玻璃 mask_K data[玻璃类型] 高钾玻璃 data.loc[mask_lead, 风化系数] ( data.loc[mask_lead, lead_barium_cols].sum(axis1) / data.loc[mask_lead, 二氧化硅(SiO2)] ) data.loc[mask_K, 风化系数] ( data.loc[mask_K, 二氧化硅(SiO2)] / data.loc[mask_K, high_K_cols].sum(axis1) ) data[风化系数].describe()这里要特别注意分子分母的方向搞反了阈值含义就完全变了。论文里对不同类型分别观察风化系数分布发现铅钡玻璃的判别阈值放在 1.0 附近高钾玻璃的阈值放在 10.0 附近。两者差了一个数量级就是因为公式设计方向不同。所以后面对比风化程度时不能跨类型直接比较系数绝对值。2.3 数据探索结论铅钡玻璃比高钾玻璃更怕风化用 pandas 的交叉表做数据探索能发现几个非常规律的结论。第一个规律纹饰 B 的玻璃文物全部是高钾玻璃、全部风化、颜色全部是蓝绿色。第二个规律同样是纹饰 A高钾玻璃全部无风化铅钡玻璃风化与无风化的比例是 11:5纹饰 C 下高钾玻璃全部无风化铅钡玻璃风化比例为 17:73。这个结论直接回答了论文问题一里的“风化与类型、纹饰、颜色的关系”同等纹饰条件下铅钡玻璃比高钾玻璃更容易风化纹饰 B 对应的埋藏环境风化效应更强连相对耐风化的高钾玻璃都全部风化。还有一个操作层面的坑有些采样点的“是否风化”标注与风化系数计算出的结果是矛盾的比如样本 20、48、30 部位 1。论文的处理方式是以风化系数为准手工修正标签。我建议在代码里直接加一列“标签修正”记录这个过程方便后面复查不要默默改原始数据。3. 问题一核心因子分析反推未风化成分数据探索做完问题一的难点在于已知风化后的成分怎么反推风化前的成分论文用因子分析建模。3.1 建模思路把风化当线性变换论文做了一个关键假设风化是一个连续过程在严重风化之前成分指标的变化基本是线性的。因此可以把某件文物的当前成分x看作未风化基础成分f经过线性变换再加上随机误差的结果公式形式是x μ A f ε。其中A是因子载荷矩阵f是公共因子向量对应未风化时的成分特征。这里用主成分分析来估计载荷矩阵A先求成分数据的相关系数矩阵再对相关系数矩阵做特征分解取前若干个特征值对应的特征向量构造A [sqrt(λ1)μ1, sqrt(λ2)μ2, ...]。因为A通常是长方阵反推f时不能直接求逆要用伪逆。3.2 代码实现主成分估计载荷矩阵并求伪逆def predict_unweathered(X, keep_ratio0.97): # X 是风化样本的成分矩阵每行一个样本每列一种化学成分 mu X.mean(axis0) sigma X.std(axis0) Z (X - mu) / sigma # 标准化消除量纲影响 R np.corrcoef(Z, rowvarFalse) # 相关系数矩阵 eigvals, eigvecs np.linalg.eigh(R) idx np.argsort(eigvals)[::-1] eigvals eigvals[idx] eigvecs eigvecs[:, idx] cum_ratio np.cumsum(eigvals) / eigvals.sum() n_comp int(np.argmax(cum_ratio keep_ratio)) 1 print(f保留 {n_comp} 个主成分累计贡献率 {cum_ratio[n_comp-1]:.2%}) A np.sqrt(eigvals[:n_comp]) * eigvecs[:, :n_comp].T # 载荷矩阵 A_pinv np.linalg.pinv(A) # 伪逆A 不是方阵时不能直接 inv F Z A_pinv.T # 公共因子得分也就是未风化成分的标准化估计 X_recover F A mu # 还原到原始成分尺度 return X_recover, n_comp代码逻辑分四步一是标准化因为氧化物含量的量纲和方差差异很大不标准化会让第一主成分完全被高含量成分主导二是对相关系数矩阵做特征分解相关系数矩阵对角线都是 1天然消除了量纲影响三是按累计贡献率自动选择保留的主成分个数比如 97% 就对应贡献率 97% 以上的主成分集合四是用伪逆从载荷矩阵反推公共因子再映射回原始尺度。偏离均值越远的样本反推后的成分变化越大。你可以通过对比预测值与原始值的方向来判断是否符合风化规律铅钡玻璃预测后二氧化硅应该上升、氧化铅应该下降高钾玻璃正好相反二氧化硅应该下降、氧化钾等碱性氧化物应该上升。如果符号方向反了说明标准化或载荷矩阵构造有误。3.3 结果解读铅钡玻璃预测成分长什么样论文给出了“铅钡玻璃-C 纹”类别下严重风化点反推未风化成分的预测结果。观察表格里的数据最直观的感受是原始风化样本二氧化硅含量低到 36%氧化铅高达 47%预测未风化后二氧化硅恢复到 51%70% 的区间氧化铅明显回落到个位数。这说明因子分析在这里是有效的预测结果基本符合铅钡玻璃未风化的成分特征。这里有一个使用细节需要留意预测值并不是真实的历史成分而是模型在“风化过程近似线性”假设下的估计。它的用途是给后续分类提供相对干净的输入特征而不是直接作为考古结论。做结果分析时最好同时输出预测值和原始值让读者能直观对比风化前后的变化方向。4. 问题二与问题三区间分类、亚类划分和敏感性分析问题二和三用的是同一套分类逻辑先用主成分分析筛选出影响较大的化学成分再对每个类别的数据去掉一个最大值和一个最小值用剩余数据的最大最小值形成分类区间。4.1 主成分筛选与区间构建论文把样本分为四类铅钡风化、铅钡未风化、高钾风化、高钾未风化。分别对每类做主成分分析选出累计贡献率达到 85% 以上的主成分对应的原始化学成分然后构建区间。比如铅钡未风化类型里二氧化硅、氧化钠、氧化钾、氧化钙、氧化镁、氧化铝六类成分累计贡献率 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]。区间构建的代码实现很直接def build_interval(sub_df, comp_cols, drop1): lo, hi {}, {} for col in comp_cols: vals np.sort(sub_df[col].values) if len(vals) 2 * drop: vals vals[drop:-drop] # 去掉最大值和最小值 lo[col], hi[col] vals.min(), vals.max() return lo, hidrop1表示去掉一个最大值和一个最小值目的是排除极端样本对区间的拉伸让分类区间更集中、更代表这个类别的主体分布。如果某个类别的样本量特别少比如高钾风化只有 7 个样本这时候去掉两个极值可能会让区间变得过窄需要根据实际样本量调整drop值或者直接不丢。判断规则是化学成分数值落在区间里面或者靠近区间就属于这个类别。实际代码里我一般会把“靠近”实现成“落在区间内的成分个数最多的类别”而不是要求所有成分都必须落在区间内因为风化数据噪声较多硬性全命中会导致很多样本分不出去。4.2 基于二氧化硅的最短距离法亚类划分主成分分析结果里二氧化硅在四类数据中的区分度都最大所以论文选择二氧化硅作为亚类划分的依据。亚类划分方法是最短距离法聚类也就是单链接聚类两类之间的距离用两类间最邻近样本的距离表示每次合并距离最近的两类直到合成目标个数。用 scikit-learn 可以一行实现from sklearn.cluster import AgglomerativeClustering sio2 df[二氧化硅(SiO2)].values.reshape(-1, 1) model AgglomerativeClustering(n_clusters2, linkagesingle) labels model.fit_predict(sio2) center_1 sio2[labels 0].mean() center_2 sio2[labels 1].mean() split (center_1 center_2) / 2 print(f两个聚点: {center_1:.2f}, {center_2:.2f}分界点: {split:.2f})linkagesingle对应最短距离法n_clusters2表示聚成两个亚类。聚点就是每个簇的均值分界点取两个聚点的平均值比如铅钡未风化类型的两个聚点是 64.82 和 52.49分界点 58.70于是第一亚类的区间是 [50.61, 58.70)第二亚类是 [58.70, 69.71]。用最短距离法而不是 KMeans是因为一维数据上单链接聚类对离群点更直观且结果稳定可复现。分组完成后要按原始数据重新计算每个亚类的二氧化硅区间形成最终的亚类区间表。4.3 表单 3 的未知样本鉴别与扰动测试问题三就是把表单 3 的未知文物成分代入问题二建好的四类区间落在哪个类别区间就判为哪类。敏感性分析的做法是给数据添加 3% 的随机扰动再重新分类看类别是否改变。论文结论是问题二的亚类划分对 3% 扰动比较敏感四类里有三类误差超过 20%但问题三用同样的 3% 扰动测试所有未知文物的分类结果没有变化说明分类结果对微小扰动不敏感。扰动实验的代码要注意“3%”的实现方式rng np.random.default_rng(42) def perturb(values, ratio0.03): noise rng.uniform(-ratio, ratio, sizevalues.shape) return values * (1 noise) perturbed_sio2 perturb(sio2)这里用的是乘法扰动也就是在原始值基础上浮动 ±3%不会产生负值。如果用原始值 随机数的加法扰动当原始值本身很小、接近 0 时扰动后的结果可能直接变成负值这在化学成分数据里是不合理的。敏感性分析的结论直接写在论文里问题二的亚类划分对扰动敏感问题三的整体分类不敏感。这个反差其实不难理解——亚类分界点附近样本本来就密集微小扰动很容易跨过边界但大类之间成分差距大比如铅钡和高钾的二氧化硅含量区间几乎不重叠扰动 3% 很难跨类。5. 问题四的灰色关联度与常见问题排查问题四分析不同类别玻璃化学成分之间的关联关系用的方法是灰色关联度分析。这类方法很适合样本量小、数据完整度一般的考古数据。5.1 灰色关联度计算代码和参数说明灰色关联度的计算分为四步确定参考数列和比较数列、无量纲化、计算关联系数、计算关联度并排序。这里参考数列可以理解为样本编号序列比较数列是各化学成分序列。无量纲化用均值化处理也就是每列除以该列的均值。关联系数公式里的 ρ 是分辨系数通常取 0.5当 ρ ≤ 0.5463 时分辨力最好。def grey_relational_degree(ref, comp, rho0.5): ref ref / ref.mean() comp comp / comp.mean() delta np.abs(ref - comp) delta_min delta.min() delta_max delta.max() coef (delta_min rho * delta_max) / (delta rho * delta_max) return coef.mean() # 示例计算铅钡未风化样本中氧化钙与参考序列的关联度 # ca data_sub[氧化钙(CaO)].values # ref np.arange(1, len(ca) 1).astype(float) # print(grey_relational_degree(ref, ca))代码里的ref / ref.mean()和comp / comp.mean()是均值化处理把数据统一压缩到 1 附近消除量纲差异。rho0.5是分辨系数越小分辨力越强但超过 0.5463 后分辨力下降所以一般取 0.5。最后返回的关联度是一个 0 到 1 之间的值数值越接近 1说明该化学成分与参考序列的整体变化趋势越一致。5.2 关联度排序不同玻璃类型的成分关系差异论文对四个类别分别做灰色关联度分析按关联度从大到小排序。铅钡玻璃的结果是二氧化硅 氧化钙 氧化钾 氧化铝 氧化镁 氧化钠。高钾玻璃的结果是二氧化硅 氧化钠 氧化钾 氧化钙 氧化镁。可以看到两类玻璃的共同点是二氧化硅的关联性始终最强说明二氧化硅在成分体系中的主导地位不受类型影响。差异点在于氧化钠在铅钡玻璃中排名最末在高钾玻璃中排名第二氧化钙在铅钡玻璃中排名第二在高钾玻璃中降到第四。这个结果和风化机理是能对应上的——铅钡玻璃风化后氧化钠流失严重高钾玻璃本身就以氧化钾、氧化钠为主要助熔剂。灰色关联度这里不用过度解读绝对数值重点是排序。排序结果验证了成分之间的主次关系给考古分析的定性结论提供了量化依据。5.3 常见问题排查五个容易翻车的细节第一不标准化直接做因子分析。现象第一主成分完全被二氧化硅主导载荷矩阵几乎只反映高含量成分的波动。原因二氧化硅原始值在 60 到 90 之间氧化铜可能在 0.5 以下量纲差异巨大。解决对相关矩阵做特征分解前必须标准化或者直接用np.corrcoef从标准化后的数据计算相关系数矩阵。第二载荷矩阵求逆报错。现象矩阵维度不是方阵np.linalg.inv(A)直接抛LinAlgError。原因因子个数小于原始成分个数时A是长方形矩阵不存在严格逆矩阵。解决改用np.linalg.pinv(A)求伪逆。这是理解因子分析建模时最容易卡住的地方。第三敏感度分析扰动方式用错。现象加水随机扰动后原本接近 0 的成分变成了负值分类结果完全崩溃。原因加法扰动对绝对量小、接近 0 的变量影响过大。解决统一用值 × (1 uniform(-0.03, 0.03))的乘法扰动保证扰动后的值始终与原始值同号。第四风化系数阈值跨类型混用。现象把铅钡玻璃的阈值 1.0 用到高钾玻璃上导致所有高钾样品都被判定为严重风化。原因两种玻璃风化方向完全相反风化系数定义公式不同阈值自然不同。解决分组计算、分组设阈值先看风化系数分布找突变点再定阈值。第五直接相信原始标签。现象样本 20、48、30 部位 1 的“是否风化”标签与成分变化规律冲突。原因可能是人工记录误差或采样点命名不规范。解决用风化系数排序辅助判断以模型计算结果为准修正标签并在论文里保留修正记录。6. 二次开发把区间法改造成可复用的分类流水线如果只想跑通论文到第 5 章就够了。但如果要拿这套方法去做其他文物数据集或者想把它接到自己的分析流程里我会建议做一次小重构把论文里的“区间规则”和“类别判断”封装成一个独立模块数据进来直接出结果。一个比较实用的改造是给区间匹配加一个“最近邻兜底”。纯区间规则的问题是如果一个样本某几个成分落在区间外直接会被判为不匹配尤其是风化数据本身噪声大边界附近的样本很容易因为 0.1% 的偏差被判错。我的做法是把区间匹配的命中数作为主判据如果两个类别命中数相同再计算该样本与两个类别均值向量的欧氏距离取距离近的类别class GlassClassifier: def __init__(self, rules, centers): self.rules rules # {类别: {成分: (low, high), ...}} self.centers centers # {类别: 各成分均值向量} def predict(self, row): scores {} for cls, rule in self.rules.items(): hit sum( 1 for col, (lo, hi) in rule.items() if lo row[col] hi ) scores[cls] hit max_hit max(scores.values()) best [c for c, s in scores.items() if s max_hit] if len(best) 1: return best[0] # 平票时用最近邻兜底 return min(best, keylambda c: np.linalg.norm( row[self.cols] - self.centers[c] ))这一段改造本身不复杂但很有用。把区间规则抽成rules字典后换数据集只需要重新生成字典模型代码不用动平票时用距离兜底能明显降低边界样本的误判率。我从那以后每次跑这种“区间敏感性分析”的建模流程都会先做一次这层封装再做扰动实验省了后面反复改判断逻辑的时间。希望这一套拆解对你跑通这道 2022 年数模 C 题有帮助。本文还有配套的精品资源点击获取