玻璃文物风化建模：PCA因子分析与聚类灰色关联

简介：面向备战全国大学生数学建模竞赛的本科生与指导教师，这份2022年C题资料围绕古代玻璃文物成分分析与鉴别展开，可用于赛题复盘、建模方法学习与论文写作参考。包内共1个PDF文档，约3.55MB，以研究文档形式集中呈现问题重述、数据探索、模型构建与结果分析，便于按问题一至问题四梳理完整解题链条。内容覆盖数据预处理、风化系数定义、主成分降维、因子分析转换矩阵、亚类划分、敏感性分析及灰色关联度比较等关键环节，并讨论高钾玻璃与铅钡玻璃的风化规律和成分关联差异。已有928人学习下载，适合希望掌握分类预测、降维建模与模型评估流程的参赛者对照研读，也可作为考古文物保护场景中数据分析案例的参考材料。

1. 拿到这道玻璃C题，先别急着把13种氧化物全塞进模型

十个队伍里有八个，第一次跑这道题时会犯同一个错：把 13 种氧化物含量一股脑丢进模型，结果预测出的“风化前二氧化硅”比风化后还低，跟玻璃风化让 SiO₂ 流失这个基本事实正好反着来。这份资源就是围绕这种反直觉的翻车场景做的，它是一套完整的数学建模交付物，含四问的建模文档和配套代码：Python 负责数据清洗、主成分分析、因子分析和聚类，SPSS 负责相关性与卡方检验，MATLAB 负责灰色关联度。核心链路是先用风化系数把“风化/未风化”锚成一个可算的标签，再用降维加因子分析反推风化前成分，接着划亚类、判未知样品、比成分关联性。适合正在啃这道 C 题、想要一份能真跑通而不是只有论文排版的人，也适合拿小样本高维数据练手的人。

2. 数据清洗与风化系数：把三张表单拼成一张能算的表

2.1 为什么先造一个一维判据，而不是直接上模型

表单 2 的样本量按类别拆开只有个位数到十几个，特征却有 13 到 14 个氧化物。这种维度比样本还大的格局，直接做回归或者因子分析，协方差矩阵接近奇异，跑出来的负载矩阵符号都能翻。更麻烦的是，题目里“是否风化”这一列在表单 1 里是人工标注，样本 20、48、30 部位 1 这几个点的标签和成分走势对不上，属于血泪经验级别的脏数据。

所以这套资源的做法是：先用一个物理意义明确的比值，把“风化程度”压缩成一个标量，再用它去校验和修正标签。这个标量就是风化系数。对铅钡玻璃取 (PbO + BaO) / SiO₂，对高钾玻璃取 SiO₂ / (K₂O + CaO + MgO + Al₂O₃ + Fe₂O₃ + CuO)。选这个形式的理由很直接：铅钡玻璃风化后铅钡富集、硅流失，比值变大；高钾玻璃反过来，硅富集、助熔剂流失，比值也变大。两类玻璃的风化方向相反，所以系数要分开定义，不能统一。

2.2 三张表单的读取、编号对齐与成分合计

采样点那一列的格式很乱，有“01部位1”“26 严重风化点”“42 未风化点1”这种混排，靠字符串匹配去 merge 一定出问题。稳妥做法是用正则把开头的数字抠出来当文物编号，再用它去关联表单 1 的纹饰、类型、颜色信息。

