简介：面向全国大学生数学建模竞赛参赛者与指导教师，这份2022年C题资料围绕古代玻璃文物的成分分析与鉴别展开，可用于赛题复盘、建模思路梳理、课程设计和方法迁移。资源包共1个PDF文件，约3.55MB，以文档形式承载赛题解析与实现说明，结构按问题顺序组织，便于检索查阅。内容覆盖数据预处理与探索、主成分分析降维、因子分析预测风化前成分、高钾玻璃与铅钡玻璃分类及亚类划分，并用敏感性分析检验模型鲁棒性，以灰色关联度分析比较化学成分关联差异。对四个问题均给出较完整的建模步骤、结果解释与模型评价，适合希望理解完整解题链条、积累数据处理与建模方法的中高级备赛读者。已有928人学习。

1. 2022 年数学建模 C 题：背景看似复杂，真正决定成败的是数据细节

2022 年数学建模竞赛的 C 题属于“背景复杂、数据规整”的典型。题干讨论古代玻璃文物的风化效应，看起来像考古研究，翻开附件却是三张结构清楚的成分百分比表。真正拉开差距的不是模型有多复杂，而是谁能把风化系数定义正确、把因子分析的负荷矩阵求逆跑通、知道随机扰动该重复多少次。

这份资料包含完整论文文档和四个问题的实现代码，覆盖四条主线：主成分分析降维、因子分析反推未风化成分、最短距离法做亚类划分、灰色关联度分析成分差异。适合两类读者：正在备战数学建模竞赛、需要一份能照着跑通的完整案例；以及手上有一批成分数据、想补充统计方法落地经验的人。下面按数据处理、建模、避坑、复现的顺序逐层展开。

2. 数据预处理与风化系数：三张表合并成一张建模表

2.1 表单 2 的成分合计与编号抽取

附件里表单 1 是文物档案（纹饰、类型、颜色、是否风化），表单 2 是每件采样点的 14 个化学成分百分比，表单 3 是待鉴别的未知样本。三张表靠“文物采样点编号”关联，但表单 2 的采样点写法不统一，有“26 严重风化点”“49 未风化点”“30 部位 1”这类后缀，直接做字符串匹配合并，通常只有一半能对上。

先把表头拉平、列名标准化，再做三件事：计算成分合计、抽取纯数字编号、把表单 1 的信息合并进来。

import pandas as pd import numpy as np RAW = "附件.xlsx" # 表单1 列序：文物采样点、纹饰、类型、颜色、表面风化 info = pd.read_excel(RAW, sheet_name="表单1") info.columns = ["编号", "纹饰", "类型", "颜色", "表面风化"] # 表单2 有两行表头，第二行才是具体成分名，header=1 跳过合并单元格那一行 comp = pd.read_excel(RAW, sheet_name="表单2", header=1) comp = comp.rename(columns={comp.columns[0]: "采样点"}) # 中文成分名统一映射成英文简写，后面写公式时不用来回翻表 ALIAS = { "二氧化硅(SiO2)": "SiO2", "氧化钠(Na2O)": "Na2O", "氧化钾(K2O)": "K2O", "氧化钙(CaO)": "CaO", "氧化镁(MgO)": "MgO", "氧化铝(Al2O3)": "Al2O3", "氧化铁(Fe2O3)": "Fe2O3", "氧化铜(CuO)": "CuO", "氧化铅(PbO)": "PbO", "氧化钡(BaO)": "BaO", "五氧化二磷(P2O5)": "P2O5", "氧化锶(SrO)": "SrO", "氧化锡(SnO2)": "SnO2", "二氧化硫(SO2)": "SO2", } comp = comp.rename(columns={k: v for k, v in ALIAS.items() if k in comp.columns}) chem_cols = [c for c in comp.columns if c not in ("采样点",)] # 1) 成分合计：判断缺失、漏检、录入错误的第一道筛子 comp["合计"] = comp[chem_cols].sum(axis=1) # 2) 从“26 严重风化点”这类字符串里抠出纯数字编号 comp["编号"] = comp["采样点"].astype(str).str.extract(r"^(\d+)").astype(float) # 3) 左连接把类型、纹饰、风化标签挂到成分表上 df = comp.merge(info, on="编号", how="left")

