简介这份2022年全国大学生数学建模C题资料包面向备赛学生与数学建模初学者围绕古代玻璃文物成分分析与鉴别问题提供赛题解析、数据探索、建模流程和结果说明。内容涉及风化效应与类型、纹饰、颜色的关系化学成分统计规律以及高钾玻璃与铅钡玻璃的分类、亚类划分和未知样品鉴别。资源包共1个PDF文件大小约3.55MB便于集中阅读和离线查阅。已有928人学习下载说明其在赛题复盘和方法参考方面具有一定认可度。读者可从中了解主成分分析降维、因子分析转换、敏感性分析、灰色关联度比较、分类预测和模型评估等完整思路并参考风化系数定义、数据预处理与结果解释适合用于赛题复盘、方法迁移和论文写作框架参考。1. 一份被风化改写过成分的玻璃表值得拆成四个问题2022 年数学建模竞赛 C 题的附件是三张表单一批古代玻璃文物的采样点、纹饰、颜色和表面风化情况一批带完整化学成分检测值的样品还有一批只给了成分、没给类型的待判别样品。难点在于古代玻璃制品的成分分析与鉴别同一批文物风化会让二氧化硅、氧化铅这类指标整体挪位拿风化后的数直接做分类等于用一张被涂改过的成绩单去排名。这道题真正要搭的是一条完整链路——先用数据探索把风化的方向和强度量化再用主成分降维压掉噪声用因子分析把成分反推回未风化状态接着做区间判别和亚类划分最后用灰色关联度解释成分之间的耦合关系。适合正在做数学建模课程设计、统计建模毕设或者需要一套完整的数据分析交付文档的人照搬结构。2. 数据预处理与探索三张表单拼成一张能算的表这道题的分水岭不在建模在预处理。三张表单的粒度不一样表单 1 是文物级别的属性表单 2 是采样点级别的成分检测表单 3 则是待判别的盲样。如果直接把 Excel 丢进模型后面所有结论都是浮的。我一般的做法是先把它们拼成一张宽表让每一行代表一个采样点左边是类型、纹饰、颜色、是否风化四个标签列右边是十三种氧化物的含量列后续无论是做相关分析还是做区间判别都只面对这一张表。2.1 表单 1 与表单 2 拼接的口径与坑位表单 2 的主键是“文物采样点”实际值形如49、30 部位1、42 未风化点1而表单 1 的编号是纯数字。这两列直接 merge 的结果必然是大量空值很多队伍跑到这一步会以为附件缺数据其实是主键没洗干净。处理顺序是先从字符串里用正则抽出数字编号再以编号为键左连接表单 1 的四个属性列最后用“拼接后标签列空值数”做一次自检。import pandas as pd import numpy as np # 表单1文物级属性表单2采样点级成分表单3待判别盲样 info pd.read_excel(附件.xlsx, sheet_name表单1) chem pd.read_excel(附件.xlsx, sheet_name表单2) # 采样点字段混着“30 部位1”“42 未风化点1”这类后缀先抽出纯数字编号 chem[采样点编号] chem[文物采样点].astype(str).str.extract(r(\d)).astype(int) info[采样点编号] info[文物采样点].astype(str).str.extract(r(\d)).astype(int) # 左连接保留全部成分行标签缺失的用“未知”占位方便后续单独筛查 merged chem.merge( info[[采样点编号, 玻璃类型, 表面风化, 颜色, 纹饰]], on采样点编号, howleft ) merged[玻璃类型] merged[玻璃类型].fillna(未知) # 自检标签缺失率超过 10% 说明主键还没洗干净回上一步查 missing_rate (merged[玻璃类型] 未知).mean() assert missing_rate 0.10, f标签缺失率 {missing_rate:.1%}先回头核对编号这段代码做三件事编号抽取、左连接、缺失自检。extract(r(\d))正则会保留第一个连续数字段对30 部位1这种写法足够howleft是刻意的选择成分表是分析主体不能因为标签缺失就被 inner join 丢掉最后的assert是给自己留的后悔药缺失率异常时脚本直接中断而不是一路算到最后才发现样本量对不上。还有一列必须加每种样品的成分合计。因为风化样品的氧化物检测值加总常常只有 60% 到 90%剩下的部分是被环境交换掉或者没检出的组分。合计列不参与建模但它是判断数据可信度的第一道关卡合计明显偏离 100% 且又标注“无风化”的样本通常意味着检测批次有问题。2.2 风化系数把“风化到什么程度”压成一个可比较的数风化是个模糊概念表单里的“风化/无风化”又是人工标注。要把这件事变成可计算的量必须自己定义一个风化系数。这里两类玻璃的方向是相反的这是整道题最容易被忽略的细节钾玻璃风化后二氧化硅占比显著升高、氧化钾显著降低铅钡玻璃风化后二氧化硅显著降低氧化铅、氧化钡显著升高。所以两个系数要分开定义铅钡玻璃用(PbO BaO) / SiO2值越大风化越重高钾玻璃用SiO2 / (K2O CaO MgO Al2O3 Fe2O3 CuO)同样是值越大风化越重。两个式子看起来不对称但方向统一都是“风化后变大”。# 铅钡玻璃助熔剂残留PbO、BaO相对骨架SiO2的富集程度 # 高钾玻璃骨架SiO2相对助熔剂组分K2O 等的富集程度 K_HIGH_K [氧化钾(K2O), 氧化钙(CaO), 氧化镁(MgO), 氧化铝(Al2O3), 氧化铁(Fe2O3), 氧化铜(CuO)] def weathering_coefficient(row): si row.get(二氧化硅(SiO2), np.nan) if row[玻璃类型] 铅钡: flux row.get(氧化铅(PbO), 0) row.get(氧化钡(BaO), 0) return flux / si if si and si 0 else np.nan flux sum(row.get(c, 0) for c in K_HIGH_K) return si / flux if flux and flux 0 else np.nan merged[风化系数] merged.apply(weathering_coefficient, axis1)阈值不是拍脑袋定的。做法是按玻璃类型和纹饰分组把风化系数升序排列画出来找那道“跳跃”铅钡玻璃 A 纹的数据在 0.879 和 1.3 之间有明显断裂C 纹在 0.8 和 1.2 之间跳变所以铅钡的阈值放在 1.0 是站得住的高钾玻璃的系数从 7.23 直接跳到 12.41阈值放在 10.0 更稳。这套“先排序找断点、再定阈值”的流程要写进论文答辩时它比阈值本身更值钱。提示用风化系数重判后会和表单里的人工标注出现少量冲突例如个别采样点标注无风化但系数明显偏高。论文里要明确说明以风化系数为准并列出调整清单不能悄悄改掉。2.3 探索阶段就能定下的三条硬结论数据探索不是走过场这道题有三条结论是直接从交叉表里读出来的后面每一步建模都在用它们。第一条纹饰 B 的样品全部是高钾玻璃、全部风化、颜色全部蓝绿同等条件下铅钡比高钾更容易风化。第二条饰纹 B 代表的环境风化效应更强连相对耐风化的高钾玻璃都被全部侵蚀。第三条两条风化路径的成分变化方向完全相反这正是问题一必须分类建模、不能合并建模的根本原因。把这三条写进论文的问题分析节读者立刻能看出你不是在套模板。很多队伍跳过这一步直接上 PCA结果后面解释主成分含义时说不清风化到底改变了什么只能干巴巴地报贡献率。3. 问题一主成分降维加因子分析反推未风化成分问题一要回答两件事风化与类型、纹饰、颜色有没有统计上的关联以及能不能从风化后的检测值反推出风化前的成分。前一半用相关分析和卡方检验后一半用主成分分析配合因子分析。整条链路的关键假设是某个时刻观测到的成分是原始成分经过连续风化演进的结果可以拆成若干公共因子的线性组合。3.1 为什么降维必须放在因子分析之前十三种氧化物之间存在明显的共线性二氧化硅降下去氧化铅往往同步升上来。如果直接对原始矩阵做因子分解特征值会集中在前一两个方向上载荷矩阵几乎被二氧化硅一个变量主导后面求出来的转换系数全是噪声。先做 PCA 把维度压到能保留 97% 以上信息的主成分个数再在低维空间里构造载荷矩阵既减少运算量也避免单个成分的检测误差被放大。from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA FEATURES [二氧化硅(SiO2), 氧化钠(Na2O), 氧化钾(K2O), 氧化钙(CaO), 氧化镁(MgO), 氧化铝(Al2O3), 氧化铁(Fe2O3), 氧化铜(CuO), 氧化铅(PbO), 氧化钡(BaO), 五氧化二磷(P2O5), 氧化锶(SrO)] # 只对“铅钡玻璃-C 纹”这一子集建模纹饰不同风化环境不同不能混在一起 subset merged[(merged[玻璃类型] 铅钡) (merged[纹饰] C)] X StandardScaler().fit_transform(subset[FEATURES].fillna(0)) # 缺失按 0 计表示未检出 pca PCA().fit(X) cum np.cumsum(pca.explained_variance_ratio_) k int(np.searchsorted(cum, 0.97) 1) # 保留 97% 以上信息所需的主成分个数 print(f保留 {k} 个主成分累计贡献率 {cum[k-1]:.2%})这段的参数选择值得说清楚fillna(0)不是随手写的成分表里的空白通常代表该组分低于检出限填 0 比填均值更符合化学含义0.97这个门槛比常见的 0.85 高因为这一步的输出要拿去反推成分压缩太狠会把风化带来的真实差异一起压掉searchsorted是把“累计贡献率达到阈值”翻译成具体个数避免手数。3.2 载荷矩阵求逆从风化成分到原始成分的转换链路因子分析模型写成矩阵形式是x μ A f ε其中x是观测到的成分向量A是载荷矩阵f是公共因子向量——在这道题里f就是未风化时的成分。要求的东西很明确把A反解出来。标准做法是对相关系数矩阵做特征分解取前 k 个特征值开方后乘对应特征向量拼成A再通过f (AᵀA)⁻¹Aᵀ (x - μ)得到反推值。import numpy as np def factor_loadings(X, k): X: n×m 标准化后的成分矩阵返回载荷矩阵 Am×k R np.corrcoef(X, rowvarFalse) # 相关系数矩阵 vals, vecs np.linalg.eigh(R) order np.argsort(vals)[::-1] vals, vecs vals[order], vecs[:, order] A vecs[:, :k] * np.sqrt(np.clip(vals[:k], 0, None)) # A [sqrt(λi)·ui] return A def backcast(X, A): 由风化后成分反推公共因子近似未风化成分 mu X.mean(axis0) W np.linalg.pinv(A.T A) A.T # 最小二乘意义下的反解 return (X - mu) W.T mu A factor_loadings(X, k) recovered backcast(X, A)这里有个必须点破的地方A是m×k的非方阵np.linalg.inv(A)会直接抛LinAlgError。原文里说“对 A 求逆”工程上落地必须是伪逆或最小二乘解这也是我第一版代码翻车的地方。pinv与显式写inv(A.T A) A.T在数值上等价后者可读性更好写进论文的公式推导也更顺。输出结果按采样点整理成表就是那张“未风化成分预测值”表格式上保留两位小数并在表注里说明反推值只对同一纹饰、同一类型的样本成立。3.3 相关性检验与卡方检验怎么落成论文里的表问题一的定性部分要用统计检验支撑“风化与类型、纹饰、颜色有关系”这个判断。做法是先把类别变量量化跑一遍相关系数矩阵再对三个属性分别做卡方检验同时报告效应量。光有 P 值不够样本量小的时候 P 值很容易失真Phi、Cramérs V、列联系数这三个指标能说明关联的强度。分析项卡方值校正卡方P 值PhiCramérs V列联系数纹饰16.39416.3940.0000.5610.5610.490颜色26.01726.0170.0010.7070.7070.577表面风化3.8612.7580.0490.2720.2720.263读表的方式要在正文说透颜色的关联强度最高说明颜色与玻璃体系高度绑定表面风化的 P 值刚过 0.05 线效应量也最小意味着单看“是否风化”这个二分标签解释力有限——这恰好是为什么要自建风化系数、把二分类变成连续量的理由。相关系数矩阵里纹饰与颜色 -0.492、纹饰与类型 -0.448、颜色与类型 0.554三个属性互相纠缠做解释时要避免把相关当因果。4. 问题二到问题四区间分类、亚类细分与灰色关联度后三个问题是一条链问题二在高钾、铅钡两个大类下继续往下切亚类问题三拿切好的规则去判盲样问题四在亚类基础上解释成分之间的耦合关系。三段共用同一套主成分结果所以代码组织上应该是一个features.py产出主成分和筛选后的成分列后面三个脚本各自读取而不是每个问题重新算一遍。4.1 用累计贡献率卡住成分筛选再做去极值区间问题二的思路是对高钾风化、高钾未风化、铅钡风化、铅钡未风化四类分别做 PCA选出累计贡献率能到 85% 以上的那几个化学成分然后对每个成分去掉一个最大值和一个最小值用剩余数据的极值构成判别区间。去极值的目的是防止单个异常采样点把区间拉得过宽导致所有盲样都落进去。def robust_interval(x): 去掉一个最大值和一个最小值后取剩余极值抗单点异常 x np.sort(np.asarray(x, dtypefloat)) if len(x) 4: x x[1:-1] return float(x.min()), float(x.max()) # 只在贡献率靠前的成分上建区间避免维度爆炸 rule {} for group, sub in merged.groupby([玻璃类型, 表面风化]): cols top_components(sub, FEATURES, threshold0.85) # 返回累计贡献率达标的成分列 rule[group] {c: robust_interval(sub[c].dropna()) for c in cols}四类样本最终落下来的判别区间大致是这样数值直接决定盲样归属论文里必须原样给出类别入选成分数累计贡献率二氧化硅区间其余成分区间节选铅钡未风化688.89%50.61 ~ 69.71Na₂O 0.80~5.74K₂O 0.15~0.32高钾未风化588.89%59.81 ~ 69.33Na₂O 2.10~2.86K₂O 7.37~12.53铅钡风化587.36%4.61 ~ 37.36Na₂O 1.38~2.22K₂O 0.14~1.05高钾风化395.41%92.35 ~ 95.02K₂O 0.59~1.01判别规则写起来很简单目标样品的成分值落在对应区间内或者紧邻区间边界就归入该类多个类别同时命中时比较各成分偏离区间中心的归一化距离取最小者。这个“紧邻”要给出量化定义比如偏离不超过区间宽度的 10%不能留一个模糊的“靠近”。4.2 最短距离法做亚类细分亚类划分的切入点是从区间表里挑差异最大的成分。观察下来二氧化硅的类内跨度最大铅钡风化从 4.61 一路铺到 37.36天然适合做二次切分。方法用最短距离法聚类得到两个聚点后取中点作为分界把每类样本切成两个亚类。from scipy.cluster.hierarchy import linkage, fcluster def split_subclass(values): 最短距离法聚成 2 类返回两类的分界值聚点中点 z linkage(np.asarray(values).reshape(-1, 1), methodsingle) labels fcluster(z, t2, criterionmaxclust) centers [np.asarray(values)[labels c].mean() for c in (1, 2)] return float(np.mean(centers)), centers # 四类样本的二氧化硅聚点与切分边界 bounds {} for group, sub in merged.groupby([玻璃类型, 表面风化]): cut, centers split_subclass(sub[二氧化硅(SiO2)].dropna()) bounds[group] (round(cut, 2), [round(c, 2) for c in centers])跑出来的聚点和边界是铅钡无风化 64.82 与 52.49切分点 58.70高钾无风化 71.22 与 62.77切分点 67.00铅钡风化 33.96 与 15.63切分点 24.80高钾风化 94.66 与 92.57切分点 93.62。methodsingle是刻意的最短距离法对链状分布的数据更敏感成分数据正好是这种形态如果用 ward 或 average边界会被类内密度带偏。划分合理性要配一张分布图来证四个类别的样本在两个亚类里的占比不能出现极端倾斜比如 90% 挤在一边那就说明切分没意义。原文用饼图呈现这个比例思路是对的只是饼图信息密度低改成并排堆叠条形图答辩时更好读。4.3 盲样判别与 3% 扰动的敏感性实验问题三把表单 3 的数据往区间表里代得到归属类型。敏感性分析的做法是给原始数据乘上1 ± 3%的随机扰动重复若干次看分类结果翻转的比例。def sensitivity(X, classify, ratio0.03, times200, seed42): X: 待判样品特征矩阵classify: 单行判别函数返回类型翻转率 rng np.random.default_rng(seed) base np.array([classify(row) for row in X]) flips [] for _ in range(times): noisy X * (1 rng.uniform(-ratio, ratio, X.shape)) pred np.array([classify(row) for row in noisy]) flips.append((pred ! base).mean()) return float(np.mean(flips)) flip_rate sensitivity(form3_values, classify_by_interval)参数上ratio0.03对应原文的 3% 扰动times200是为了让翻转率的标准误压到 3% 以内seed固定保证结果可复现——这三项都要写进论文的实验设置。这里有个必须诚实面对的结论差异问题三的盲样在 3% 扰动下类型基本不变可以说鲁棒性不错但问题二的亚类划分在同样扰动下有若干组误差超过 20%敏感性偏强。原因是亚类边界本质上是按数据分位数切的边界附近的样本天然经不起扰动。处理办法不是把结论抹平而是分层表述类别归属层面稳定亚类归属层面对扰动敏感并给出“样本到边界的距离”这个指标来说明哪些个体不可靠。这种分层结论比一句“鲁棒性较强”有信息量得多。4.4 灰色关联度五步法与四类样本的排序结果问题四要回答的是成分之间的耦合关系。灰色关联度适合这种样本量小、分布未知的场合步骤固定确定参考数列和比较数列、无量纲化、算关联系数、取均值得到关联度、排序。def grey_relational(ref, comp, rho0.5): ref: 1×n 参考数列comp: m×n 比较数列每行一个成分均已无量纲化 delta np.abs(ref - comp) dmin, dmax delta.min(), delta.max() xi (dmin rho * dmax) / (delta rho * dmax) # 关联系数 return xi.mean(axis1), xi # 返回关联度与关联系数矩阵 # 均值化处理消除不同氧化物含量量级差异 def mean_norm(M): return M / M.mean(axis1, keepdimsTrue)rho是分辨系数取 0.5 是通行做法理论上rho ≤ 0.5463时分辨力最好这一点写进论文能显示你不是照抄公式。无量纲化用均值化而不是初值化因为成分数据没有天然的时间起点初值化会把第一个采样点的偶然误差放大到整个序列。四类样本跑出来的排序是铅钡未风化二氧化硅 氧化钙 氧化钾 氧化铝 氧化镁 氧化钠高钾未风化二氧化硅 氧化钙 氧化钾 氧化镁 氧化钠高钾风化二氧化硅 氧化钠 氧化钾 氧化钙 氧化镁。共同点是二氧化硅始终排第一说明骨架成分与其他组分的关系最紧密差异点在铅钡体系里氧化铝进了前四高钾体系里氧化钠的位置更靠前这跟两类玻璃使用的助熔剂不同是对得上的。差异性分析就围绕这两点展开比笼统地说“关联度存在差异”具体得多。5. 避坑与排查五处最容易翻车的地方这一节是我带过几轮之后攒下来的基本每一处都在实际交付里出现过。写在这里是想让后来人少走一遍弯路——这五个坑不解决模型跑得通但结论站不住。5.1 数据层的两个坑第一个坑主键没洗直接 merge拼出来一片空值。现象是连接后标签列大面积缺失跑不下去。原因是表单 2 的采样点字段带着“部位1”“未风化点1”这类后缀和表单 1 的纯数字对不上。解决办法是先用正则抽编号再连接并且加一条缺失率断言超过 10% 就中断。我现在每个需要连接的项目都会写这条断言它拦下来的问题比任何模型检查都多。第二个坑两类玻璃共用一套风化系数。现象是算出来的系数全挤在一个窄区间里找不出阈值断点。原因是铅钡玻璃的风化表现为助熔剂富集系数该用(PbOBaO)/SiO2高钾玻璃的风化表现为骨架富集系数该用SiO2/(K2OCaO…)两者方向相反混用等于把信号抵消掉。解决方式是先按类型分组再分别排序画散点找断点阈值 1.0 和 10.0 就是这么定出来的不是查文献查来的。5.2 模型层的两个坑第三个坑对载荷矩阵直接求逆。现象是LinAlgError: Last 2 dimensions of the array must be square或者强行加上.I之后维度对不上。原因是载荷矩阵是m×k的长方阵k 远小于 m。解决办法是用伪逆pinv或者显式写最小二乘解(AᵀA)⁻¹Aᵀ。顺带说一句伪逆在这里不是妥协是唯一正确的解法——原始变量数多于公共因子数方程组本来就超定。第四个坑敏感性结论前后打架。现象是摘要写“模型对扰动不敏感、鲁棒性较强”正文问题二那节又写“有三类误差超过 20%、敏感性较强”。原因是两处测的对象不同问题三测的是大类归属问题二测的是亚类归属边界附近的样本天然容易被 3% 的扰动推过去。解决办法是把结论分层写清楚并且补一个“样本到边界的归一化距离”指标明确指出哪些个体的判别结果仅供参考。评审最反感的不是模型有弱点是同一篇文档里自相矛盾。5.3 论文层的一个坑第五个坑灰色关联度当成相关系数解释。现象是关联度排出来一堆 0.8、0.9看起来所有成分关系都很强排序失去区分度。原因有两个一是忘了做无量纲化量级大的成分天然占优二是灰色关联度衡量的本来就是序列几何形状的相似性不是线性相关强度数值偏高是正常的不能按相关系数的标准去解读。解决办法是老老实实做均值化处理在正文里明确写出“关联度反映的是变化趋势的贴近程度”并且只做排序比较不给绝对值下结论。注意所有阈值、区间、聚点都要在论文中给出确定来源——要么是数据断点要么是贡献率门槛要么是文献方法不能出现“经多次尝试取 1.5 效果较好”这种表述。6. 验收与答辩怎么证明这套模型站得住交付前的自检我一般分三张表走。第一张是流程完整性表逐条核对数据预处理有没有拼接自检、风化系数有没有分组计算、阈值有没有给出断点依据、问题一到问题四有没有各自的输入输出、敏感性实验有没有固定随机种子。缺一项就回去补不用等评审提问。第二张是数字一致性表把论文正文、图表、代码输出里的关键数字对齐铅钡的阈值 1.0、高钾的阈值 10.0、四个类别的硅含量区间、四个切分点 58.70 / 67.00 / 24.80 / 93.62、四类关联度排序。这几组数字在摘要、正文、附录里出现的次数超过十次错一个地方整篇的可信度就打折。我习惯在代码里把所有阈值集中到一个config.py正文数字全部从代码输出复制不手打。第三张是结论边界表写清楚每条结论的适用范围。比如反推未风化成分的模型只在同类型、同纹饰的子集内成立跨纹饰套用没有依据亚类划分的切分点依赖当前样本分布换成新一批文物需要重新标定灰色关联度的排序结论只在四类样本内部可比不能跨类别比绝对值。把这些边界写在结论节反而是加分项。答辩现场的高频追问集中在这几个方向。问“阈值为什么是 1.0 和 10.0”标准回答是排序后观察断点位置而不是引用某篇文献的推荐值同时说明断点两侧的样本间隔分别是多大。问“既然做了 PCA 为什么还要做因子分析”要讲清两者目标不同PCA 是压缩维度、去掉共线性噪声因子分析是建立可反解的生成模型前者是后者的前置步骤。问“3% 扰动的依据”回答是成分检测的相对误差量级配合固定随机种子和 200 次重复来说明结果稳定同时主动交代亚类层面的敏感性。# 交付前跑一遍完整链路任何一步失败都不允许打包 python src/preprocess.py # 输出 merged.csv含拼接自检日志 python src/q1_factor.py # 输出未风化成分预测表与载荷矩阵 python src/q2_interval.py # 输出四类判别区间与亚类切分点 python src/q3_predict.py # 输出盲样归属与扰动翻转率 python src/q4_grey.py # 输出四类关联度排序 python src/check_consistency.py # 比对 config.py 与论文数字是否一致这套检查脚本是我在一次答辩前夜临时写的当时正文改了一版阈值却忘了同步附录表格被问到差点答不上来。从那以后我每次交付前都强制把这条链路跑一遍宁可多花二十分钟也不要在现场解释数字为什么对不上。这套文档加代码我按问题一到问题四分好了目录附件表单放进去就能直接跑需要的话按文末方式取。希望帮到你。本文还有配套的精品资源点击获取