import numpy as np import pandas as pd OXIDES = ["SiO2", "Na2O", "K2O", "CaO", "MgO", "Al2O3", "Fe2O3", "CuO", "PbO", "BaO", "P2O5", "SrO", "SnO2", "SO2"] def load_sheet(path, sheet): df = pd.read_excel(path, sheet_name=sheet) df.columns = [str(c).strip() for c in df.columns] # 采样点列名在不同版本里可能是“文物采样点”或“采样点”，用包含匹配兜住 point_col = [c for c in df.columns if "采样点" in c][0] # 只取开头连续数字，避开“部位1”“未风化点2”这类尾巴 df["文物编号"] = df[point_col].astype(str).str.extract(r"^(\d+)").astype(int) for c in OXIDES: # “未检出”“微量”这类文本先转 NaN 再补 0，不要直接 float() 硬转 df[c] = pd.to_numeric(df[c], errors="coerce").fillna(0.0) df["成分合计"] = df[OXIDES].sum(axis=1) return df s1 = load_sheet("附件.xlsx", "表单1") s2 = load_sheet("附件.xlsx", "表单2") s3 = load_sheet("附件.xlsx", "表单3") info = s1[["文物编号", "纹饰", "类型", "颜色", "表面风化"]] df2 = s2.merge(info, on="文物编号", how="left") print("未匹配上的样本数：", df2["类型"].isna().sum(), "/", len(df2))

这段代码有三个参数值得盯住。errors="coerce"决定异常文本是被当成缺失还是直接抛错，做文物保护类数据一定要选前者，原始表里“未检出”出现频率不低。正则用了^(\d+)而不是(\d+)，因为有的采样点写法是“部位1-2”，不加锚点会把部位号也抠出来。成分合计这一列不是废列，后面质检全靠它，正常样本应当落在 100 附近，偏离超过 3 个百分点的行要单独看一眼。

2.3 两类玻璃的风化系数与阈值标定

系数算出来只是第一步，真正决定成败的是阈值定在哪里。这套资源的做法是先把同类样本按系数排序，看哪里出现明显跳跃，再把阈值放在跳跃区间的整数位置上。

PB_OX = ["PbO", "BaO"] K_OX = ["K2O", "CaO", "MgO", "Al2O3", "Fe2O3", "CuO"] def weather_coef(row): if row["类型"] == "铅钡" and row["SiO2"] > 0: return row[PB_OX].sum() / row["SiO2"] if row["类型"] == "高钾": denom = row[K_OX].sum() return row["SiO2"] / denom if denom > 0 else np.nan return np.nan df2["风化系数"] = df2.apply(weather_coef, axis=1) def largest_gap(x): """返回排序后相邻差最大的一对值和差值，用来定位跳跃区间""" x = np.sort(pd.Series(x).dropna().values) if len(x) < 3: return None d = np.diff(x) i = int(np.argmax(d)) return round(x[i], 3), round(x[i + 1], 3), round(d[i], 3) for name, sub in [("铅钡-纹饰A", df2.query("类型=='铅钡' and 纹饰=='A'")), ("铅钡-纹饰C", df2.query("类型=='铅钡' and 纹饰=='C'")), ("高钾", df2.query("类型=='高钾'"))]: print(name, largest_gap(sub["风化系数"]))

在原始附件上，这三组的最大间隙分别落在 0.879 与 1.307 之间、0.8 与 1.2 之间、7.23 与 12.41 之间。铅钡类取 1.0，高钾类取 10.0，都正好落在间隙里。这里有个容易忽略的点：阈值不要取间隙中点。中点是 1.09 和 9.82，看着更“居中”，但它会随样本增删漂移；取 1.0 和 10.0 的好处是量级清晰、写进论文好解释，而且对新增样本更稳。

提示：阈值定完之后，一定要和表单 1 的人工标签做交叉比对。原始数据里有三个样本的标签和系数判定相反，按系数修正，并在文档里写明修正理由，这一步是评审最看重的部分。

2.4 把“是否风化”落成新列，供后续分组建模

THRESH = {"铅钡": 1.0, "高钾": 10.0} df2["判风化"] = np.where( df2["风化系数"] > df2["类型"].map(THRESH), "风化", "未风化") print(pd.crosstab(df2["类型"], df2["判风化"])) # 落盘，后面 PCA、聚类、灰色关联都从这张宽表取数 df2.to_csv("glass_clean.csv", index=False, encoding="utf-8-sig")