代码里最值得说的是合计这一列。14 个成分是百分比，理论上应该接近 100，但原始数据里有不少行只有 70 多，原因包括某些成分未检出被留空，或者录入时丢了列。检查分布比直接归一化更安全：

print(df["合计"].describe()) bad = df[(df["合计"] < 85) | (df["合计"] > 110)] print(bad[["采样点", "合计", "类型"]])

参数上，85 和 110 这两个边界是根据这份数据定的经验值：低于 85 的行缺失成分太多，风化系数会失真；高于 110 基本是录入串列，要回原表核对。对低于 85 的行，一般做法是不做归一化，把缺失成分按未检出填 0。原因是归一化会把“没有检出”变成“检出一点”，对风化系数的影响方向不可控，尤其在铅钡玻璃里 PbO 占比很高，稍微缩放就会改变判定结果。

2.2 两个风化系数公式：铅钡与高钾玻璃不能共用分母

风化系数的目的是把“这件文物风化到什么程度”压缩成一个可排序的数。难点在于两类玻璃的风化机理方向相反：铅钡玻璃风化后 PbO、BaO 相对富集，SiO2 流失，所以用“铅钡氧化物之和除以二氧化硅”；高钾玻璃反过来，SiO2 富集、助熔剂成分流失，所以用“二氧化硅除以助熔剂之和”。

铅钡玻璃：

$K = \dfrac{w(\mathrm{PbO}) + w(\mathrm{BaO})}{w(\mathrm{SiO_2})}$

高钾玻璃：

$K = \dfrac{w(\mathrm{SiO_2})}{w(\mathrm{K_2O}) + w(\mathrm{CaO}) + w(\mathrm{MgO}) + w(\mathrm{Al_2O_3}) + w(\mathrm{Fe_2O_3})}$

def weathering_coef(row): if row["类型"] == "铅钡": sio2 = row.get("SiO2", np.nan) # 分母为 0 或缺失时返回 NaN，不要让 inf 混进排序 if not np.isfinite(sio2) or sio2 <= 0: return np.nan return (row.get("PbO", 0) + row.get("BaO", 0)) / sio2 if row["类型"] == "高钾": flux = sum(row.get(k, 0) for k in ["K2O", "CaO", "MgO", "Al2O3", "Fe2O3"]) if flux <= 0: return np.nan return row.get("SiO2", 0) / flux return np.nan # 表单3 没有类型标签，跳过 df["风化系数"] = df.apply(weathering_coef, axis=1)

这两个公式可以用原始数据直接验算，验算通过才说明列名映射没串。高钾玻璃某风化样本 SiO2 为 96.77，K2O、CaO、MgO、Al2O3、Fe2O3 分别是 0.92、0.21、0.81、0.26、0.84，分母合计 3.04，96.77 / 3.04 得到 31.83，和论文表格里的数值完全一致。铅钡玻璃取一件未风化点：SiO2 为 45.02，PbO 为 30.61，BaO 为 6.22，(30.61 + 6.22) / 45.02 得 0.818，也对得上。

提示：验算至少取风化、未风化各一行，两类玻璃各一组。只验一行很容易因为列错位而“看起来对”。

分母保护那句不能省。原始表里有几行 SiO2 是空值，直接做除法会得到 inf，排序时 inf 会把阈值判断全部带偏，而且 pandas 在做 describe 时不会报错，属于典型的静默失败。

2.3 阈值定在 1.0 和 10.0：看跳变，不要凭空设定

有了风化系数，下一步是把它转成“风化/未风化”的二值标签。常见的错误是凭直觉设 0.5 或者 1.0，然后发现高钾玻璃全部被标成未风化。正确做法是按类型、纹饰分组后排序，找相邻样本之间的最大间隔。

def find_threshold(values, lo=0.0, hi=None): """在 [lo, hi] 区间内找相邻排序值的最大跳变，返回跳变中点""" v = np.sort(values[np.isfinite(values)]) hi = hi if hi is not None else v[-1] v = v[(v >= lo) & (v <= hi)] if len(v) < 2: return np.nan, 0.0 gaps = np.diff(v) i = int(np.argmax(gaps)) return (v[i] + v[i + 1]) / 2, float(gaps[i]) groups = df.groupby(["类型", "纹饰"]) for (typ, patt), g in groups: thr, gap = find_threshold(g["风化系数"]) print(f"{typ}-{patt}纹 样本数 {len(g):2d} 阈值 {thr:.2f} 跳变 {gap:.2f}")

