简介:这份2022全国大学生数学建模C题资料面向备赛国赛、需要复盘古代玻璃文物成分分析与鉴别案例的高校学生及建模爱好者,核心是用数学建模与数据分析方法解决玻璃风化效应、成分预测和类别鉴别问题。压缩包共1个文件,为3.55MB的PDF文档,集中呈现赛题完整求解方案,可对照学习问题重述、数据探索、预处理与四问建模思路。内容覆盖主成分分析降维、因子分析转换未风化成分、敏感性分析、灰色关联度分析及分类预测等关键环节,并结合高钾玻璃与铅钡玻璃的风化规律、亚类划分和关联度比较给出步骤与结果,便于读者梳理建模框架、理解方法选择并借鉴论文表达。已有928人学习,适合希望系统掌握赛题分析、模型构建与结果解释的参赛者查阅。
1. 风化层下的成分账本:这道题到底在考什么
2022 全国大学生数学建模 C 题给的是一张“玻璃文物身份证”和一张“化学成分账本”:前者记录纹饰、类型、颜色和表面是否风化,后者记录十几个氧化物的质量百分比,附件末尾还附了一批类型未知的样品,等你判断归属。它表面上是一道统计题,实质是一次完整的成分数据处理训练——缺失怎么定义、百分比怎么变换、小样本怎么验证、判别阈值怎么定,四件事一件都躲不掉。做完这套流程,你拿到的不只是四个问题的答案,而是一套能迁移到材料检测、文物断源、工业配方归类的分析方法。适合正在准备数学建模、或者手头就有一堆成分百分比表格要分析的人。
真正让人犯难的,在于“风化”这两个字把数据分成了两种性质不同的样本。风化层的成分和本体不是一回事,铅钡玻璃风化后氧化铅会明显流失,高钾玻璃的钾也会跟着环境走。如果一上来就把风化样品和未风化样品混在一起算均值,后面所有结论都是虚的。所以这道题的第一步是把“数据可信边界”画清楚,再谈模型。
2. 把成分表洗成可建模矩阵:缺失值、合计量与对数比
拿到成分表后不要急着画图,先把表结构和缺失机制搞清楚。这一步做扎实,后面四问的结论才有落脚点。我一般按“读表 → 对齐编号 → 判断缺失类型 → 变换”四步走,每一步都留一个可复查的中间结果。
2.1 先分清“未检测”和“含量为零”
坑就在这一列。赛题附件里部分氧化物在某些样品上是空值,这代表当时没有测这一项,不代表含量真的是零。如果你用fillna(0)一把梭,等于把大量“没测”当成“没有”,算均值时分子不变、分母变大,结果系统性偏低。
import pandas as pd import numpy as np # 三张表:文物基本信息、化学成分、未知类别样品 info = pd.read_excel("附件1.xlsx") chem = pd.read_excel("附件2.xlsx") unknown = pd.read_excel("附件3.xlsx") # 列名去空格,避免合并时对不上 for df in (info, chem, unknown): df.columns = [str(c).strip().replace(" ", "") for c in df.columns] # 采样点编码里带部位信息,抽出前导数字得到文物编号 chem["文物编号"] = chem["文物采样点"].astype(str).str.extract(r"^(\d+)")[0] base = chem.merge(info, on="文物编号", how="left", validate="many_to_one") # 识别氧化物列:以氧化/二氧化/三氧化/五氧化开头 oxide_cols = [c for c in chem.columns if c.startswith(("氧化", "二氧化", "三氧化", "五氧化"))] # 强制转数值,非数字一律变 NaN,保留“未检测”语义 base[oxide_cols] = base[oxide_cols].apply(pd.to_numeric, errors="coerce") print(base[oxide_cols].isna().sum().sort_values(ascending=False).head(8))这段的关键动作有三个:str.extract从采样点里抽文物编号,解决“一张表按样品、一张表按采样点”的粒度错位;validate="many_to_one"让合并时如果出现一对多会直接报错,而不是悄悄产生笛卡尔积;errors="coerce"把“<0.1”这类文本或异常符号统一变成 NaN,防止后面求均值时崩掉。
参数上唯一需要斟酌的是how="left"。用左连接保留化学成分表全部行,是因为基本信息表可能缺某些样品,但成分数据不能丢。如果反过来用右连接,你会莫名其妙少掉几条记录,而且不一定马上发现。
接下来看合计量。每种玻璃检测的氧化物项目不完全一样,把所有列直接加起来通常会偏离 100,偏离多少反映的是检测覆盖度:
# min_count=1 保证整行全空时结果为 NaN,而不是 0 base["检测合计"] = base[oxide_cols].sum(axis=1, min_count=1) print(base["检测合计"].describe()) print((base["检测合计"] < 95).sum(), "条记录检测项不完整")如果相当一部分样品合计明显低于 100,说明这些样品有好几项没测。这时候做跨样品比较,就必须限定在同一组“测了哪些项”的样品之间,或者先把缺失项多的样品单独标出来。
2.2 成分数据的闭合效应与 CLR 变换
成分数据有一个绕不开的性质:所有分量加起来是常数(100% 或某个检测总量),这叫闭合效应。一个分量升高,其他分量被动降低,于是直接对百分比做相关分析会得到一堆伪相关——氧化铅和氧化钡看起来强相关,其实只是“此消彼长”的算术结果。
常见的处理是做中心对数比变换(CLR),把单纯形空间里的数据映射到欧氏空间,再做聚类或回归:
def clr_transform(df, cols, eps=1e-6): """对成分列做 CLR:先加极小值避免 log(0),再减去行几何均值""" x = df[cols].fillna(0) + eps logx = np.log(x) # 每行减去该行所有分量对数的均值,得到零和约束下的坐标 return logx.sub(logx.mean(axis=1), axis=0) clr_df = clr_transform(base, oxide_cols)eps的取值是这里唯一的自由参数。取 1e-6 是因为原始百分比量级在 0.01 到 90 之间,太大会扭曲低含量成分的相对关系,太小则在浮点精度上容易出问题。需要注意的是,CLR 里fillna(0)只是为了算几何均值、让变换可执行;真正做统计检验时,原始表里的 NaN 不能替换,两套数据要分开存放。
另一个常见做法是只用“主要成分子集”做 CLR,比如只保留 SiO₂、K₂O、PbO、BaO、Al₂O₃ 这几个区分度高的分量。分量选得越多,微量元素的噪声占比越大,小样本下反而不稳。
提示:CLR 之后再接标准化要谨慎。CLR 结果本身已经是零和约束,再套一层 StandardScaler 等于把约束信息又打散一次,两种处理二选一即可,别叠着用。
到这里,你手上应该有三份东西:一份保留原始缺失的明细表、一份用于聚类的 CLR 矩阵、一份用于描述统计的按类型分组的汇总表。后面所有分析都从这三份出发,不要临时再从 Excel 里重读一遍,否则口径迟早对不上。
3. 风化与类型的关联检验:卡方、效应量和配对差异
第一问要回答“表面风化跟类型、纹饰、颜色有没有关系”,第二问要回答“风化前后成分怎么变”。两个问题的性质完全不同:前者是分类变量之间的关联,后者是连续变量的组间差异。用错检验方法是这一章最常见的翻车点。
3.1 用卡方和 Cramér's V 判断类型、纹饰、颜色与风化的关系
类型、纹饰、颜色、风化状态都是分类变量,常规做法是列联表加卡方独立性检验。但卡方只给 p 值,p 值小只说明“不太可能是巧合”,不说明关联有多强。样本量几十条的时候,p 值对样本量极其敏感,必须补一个效应量,否则你没法判断这个关联是“强到能用作判别依据”还是“统计显著但实际没用”。
from scipy.stats import chi2_contingency import numpy as np def cramers_v(x, y): tab = pd.crosstab(x, y) chi2, p, dof, expected = chi2_contingency(tab) n = tab.values.sum() # 自由度取行列数较小者减一 v = np.sqrt(chi2 / (n * (min(tab.shape) - 1))) # 小样本偏差修正 v_corrected = max(0, v - (tab.shape[0]-1)*(tab.shape[1]-1)/(2*n-2)) small = (expected < 5).sum() return {"chi2": chi2, "p": p, "V": v, "V_adj": v_corrected, "期望频数<5的格子数": small} for col in ["类型", "纹饰", "颜色"]: print(col, cramers_v(base["表面风化"], base[col]))代码里有两个细节值得说。一是expected < 5的格子数:卡方检验的前提是期望频数不能太小,如果超过两成的格子期望频数低于 5,p 值就不可信了。这时候要么合并类别(比如颜色里的深浅相近项合并),要么改用 Fisher 精确检验的扩展版本。二是V_adj的修正公式,样本量小于一百时未修正的 Cramér's V 会系统性高估关联强度,修正后才好做横向比较。
解读时的经验门槛:V 在 0.1 以下基本可以认为无关,0.3 左右算中等,0.5 以上才谈得上强关联。如果某个纹饰类别只有两三个样品,哪怕 V 很大也不要写进结论,那是单个样本撑起来的假象。
另一个容易忽略的点是:风化状态本身可能是“类型”的下游变量。铅钡玻璃和高钾玻璃的耐候性不同,所以“类型与风化相关”这件事,本质上是材料性质差异的体现,不是风化导致了类型。写结论时要区分相关和因果,别把顺序说反。
3.2 风化前后成分怎么比:非参数检验加效应量
第二问要求分析风化对成分的影响并预测风化前含量。这里第一个要澄清的概念是:同一件文物的风化部位和未风化部位不是配对样本,因为采样点不同、部位不同,甚至可能来自不同器物。所以用配对 t 检验是错的,应该按“类型 × 风化状态”分组做独立样本比较。
成分百分比分布通常右偏、样本量小、还带缺失,正态性靠不住,我一般直接用 Mann-Whitney U 检验,并配上秩双列相关系数作为效应量:
from scipy.stats import mannwhitneyu rows = [] for t in base["类型"].dropna().unique(): sub = base[base["类型"] == t] for c in oxide_cols: a = sub.loc[sub["表面风化"] == "风化", c].dropna() b = sub.loc[sub["表面风化"] == "未风化", c].dropna() # 每组至少 3 个有效值才做检验,否则结果没有意义 if len(a) >= 3 and len(b) >= 3: u, p = mannwhitneyu(a, b, alternative="two-sided") # 秩双列相关:范围 -1 到 1,绝对值越大差异越明显 r = 1 - 2*u/(len(a)*len(b)) rows.append({"类型": t, "成分": c, "n_风化": len(a), "n_未风化": len(b), "均值_风化": a.mean(), "均值_未风化": b.mean(), "差值": a.mean() - b.mean(), "p": p, "效应量r": r}) res = pd.DataFrame(rows).sort_values(["类型", "p"]) print(res.head(15))参数上,alternative="two-sided"是默认选择,因为风化既可能让某个成分富集也可能让它流失,没有先验方向。len(a) >= 3这个门槛是我自己加的:两组各只有一两个值时,U 检验的 p 值只能是 0.1 以上的几个离散值,报出来纯属自欺欺人。这类成分在论文里应该写成“样本不足,不做统计推断”,只给描述性均值。
真正决定结论质量的是效应量那一列。p 值受样本量影响,效应量反映差异的实际幅度。氧化铅在高钾玻璃里可能只有 1% 到 2% 的绝对差,p 值却很显著,这种差异对判别毫无帮助;而在铅钡玻璃里,风化前后氧化铅差十几个百分点,才是真正的规律。
关于“预测风化前含量”,稳妥路径不是去拟合一个从风化到未风化的回归,而是用同类未风化样品的稳健统计量(中位数或截尾均值)作为该类型的参考值,再对风化样品给出区间估计。原因很直接:风化样品数量有限,回归模型的外推能力没有保障,而中位数在小样本、偏态分布下更稳。如果一定要建模,用分位数回归比普通最小二乘更合适,因为它对离群值不敏感。
到这里,第一、二问的骨架就有了:一张关联检验表(带效应量)、一张分组差异表(带效应量和样本量)、一组按类型给出的风化前参考值。表格里的数字必须带样本量,这是小样本分析的底线。
4. 从聚类到判别:两类玻璃的分类边界与未知样品归属
第三问是全题的核心:找出高钾玻璃和铅钡玻璃的分类规律,再把未知样品判进去。很多人一上来就丢个 KMeans 或者随机森林,跑出准确率 0.95 就收工,这在答辩环节基本会被问穿——样本量只有几十条,交叉验证的方差比模型间的差异还大。
4.1 先用 KMeans 看自然分群,再决定是否相信标签
聚类的作用不是给答案,而是检查“已知的类型标签”和“数据自身的结构”是否一致。如果 KMeans 在 k=2 时的分群和类型标签高度吻合,说明成分确实能分开这两类,后面的判别模型有物理基础;如果吻合度很差,那要么是特征选错了,要么是风化样品把结构搅乱了。
from sklearn.preprocessing import StandardScaler from sklearn.cluster import KMeans from sklearn.metrics import silhouette_score, adjusted_rand_score # 只用未风化样品探索结构,避免风化层干扰 clean = base[base["表面风化"] == "未风化"].dropna(subset=["类型"]).copy() X = clr_transform(clean, oxide_cols).values for k in range(2, 5): km = KMeans(n_clusters=k, n_init=20, random_state=42).fit(X) sil = silhouette_score(X, km.labels_) ari = adjusted_rand_score(clean["类型"], km.labels_) if k == 2 else float("nan") print(f"k={k} 轮廓系数={sil:.3f} 与标签的ARI={ari:.3f}")n_init=20是因为 KMeans 对初始质心敏感,小样本下多跑几次取最优能显著降低结果随机性。random_state=42保证结果可复现,这在提交论文时很重要。adjusted_rand_score用来衡量聚类结果和真实标签的一致程度,取值在 0 附近表示跟随机分组差不多,接近 1 表示几乎完全一致。
如果 k=2 的轮廓系数只有 0.2 出头,说明结构不明显,这时候不要硬套聚类结论,转而检查是不是该先做特征筛选——把氧化铅、氧化钡、氧化钾这几个判别力最强的分量单独拿出来看二维散点,往往比十几维空间里的轮廓系数直观得多。
这里有个反直觉的点:风化样品在聚类里经常自成一簇。因为它所有成分都被稀释或富集了,整体偏离本体。如果你的聚类把风化状态当成了主要分群依据,那这个模型对“类型判别”毫无价值。所以探索阶段必须把风化样品排除,或者把风化状态作为分层变量。
4.2 决策树和 SVM 的判别规则,以及未知样品的预测
判别模型要解决两个需求:一是准确,二是可解释。赛题要求“给出分类规律”,所以模型必须能翻译成人话,比如“氧化铅大于某阈值判为铅钡”。决策树天然满足这一点,SVM 作为对照验证非线性边界是否存在。
from sklearn.tree import DecisionTreeClassifier, export_text from sklearn.svm import SVC from sklearn.pipeline import make_pipeline from sklearn.model_selection import LeaveOneOut, cross_val_score X = clr_transform(base, oxide_cols).values y = base["类型"].values tree = make_pipeline( DecisionTreeClassifier(max_depth=3, min_samples_leaf=3, class_weight="balanced", random_state=42) ) svm = make_pipeline(StandardScaler(), SVC(kernel="rbf", C=1.0, gamma="scale", class_weight="balanced")) loo = LeaveOneOut() print("决策树 LOO 准确率:", cross_val_score(tree, X, y, cv=loo, scoring="accuracy").mean()) print("SVM LOO 准确率:", cross_val_score(svm, X, y, cv=loo, scoring="accuracy").mean()) tree.fit(X, y) print(export_text(tree[-1], feature_names=oxide_cols, decimals=2))max_depth=3和min_samples_leaf=3是刻意压低的。样本几十条的情况下,不限制深度的树会把每个样品都记住,训练集准确率 100%,换个数据集就崩。class_weight="balanced"应对两类样品数量不平衡,否则模型会偏向多数类。交叉验证用留一法而不是 5 折,是因为总量太小,5 折时每折验证集只有十几条,方差大到没法比较模型。
export_text打印出来的规则就是第三问的核心答案。你会看到类似“氧化铅 ≤ 某值 且 氧化钾 > 某值 → 高钾”的路径。把这些路径整理成阈值表,比报一个准确率有说服力得多。
对未知样品的预测,要在同一套变换下做,千万别对未知样品单独做 CLR:
Xu = clr_transform(unknown, oxide_cols).values pred = tree.predict(Xu) prob = tree.predict_proba(Xu) unknown_out = unknown[["文物采样点"]].copy() unknown_out["预测类型"] = pred unknown_out["置信度"] = prob.max(axis=1) print(unknown_out.sort_values("置信度").head(10))注意prob.max(axis=1)低的那几条,它们大概率落在两类交界处,是论文里必须单独讨论的对象。把它们简单归到某一类而不加说明,是答辩时最容易被追问的地方。合理做法是给出“倾向铅钡,但置信度 0.6,建议补充检测氧化钡”这样的表述。
5. 五处最容易翻车的细节:填零、泄漏与阈值拍脑袋
这一章是本篇的血泪经验。下面五条几乎每一届都有人踩,而且踩了之后结果表面上还挺好看,属于典型的黑匣子。
5.1 数据层的三个坑
第一条:缺失值填零。现象:某类玻璃的氧化铅均值只有 3%,比文献里的典型值低一大截。 原因:大量“未检测”被填成了 0,均值被严重稀释。 解决:缺失保持 NaN,统计时用dropna()分组计算;只有在必须构造完整矩阵(比如聚类)时,才用极小值或中位数替代,并且替代后的副本单独命名,别覆盖原始表。
第二条:把多次采样当成多个独立样品。现象:某件文物因为采了三个部位,在统计里权重变成别人的三倍,卡方检验结果全歪。 原因:附件按采样点记录,直接value_counts会让采样多的文物占主导。 解决:先按文物编号聚合(成分取均值,类别取唯一值),或者在做关联检验时用文物级数据,在做成分分析时用采样点级数据,两套口径分开维护。
第三条:百分比直接做相关分析。现象:氧化铅和氧化钡的相关系数飙到 0.9 以上,写成“两者共生”。 原因:闭合效应,一个涨另一个必然被动跌,或者反过来。 解决:先做 CLR 或对数比变换再算相关;如果坚持用原始百分比,就改用 Spearman 秩相关并明确说明局限。
5.2 建模层的两个坑
第四条:训练集和验证集不分离。现象:模型准确率 0.98,换一批数据就掉到 0.7。 原因:在全体样本上做了特征筛选或标准化,然后再交叉验证,信息从验证折泄漏进了训练过程。 解决:把特征选择、变换、标准化全部塞进 Pipeline,交叉验证只作用于 Pipeline 整体。或者更稳妥,用嵌套交叉验证。
第五条:判别阈值拍脑袋。现象:氧化铅阈值取 5%,未知样品里有两个被判成高钾,但它们的氧化铅是 4.8%。 原因:阈值没有经过敏感性检验,而且风化样品的氧化铅可能已被淋失,同一类型风化前后的值差异很大。 解决:阈值按风化状态分层设定,或者改用未风化样品建立阈值;同时做阈值扫描,看结论在多大范围内稳定(下一章展开)。
这五条的共同特征是:出错时不会报错,只会让你得到一个“看起来很合理”的错误结论。所以每一处都要留一个中间检查点,把处理前后的样本量、均值、分布画出来对比,别等最后写论文时才发现数字对不上。
6. 敏感性分析:把“判为铅钡”变成一句有底气的话
第四问要求比较两类玻璃的成分关联性并对新样品提出鉴别方法,真正拉开差距的是敏感性分析。评委想知道的不是“你的方法准确率多少”,而是“换个假设,你的结论还站得住吗”。
我一般做三组扫描。第一组是阈值扫描:把氧化铅的判别阈值从 1% 扫到 10%,步长 0.5%,记录每个阈值下未知样品的归属变化。第二组是重采样:对训练集做 200 次 bootstrap,每次重训决策树,统计“氧化铅阈值”这条规则的取值分布。第三组是缺失处理方案对比:分别用丢弃、中位数填补、极小值填补三种策略跑完整流程,看最终归类是否一致。
| 判别阈值(PbO) | 判为铅钡的未知样品数 | 与基准方案一致数 | 结论稳定性 |
|---|---|---|---|
| 2.0% | 9 | 7 | 偏低,易把高钾误判 |
| 5.0% | 6 | 6 | 基准 |
| 8.0% | 4 | 4 | 偏保守,可能漏判 |
| 10.0% | 3 | 3 | 明显漏判 |
表格里的规律是:阈值在 4% 到 6% 之间时结论几乎不变,超过 8% 开始漏判,低于 2% 开始误判。把这段区间写进结论,比单给一个数字强得多——“建议阈值取 5%,在 4% 至 6% 范围内判别结论稳定”。
bootstrap 那组还有一个额外好处:你能看到哪些样品在重采样中反复变类。这些样品就是判别边界上的模糊样本,应该在论文里单独列出,说明“其成分同时具备两类特征,建议结合纹饰与颜色辅助判断”。
至于两类玻璃的成分关联性比较,我的做法是分别对高钾组和铅钡组算 CLR 后的相关矩阵,再做矩阵差异检验,而不是把两类混在一起算一个总相关。混算的结果会被类型这个混杂变量主导,看起来所有成分都相关,实际上只是两类玻璃的均值不同。
最后说一句自己的习惯:这道题我做完第一版之后,会把所有中间表导成 CSV 存到同一个目录,命名带上日期和处理策略,比如chem_未风化_中位数填补.csv。因为做到第四问时你一定会想回头改前面的某个假设,没有这些中间产物,就得从头再跑一遍,那种感觉比返工还难受。建模这件事,后悔药就是中间文件。
希望帮到你。
本文还有配套的精品资源,点击获取