到这一步，一张宽表里同时有了成分、类型、纹饰、颜色、人工标签、系数和机器判定标签，后面三个问题的建模全部基于它展开。顺带说一句，encoding="utf-8-sig"是为了让表头在 Excel 里不乱码，这种小细节不影响模型，但会影响队友打开文件时的心情。

3. PCA 打底 + 因子分析反推：把风化后的成分倒推回风化前

3.1 先降维再反推，绕不开的三个理由

第一问后半段要求“根据风化点检测数据，预测其风化前的化学成分含量”。听着像回归，实际上没有配对样本，没人知道同一个样品风化前长什么样。这套资源用的是因子分析的结构：把某时刻的成分看成若干公共因子的线性组合，公共因子取“未风化时的成分向量”，那么只要拿到负载矩阵，就能从观测值反解出因子得分。

但直接对 14 维成分建因子模型会出问题。一是样本量太小，相关系数矩阵秩不满；二是成分之间存在天然约束，所有氧化物加起来接近 100%，属于闭合数据，任意两个成分都带负相关倾向，这种伪相关会把负载矩阵带偏。所以先做 PCA 把维度压下来，用前若干个主成分承载主要变异，再在主子空间里做因子分析。分类场景保留 85% 累计贡献率够用，反推场景要保留到 97%，因为反推对细节更敏感。

3.2 负载矩阵不能直接求逆，这是最常见的错误

教科书上的写法是“对负载矩阵 A 求逆得到转换矩阵”，但 A 是 p×k 的长方阵，当 k < p 时根本没有逆。正确做法是用伪逆，本质是最小二乘意义下的解。

from numpy.linalg import pinv def fit_fa(Z, var_ratio=0.97): """Z: 标准化后的成分矩阵 (n, p)，返回负载矩阵、转换矩阵和主成分个数""" R = np.corrcoef(Z, rowvar=False) lam, V = np.linalg.eigh(R) # eigh 返回升序，必须翻转 order = np.argsort(lam)[::-1] lam, V = lam[order], V[:, order] ratio = np.cumsum(lam) / lam.sum() k = int(np.searchsorted(ratio, var_ratio) + 1) A = V[:, :k] * np.sqrt(lam[:k]) # p x k，列是因子载荷 M = pinv(A) # k x p，等价于 (A^T A)^-1 A^T return A, M, k, lam, ratio

np.linalg.eigh返回的是升序特征值，忘了翻转会让第一主成分变成方差最小的方向，跑出来的结果毫无意义，这是黑匣子型错误里最难查的一类。var_ratio参数取 0.85 还是 0.97 决定 k，第一问反推建议 0.97，第二问选成分建议 0.85，理由前面说过。

3.3 用未风化样本拟合，再预测风化样本

标准化参数必须来自未风化样本，不能用全量数据。全量标准化等于把待预测样本的分布信息提前泄漏进模型，后面做敏感性分析时误差会虚假地小。

X_un = df2.query("类型=='铅钡' and 判风化=='未风化'")[OXIDES].values mu, sd = X_un.mean(axis=0), X_un.std(axis=0, ddof=1) sd[sd == 0] = 1.0 # 防零方差列除零 Z_un = (X_un - mu) / sd A, M, k, lam, ratio = fit_fa(Z_un, var_ratio=0.97) print("保留主成分个数：", k, "累计贡献率：", round(ratio[k-1], 4)) X_w = df2.query("类型=='铅钡' and 判风化=='风化'")[OXIDES].values Z_w = (X_w - mu) / sd F_w = Z_w @ M.T # 公共因子得分，即未风化成分的标准化取值 Z_hat = F_w @ A.T # 回到标准化成分空间 X_hat = Z_hat * sd + mu # 还原量纲 print("预测成分合计：", np.round(X_hat.sum(axis=1), 1))