跑出来三组结果：

阈值定完之后还要处理一次标签冲突。原始表单 1 已经有“表面风化”这一列人工标注，拿它和系数判定结果对比，会发现个别样本对不上，比如第 20、48 号和“30 部位 1”。这类冲突不要一律以人工标注为准，先看合计列：合计明显偏低的那几件，通常是成分缺失导致系数偏大，应保留人工标签；合计正常的，以系数为准并做记录。论文里把调整过程写清楚，比悄悄改数据更可信。

3. PCA 加因子分析：负荷矩阵求逆那一步在算什么

3.1 先 PCA 再因子分析的理由

14 个成分百分比天然存在一个线性约束：它们加起来约等于 100，这意味着相关系数矩阵的秩最多 13。再叠加样本量问题——铅钡玻璃风化组只有十几条记录，不到 14 个特征数——直接对协方差矩阵求逆会非常不稳定，算出来的转换矩阵里会出现 ±10 这样的系数，预测值直接飞出 0 到 100 的合理范围。

常见做法是先做主成分分析，把有效维度压到 6 到 8 维，再在低维空间里做因子分解。论文里给出的累计贡献率可以作为压缩依据：

高钾风化组三个成分就到 95.41%，是因为这一组 SiO2 占比普遍在 92 以上，方差几乎全被它一个变量占满。这种“一个变量吃掉大部分方差”的结构在百分比数据里很常见，降维本质上是在压制它的支配地位，让 PbO、BaO、K2O 这些真正携带风化信息的成分有机会进入模型。

3.2 负荷矩阵 A 与转换矩阵

因子分析模型写成矩阵形式是 $x = \mu + Af + \varepsilon$，其中 $A$ 是负荷矩阵，$f$ 是公共因子向量。$A$ 的构造方式是对相关系数矩阵做特征分解，取前 $k$ 个特征值开方后乘对应特征向量：

$A = [\sqrt{\lambda_1}u_1,\ \sqrt{\lambda_2}u_2,\ \dots,\ \sqrt{\lambda_k}u_k]$

要由观测到的风化成分反推因子，理论上要对 $A$ 求逆。但 $A$ 是 $p \times k$ 的长方阵（$p=14$，$k=5$ 或 $6$），根本没有严格意义的逆矩阵。工程上的处理是用伪逆，或者改用最小二乘解：

from sklearn.preprocessing import StandardScaler def factor_transform(X, n_components=6): """X: (n_samples, n_features) 同一组风化样本的成分矩阵 返回：重构后的成分矩阵、负荷矩阵、特征值 """ scaler = StandardScaler() Z = scaler.fit_transform(X) # 标准化，消除百分比量纲差异 R = np.corrcoef(Z, rowvar=False) # 相关系数矩阵 (p, p) eigval, eigvec = np.linalg.eigh(R) # 对称矩阵用 eigh，比 eig 稳 order = np.argsort(eigval)[::-1] eigval, eigvec = eigval[order], eigvec[:, order] k = min(n_components, int((eigval > 1e-8).sum())) A = eigvec[:, :k] * np.sqrt(np.clip(eigval[:k], 0, None)) # 负荷矩阵 (p, k) # 因子得分：用伪逆而不是 inv，k < p 时 inv 必报奇异 F = Z @ np.linalg.pinv(A).T # (n, k) Z_hat = F @ A.T # 重构回原始维度 (n, p) return scaler.inverse_transform(Z_hat), A, eigval[:k]

参数说明：n_components按累计贡献率 85% 这条线选，铅钡无风化取 6，高钾风化取 3；np.linalg.eigh专门用于对称矩阵，特征向量正交性比通用eig好，能避免复特征值；pinv的截断阈值默认是 1e-15，如果特征值里有接近 0 的项，可以手动传rcond=1e-6把病态方向直接丢掉。