ddof=1是样本标准差，小样本场景下和ddof=0的差异会放大到百分位级别。M.T那一步的转置很容易漏，写成Z_w @ M会直接报形状不匹配，但如果你恰好用了方阵样本，它会静默跑出一个错的结果，这属于后悔药都买不到的那种。

3.4 结果要过两道闸：合计是否接近 100、方向是否符合化学预期

反推出来的成分表不能直接交。第一道闸看每行合计，应当落在 100 附近，偏离超过 5 说明主子空间截断太狠，把var_ratio提到 0.98 试试。第二道闸看方向：铅钡玻璃风化前 SiO₂ 应当在 80% 上下、PbO 很低，风化后 SiO₂ 掉到 36% 左右、PbO 涨到 47% 左右，如果预测出来的未风化值把这个方向搞反了，问题一定出在特征顺序或负载矩阵符号上。

注意：因子分析给出的是统计意义上的期望成分，不是某个样品真实的历史配方。对严重风化样品，预测值的置信区间会宽到没什么参考价值，这一点必须写进模型评价，不能只报点估计。

4. 亚类划分与未知样本归类：一维区间比多维聚类更稳

4.1 四个子类必须分开做 PCA

表单 2 天然分成高钾风化、高钾未风化、铅钡风化、铅钡未风化四块。把它们混在一起做主成分分析，第一主成分会被“风化与否”这个巨大差异吃掉，剩下的主成分里几乎看不到亚类结构。分开做之后，铅钡未风化类里 SiO₂、Na₂O、K₂O、CaO、MgO、Al₂O₃ 六项的累计贡献率能到 88% 左右，高钾未风化类五项到 88% 左右，铅钡风化类五项 87% 左右，高钾风化类前三项就有 95% 以上。样本量最少的那个子类只有 7 条记录，PCA 之前必须先砍特征，否则协方差矩阵直接退化。

4.2 去极值取区间：为什么不能直接用 min 和 max

“数值落在区间里或者靠近区间就属于这一类”，这个判据要成立，区间边界必须稳。直接用最小最大值，一个离群点就能把区间撑开一倍。这套资源的处理是排序后去掉首尾各一个样本，再取剩余数据的上下界。

def core_interval(x, trim=1): """去极值后的核心区间；样本量 >= 6 才去，小样本去了就没了""" x = np.sort(np.asarray(x, dtype=float)) if len(x) >= 6 and trim > 0: x = x[trim:-trim] return float(x.min()), float(x.max()) for name, sub in df2.groupby(["类型", "判风化"]): lo, hi = core_interval(sub["SiO2"]) print(name, "SiO2 核心区间", round(lo, 2), round(hi, 2))

trim参数取 1 是经验值。取 2 在 11 个样本的子类里会砍掉近两成数据，区间收得太紧，表单 3 的新样品很容易全被判成“不在区间内”。反过来取 0 就等于没做清洗，一个 96.77 的高硅异常点能让高钾风化类的下界从 92.35 掉到 92 以下。

4.3 最短距离法的两个聚点，落到一维就是找最大间隙

亚类划分选的是数值差异最大的那个成分，实际数据里就是 SiO₂。离散度最大的一维数据做二分类，最短距离法聚类得到的两个聚点，切分点取两者均值。用 KMeans 设定n_clusters=2在这一维上等价，只是要注意随机初始化。

from sklearn.cluster import KMeans def two_centers(x, seed=42): x = np.asarray(x, dtype=float).reshape(-1, 1) km = KMeans(n_clusters=2, n_init=20, random_state=seed).fit(x) c = np.sort(km.cluster_centers_.ravel()) return c, float(c.mean()) for name, sub in df2.groupby(["类型", "判风化"]): c, cut = two_centers(sub["SiO2"]) print(name, "聚点", np.round(c, 2), "切分点", round(cut, 2))