这里必须说清一个边界问题：上面这段代码得到的Z_hat是把样本投影到因子空间再投影回来的重构值，物理含义是“去掉次要成分噪声后的成分”，不是严格意义上“风化前的真实成分”。论文把它解释为未风化预测值，成立的前提是假设“风化过程对成分的影响是线性的，且同组样本共享同一因子结构”。答辩时如果被追问这一点，能讲清这是线性假设下的近似，比硬说成“还原”要稳。

3.3 结果校验：预测值的成分合计与变化方向都要检查

因子分析跑完不能直接抄进论文，至少做两个方向的校验。第一是成分合计：预测出来的一行 14 个成分加起来还在 90 到 110 之间，说明变换没有把数据尺度搞坏。第二是变化方向：铅钡玻璃风化后 SiO2 流失，反推的未风化值里 SiO2 均值应该高于观测值；高钾玻璃相反，SiO2 应该低于观测值。这两个断言写进代码里，跑不过就回头查列序。

X_bly = df[(df["类型"] == "铅钡") & (df["纹饰"] == "C") & (df["风化系数"] > 1.0)][chem_cols].fillna(0).values pred, A, eig = factor_transform(X_bly, n_components=6) j_si = chem_cols.index("SiO2") # 铅钡玻璃风化后 SiO2 流失，预测的未风化值应整体偏高 assert pred[:, j_si].mean() > X_bly[:, j_si].mean(), "方向反了，检查列序或符号" tot = pred.sum(axis=1) assert tot.min() > 80 and tot.max() < 120, f"成分合计越界：{tot.min():.1f}~{tot.max():.1f}" print("SiO2 观测均值", X_bly[:, j_si].mean().round(2), "预测均值", pred[:, j_si].mean().round(2))

实测下来铅钡 C 纹的 SiO2 观测均值在 20 出头，预测均值能到 55 附近，方向对得上；高钾风化组反过来，SiO2 从 93 左右被拉回到 88 上下。方向对了，再去看具体数值是否落在论文表格的合理区间里。如果预测值里出现负数或者超过 100 的单项，多半是n_components取太多，把噪声方向也当成信号放进来了，往回减一到两维通常就能压住。

4. 亚类划分与未知鉴别：去极值建区间，最短距离切两刀

4.1 去极值后的区间构造

第二问要求在高钾、铅钡两大类内部再划亚类。论文用的方法很朴素：先用 PCA 选出贡献大的几个成分，对每个成分去掉一个最大值和一个最小值，用剩余数据的上下界构成区间，数值落进区间或者靠近区间就归入该类。

def build_interval(x, trim=1, tol=0.05): """去掉首尾各 trim 个极值后构建区间，tol 是容差比例""" x = np.sort(np.asarray(x, dtype=float)) x = x[np.isfinite(x)] if len(x) <= 2 * trim: return None core = x[trim:len(x) - trim] lo, hi = core.min(), core.max() pad = (hi - lo) * tol # 容差，避免边界样本被判出去 return round(lo - pad, 2), round(hi + pad, 2)

trim=1对应去掉一个最高一个最低，这是论文里的做法；tol=0.05是我加的容差，因为纯用 max/min 建区间时，新样本只要稍微超出一点点就被判为“不属于”，类别覆盖率会莫名其妙掉到 70% 以下。加上 5% 的区间宽度作为缓冲，覆盖率和误判率能平衡得比较好。

四个类别的主要成分区间跑出来是这样：

表里高钾风化组的 Na2O 区间是 0，这一列所有样本都是 0，去极值之后区间退化成一个点。这种退化区间在判别时等于“要求新样本 Na2O 必须精确等于 0”，实际数据里只要有一点点检出就会被判出去。处理办法是把这类零方差成分直接从判别条件里剔除，只保留有区分度的列。

4.2 最短距离法聚出两个聚点，切点取中点

亚类划分选的是 SiO2 这一列，因为四类数据里它的方差最大、区分度最好。用最短距离法（single linkage）聚成两类，取两类的中心作为聚点，切点就是两个聚点的算术平均。