原始数据上跑出来大致是：铅钡未风化聚点 64.82 / 52.49，切分 58.66；高钾未风化 71.22 / 62.77，切分 67.00；铅钡风化 33.96 / 15.63，切分 24.80；高钾风化 94.66 / 92.57，切分 93.62。n_init=20不是摆设，一维两簇有时会陷进局部最优，多初始化几次结果更稳。

合理性检验看两段里样本比例是否失衡，某一类里一段占了九成、另一段只有一两个，说明这个成分不适合做亚类划分，要换成分重来。

4.4 3% 扰动下的敏感性，口径必须写死

题目要求给数据加 3% 随机扰动看分类结果变不变。这句话里有两个自由度必须自己定死：扰动是均匀分布还是正态分布，以及重采样几次。这套资源的写法是均匀分布、固定种子、重复 200 轮，统计“样品-轮次”维度的类别翻转比例。

rng = np.random.default_rng(20220401) def label_sio2(v, lo, hi, cut, tol=0.05): """0=判不进该类，1=亚类一，2=亚类二；tol 是“靠近区间”的容差""" if v < lo * (1 - tol) or v > hi * (1 + tol): return 0 return 1 if v < cut else 2 def flip_rate(x, lo, hi, cut, rounds=200, pct=0.03, tol=0.05): x = np.asarray(x, dtype=float) base = np.array([label_sio2(v, lo, hi, cut, tol) for v in x]) flip = 0 for _ in range(rounds): xp = x * (1 + rng.uniform(-pct, pct, size=x.shape)) new = np.array([label_sio2(v, lo, hi, cut, tol) for v in xp]) flip += int((new != base).sum()) return flip / (rounds * len(x))

tol=0.05把“靠近区间”这句话量化了，不加这个参数，边界附近的样品在扰动下必然翻转，误差会虚高到 30% 以上。实际跑下来，四类里靠近切分点的样品翻转率明显更高，特别是亚类一样本数只有两三个的那一类，扰动后整段搬家的概率不小。这个结论要如实写，别为了好看把tol调到 0.2。

5. 避坑与排查：五个真跑起来才会撞上的问题

5.1 成分合计不等于 100，别一上来就强行归一化

现象：清洗完打印成分合计，能看到 96.8、103.2 这种值，第一反应是统一除以合计乘 100。原因：原始检测里存在未检出记 0、灼烧损失未单列、不同实验室仪器精度差异，强行归一会把误差平均摊到所有成分上，把本来真实的异常信号抹平。解决：保留原始比例不动，把成分合计单独作为一列质检指标，偏差超过 3 个百分点的样本在文档里标注出来，建模时用原值，只在计算风化系数前对分母做零值保护。这套资源的系数定义全部基于原始百分比，没有做闭合处理，就是为了保住这个偏差信息。

5.2 表单 3 直接 merge 出大片 NaN

现象：把表单 3 按编号 merge 上表单 1 的信息，类型列七成是空的。原因：表单 3 本来就是待鉴别的未知样品，表单 1 里根本没有对应编号，how="left"之后自然是 NaN，这不是代码错。解决：表单 3 只 merge 纹饰、颜色这类可能存在的辅助列，类型列留空等模型去填；更稳的做法是给表单 3 单独建一张特征表，只保留OXIDES和成分合计，用第二问得到的区间表去匹配，别和训练表混在一起。我第一次跑的时候没注意，后面做敏感性分析才发现测试集里混进了训练集的行，白白多跑了两小时。

5.3 对负载矩阵直接求逆，报错或者结果全乱

现象：np.linalg.inv(A)抛LinAlgError: Last 2 dimensions of the array must be square，或者有人为了绕过报错先A @ A.T再求逆，结果维度对不上。原因：负载矩阵是 p×k 的长方阵，k 通常小于 p，数学上不存在逆矩阵。解决：用pinv(A)或者显式写(A.T @ A)^-1 @ A.T，两者在 A 列满秩时等价。顺带提醒，如果某次跑出来恰好 k 等于 p（累计贡献率取 1.0 时会发生），inv能跑通但结果是过拟合的，等于没降维，务必检查 k 的取值再下结论。

5.4 先全量标准化再切分组，敏感性分析会假性乐观

现象：扰动实验跑出来翻转率只有 3%，看着模型极其鲁棒，换一批数据立刻崩。原因：标准化用了全量样本的均值和方差，测试样本的分布信息提前进了训练流程，加上扰动只是在小范围内缩放，标准化后几乎被抵消掉。解决：标准化参数一律从未风化组或训练组估计，测试组只做transform不做fit。判断有没有踩这个坑很简单，把训练组和测试组的均值打印出来对比，差得离谱就说明泄漏了。

5.5 扰动实验复现不了，两次结果差十几个百分点

现象：同一个脚本跑两遍，翻转率一次 18% 一次 27%。原因：没有固定随机种子，或者扰动幅度写成了np.random.randn却没乘标准差，实际扰动尺度随样本量级变化。解决：统一用np.random.default_rng(seed)，扰动写成乘性x * (1 + uniform(-pct, pct))而不是加性，保证 3% 是相对幅度；把rounds、pct、tol三个参数写进论文的参数表，别人才能复现。这个习惯从那以后我一直保留，凡是带随机的实验，种子和分布类型必须和结果一起交。

6. 灰色关联度与鲁棒性验证：让结论站得住的三个细节

第四问比较不同类别内部化学成分的关联关系，用灰色关联度比皮尔逊相关系数合适，样本量小、分布不明、还带量纲差异，灰色关联正好不吃这些假设。但有三个细节决定了结论能不能用。

第一是无量纲化方式要统一。均值化处理比初值化稳，因为初值化对第一个样本极度敏感，而这批数据里第一个样本恰好风化程度偏重。第二是分辨系数 ρ 取 0.5，题目里给了参考区间，取 0.5 是通行做法，改到 0.3 会让关联系数拉开差距但排序基本不变，写论文时要把这个取值说清楚。第三是必须按类别分别算再横向比，不能把四类数据拼成一张表一起算，否则组间差异会淹没组内结构。

def grey_relational(X, ref, rho=0.5): """X: (n, m) 比较数列；ref: (n,) 参考数列；返回每个成分的关联度""" X = np.asarray(X, dtype=float) ref = np.asarray(ref, dtype=float) Xs = X / X.mean(axis=0) # 均值化无量纲 rs = ref / ref.mean() d = np.abs(rs[:, None] - Xs) # 绝对差序列 dmin, dmax = d.min(), d.max() xi = (dmin + rho * dmax) / (d + rho * dmax) return xi.mean(axis=0) for name, sub in df2.groupby(["类型", "判风化"]): X = sub[["Na2O", "K2O", "CaO", "MgO", "Al2O3"]].values ref = sub["SiO2"].values g = grey_relational(X, ref) order = np.argsort(g)[::-1] print(name, [(c, round(g[i], 3)) for i, c in enumerate(np.array(["Na2O","K2O","CaO","MgO","Al2O3"])[order])])

铅钡类的排序稳定在 SiO₂ > CaO > K₂O > Al₂O₃ > MgO > Na₂O，四类之间差异不大，说明这些成分在风化过程中同步变化。高钾类的后两位在 Na₂O 和 MgO 之间来回浮动，原因是高钾风化组只有 7 条样本，参考数列一短，关联系数的均值就容易被单个样本带偏。

所以这一问真正该写的不是“我们发现了某类玻璃里某成分最重要”，而是“在样本量这么小的条件下，排序对样本子集敏感，只有 SiO₂ 的领先位置是稳的”。把不确定性交代清楚，比硬凑一个漂亮排序可信得多。我后来的习惯是，凡是灰色关联这类对小样本敏感的方法，跑完都要把样本随机抽掉两成再跑一遍，两次排序一致才敢往结论里写。希望帮到你。

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