from scipy.cluster.hierarchy import linkage, fcluster def two_clusters(x): """对一维数据做最短距离法聚类，返回升序的两个聚点与标签""" x = np.asarray(x, dtype=float).reshape(-1, 1) Z = linkage(x, method="single") # 最短距离法 lab = fcluster(Z, t=2, criterion="maxclust") centers = sorted(float(x[lab == k].mean()) for k in np.unique(lab)) return centers, lab for name, sub in df.groupby(["类型", "表面风化"]): c, _ = two_clusters(sub["SiO2"].dropna().values) print(name, "聚点", [round(v, 2) for v in c], "切点", round(sum(c) / 2, 2))

跑出来的结果和论文表格能对上：铅钡无风化聚点 64.82 和 52.49，中点 58.66，四舍五入后正好是论文里的 58.70；高钾无风化 71.22 和 62.77，中点 67.00；铅钡风化 33.96 和 15.63，中点 24.80；高钾风化 94.66 和 92.57，中点 93.62。这说明论文里的切点就是简单平均，没有用加权或者其他更复杂的规则。这一点对复现很关键，很多人卡在切点对不上，其实只是聚类方法和取中心的规则不一致。

4.3 表单 3 判别与 3% 扰动敏感性测试

表单 3 的样本没有类型标签，判别逻辑是逐条比对四个类别的区间表，命中哪个类就判为哪个类。多个区间同时命中时，按“区间归一化距离”取最小的那个，比简单投票稳定。

def classify(sample, intervals): """sample: dict{成分: 数值}; intervals: {类别: {成分: (lo, hi)}}""" best, best_score = None, np.inf for cls, rules in intervals.items(): score, miss = 0.0, 0 for chem, (lo, hi) in rules.items(): v = sample.get(chem, np.nan) if not np.isfinite(v): miss += 1 continue width = max(hi - lo, 1e-6) # 落在区间内记 0，落在区间外按偏离宽度归一化 score += max(0.0, (lo - v) / width, (v - hi) / width) score += miss * 0.5 # 缺失成分适当惩罚 if score < best_score: best, best_score = cls, score return best

敏感性测试的做法是给每个成分乘一个 ±3% 的随机扰动，重跑判别，看类别是否改变。这里有个容易出问题的地方：只跑一次就下结论。单次随机的结果波动很大，写进论文的数字别人复现不出来。

def perturb_stability(sample, intervals, n_iter=200, amp=0.03, seed=42): rng = np.random.default_rng(seed) hits = [] for _ in range(n_iter): noisy = {k: v * (1 + rng.uniform(-amp, amp)) for k, v in sample.items()} hits.append(classify(noisy, intervals)) return pd.Series(hits).value_counts(normalize=True).round(3)

跑 200 次取类别分布的占比，比单次结论更有说服力。论文里的两个敏感性结论看起来矛盾——亚类划分在 3% 扰动下有三类误差超过 20%，而表单 3 的类别鉴别在同样扰动下几乎不变——原因在于两者的判定边界宽度不同：亚类切点落在样本密集区，稍微一动就跨类；大类鉴别的区间宽得多，3% 的扰动根本推不过去。这个差异本身就是一条可以写进论文讨论的结论，比硬说“模型鲁棒”要诚实。

注意：扰动之后各成分的合计会偏离 100，如果要严格保持百分比约束，扰动后需要按合计重新归一。但归一化会互相抵消一部分扰动幅度，实际灵敏度会偏低，报告结果时最好注明是否做了归一。

5. 避坑排查：五个让论文结论翻车的细节

5.1 风化系数阈值跨类型套用

现象：高钾玻璃用 1.0 做阈值，全部被判成“风化”，或者反过来全部判成“未风化”，第二问的四个类别直接少一半。

原因：两类玻璃的系数尺度差一个数量级。铅钡玻璃的系数集中在 0.2 到 3 之间，高钾玻璃集中在 1.7 到 31 之间，用同一个阈值必然错。

解决：先按类型和纹饰分组，再对每组单独找排序后的最大跳变。铅钡定 1.0，高钾定 10.0，两组阈值分开写进论文，并在表格里注明依据是“相邻样本系数间隔”。

5.2 成分合计不查就归一化

现象：PCA 第一主成分的方差贡献率超过 70%，载荷几乎全压在 SiO2 上，后面几个成分的区分度全被压没了。

原因：部分样本的成分合计只有 70 多，归一化把这些行的每个成分都按比例放大，SiO2 的绝对占比被进一步推高，组间差异被稀释。

解决：先看合计分布，低于 85 的行标记出来单独处理，未检出成分按 0 保留，不做归一化。归一化只对合计在 85 到 110 之间的行使用。

5.3 因子分析直接对负荷矩阵求逆

现象：numpy.linalg.LinAlgError: Singular matrix，或者求出来的转换矩阵系数在 ±50 以上，预测值出现负数和超过 100 的项。

原因：负荷矩阵是 14 乘 5 或 14 乘 6 的长方阵，不存在严格逆矩阵；样本量小于特征数时相关系数矩阵本身就是奇异的。

解决：改用np.linalg.pinv求伪逆，或者用最小二乘解np.linalg.lstsq。同时把主成分数压到累计贡献率 85% 对应的维数，多余的方向直接截断。

5.4 区间退化成一个点

现象：某个类别在判别时永远命不中，表单 3 的样本全被判到另一个类里。

原因：某成分在该类所有样本中取值都是 0（比如高钾风化组的 Na2O），去极值后区间变成 [0, 0]，任何非零值都被判为区间外。

解决：建模前统计每个成分的方差，方差为 0 的列直接从判别规则里剔除，不要让它参与打分。判别函数里对“区间宽度接近 0”的情况也要做保护，否则归一化距离会算出 inf。

5.5 灰色关联度的参考数列选错

现象：关联度排序和肉眼从数据里看出的趋势完全相反，SiO2 排到最后。

原因：把“样本序号”当成了参考数列，或者对已经归一化的数据又做了一次均值化，导致量纲处理重复。

解决：明确参考数列是待比较的那一列指标，比较数列是其余成分列；无量纲化只做一次，均值化和初值化二选一，不要叠加。分辨系数取 0.5，这是分辨力较好的常用值，低于 0.5463 时分辨力更好但结果对噪声更敏感。

6. 灰色关联度排序与复现验收：ρ 取 0.5 之后怎么把结论讲圆

第四问要在四个类别内部比较化学成分之间的关联强度。灰色关联度的计算量不大，关键在于参考数列定义和无量纲化方式要固定下来。

def grey_relational(X, ref_idx=0, rho=0.5): """X: (n_samples, n_indicators)，按列均值化 ref_idx: 参考数列所在列索引 返回：各指标与参考数列的关联度 """ X = np.asarray(X, dtype=float) X = X / X.mean(axis=0) # 均值化，只做一次 ref = X[:, ref_idx] diff = np.abs(X - ref[:, None]) # 绝对差序列 dmin, dmax = diff.min(), diff.max() xi = (dmin + rho * dmax) / (diff + rho * dmax) return xi.mean(axis=0) # 关联系数取均值即关联度

参数只有两个：ref_idx决定以哪个成分为基准，rho是分辨系数。跑完四类数据，排序结果是：铅钡未风化为 SiO2 > CaO > K2O > Al2O3 > MgO > Na2O；高钾未风化为 SiO2 > CaO > K2O > MgO > Na2O；铅钡风化为 SiO2 > CaO > K2O；高钾风化为 SiO2 > Na2O > K2O > CaO > MgO。

真正值得写进论文的不是排序本身，而是两类玻璃的差异位置：铅钡体系里 CaO 排第二，高钾体系风化后 Na2O 挤到第二。这个位次变化正好对应风化机理——铅钡玻璃风化时 CaO 随 BaO 一起相对富集，高钾玻璃风化时 Na2O 的流失比 K2O 更快。把排序和机理对上，结论才立得住；只贴一张关联度表格，评阅时很难给高分。

复现验收我一般按这几条走一遍：风化系数用原始数据验算两组、阈值分组独立确定、PCA 保留维数按累计贡献率 85% 卡、因子分析方向断言通过、四类区间逐列检查零方差、扰动实验跑 200 次取分布。任何一条不过，先停下来查数据，不要急着往下写模型。

血泪经验是早年间做类似题，直接拿别人跑出来的区间表往自己数据上套，结果类别覆盖率只有六成，回头查了半天才发现是列名映射串了两位，SiO2 列里装的是 Na2O 的数。从那以后我每次做成分数据建模，第一步都强制把列名和原始表头逐行对一遍，跑完再打印前五行人工扫一眼。希望帮到你。

本文还有配套的精品资源，点击获取