古代玻璃文物成分分析建模复现：风化系数、PCA与灰色关联度

简介：面向参加全国大学生数学建模竞赛、希望系统复盘C题解法的学生与指导教师，这份资料围绕古代玻璃文物的成分分析与鉴别展开，覆盖风化效应、类型判别、未风化成分预测及化学成分关联性分析等核心问题。包内仅有1个PDF文件，压缩包约3.55MB，以文档形式集中呈现赛题解读、建模思路、数据处理过程和结果分析，便于按问题顺序查阅。已有928人学习下载，说明其在赛题复盘与建模方法参考方面具有一定关注度。读者可从中了解数据预处理与探索、主成分分析降维、因子分析转换、敏感性检验、灰色关联度分析以及分类预测与模型评估的具体做法，并借助完整赛题方案梳理从特征选择到模型解释的链路，适合用作竞赛复盘、方法对照和论文写作参考。

1. 从一份对不上号的表格说起

先说场景。你拿到一套数学建模赛题 C 题的交付包，里面有论文、三张 Excel 附件、若干 Python 脚本，还有几段 MATLAB 代码。你想复现论文表 7 里那列“未风化成分预测值”，跑完脚本发现数字和纸上印的差了一位小数，于是开始怀疑自己环境和数据是不是搞错了。这套资源的价值就在这里：它把古代玻璃文物的成分分析拆成四条能独立验证的链路——风化系数标定、主成分降维加因子分析反推、区间划分做亚类、灰色关联度做成分关联。适合想完整跑通一次建模流程的新手，也适合准备把赛题代码改成自己项目模板的老手。

2. 技术路线与文件构成：四个问题怎么落到一份可跑的交付物上

2.1 四个问题其实是四条独立技术路线

把论文目录翻开看，四个问题并不是一条线串下来的，而是四套方法各管一段，互相之间只通过“特征列”和“中间结果”耦合。先把这个映射关系理清，后面读代码就不会迷路。

问题一处理的是“风化前后成分怎么变”。输入是表单 1 的文物属性和表单 2 的成分含量，先做相关性分析和卡方检验，确认风化与类型、纹饰、颜色有没有统计关联；再对风化样本做主成分降维，用因子分析的载荷矩阵求逆，反推出未风化时的成分向量。产物是一张“预测未风化成分”表，论文里对应表 7。

问题二处理的是“同类玻璃内部还能不能再分”。做法是先按类型和风化状态把表单 2 切成四组，每组做主成分分析挑出贡献率高的几个成分，去掉最大最小值后取剩余极值构成区间；再挑差异最大的二氧化硅做最短距离法聚类，得到两个聚点，取中点作为亚类分界。产物是四组区间表和四个分界点。

问题三是最省事的一问。把表单 3 的未知样品逐条代入第二问的区间表，落在区间内或靠近区间就判为该类；然后给数据乘一个 3% 的随机扰动，看判别结果翻不翻，翻得多就说明模型敏感。产物是一张类别归属表加一组扰动前后的对比。

问题四处理“成分之间谁跟谁更相关”。在第二问的四个亚类基础上，对每个类别用灰色关联度分析，把某个成分当参考序列、其余成分当比较序列，算出关联度并排序。产物是四条排序链，论文里给出的结论是二氧化硅始终排第一，其余成分的相对位置在铅钡和高钾之间有差别。

提示：四条路线共享同一份清洗后的宽表，所以第 3 章的预处理代码是后面所有分析的地基，跑不通它，后面全是空中楼阁。

2.2 文件构成与目录组织

压缩包解开之后通常是“论文 + 数据 + 脚本”三块。论文是 PDF 或 Word，里面表 1 到表 7、图 1 到图 16 是核对脚本输出的唯一基准；数据是三个 xlsx，分别对应文物基本信息、已分类样品的成分含量、待鉴别的未知样品；脚本则按问题编号分散存放，Python 和 MATLAB 混着用。

实际包里的命名可能和下面不一样，按“问题编号 + 方法名”去找基本不会错：

competition-c-glass/ ├── paper/ │ └── main.pdf # 论文正文，表 1~7、图 1~16 是核对基准 ├── data/ │ ├── form1.xlsx # 类型、纹饰、颜色、表面风化 │ ├── form2.xlsx # 已分类样品的 14 项化学成分含量 │ └── form3.xlsx # 待鉴别的未知类别样品 ├── code/ │ ├── q1_factor.py # 相关性检验 + PCA + 因子分析反推 │ ├── q2_pca_cluster.py # 四组区间生成 + 最短距离法聚类 │ ├── q3_predict.py # 表单 3 代入判别 + 3% 扰动 │ └── q4_grey.m # 灰色关联度，MATLAB 实现 └── README.md

这里有个现实问题：预处理那一段，论文里写的是“借助 Excel 功能对表单 2 进行预处理”，也就是说加列、拼表、排序这些活是在表格软件里手工做的，脚本里并不包含。这就导致你直接跑 q1_factor.py，它读进来的数据可能已经被人改过一轮，和你手上的原始 xlsx 不是同一份。我一般的做法是先把 Excel 里那几步用 pandas 重写一遍，让从原始表单到分析宽表的每一步都能脚本复现，再往下接建模代码。

2.3 技术栈与运行环境

这套资源的技术栈是典型的“三件套拼接”：Excel 做手工预处理和透视，SPSS 做卡方检验和相关性分析，Python 做降维、聚类、区间判别，MATLAB 做灰色关联度。论文表 3 到表 5 的相关系数、P 值、Phi、Cramér's V、lambda 都是 SPSS 的输出格式，表 6 的转换系数矩阵和表 7 的预测值来自 Python，图 13 到图 16 的关联系数来自 MATLAB。

复现时不必把三套工具都装上。SPSS 那部分用 scipy 的chi2_contingency加手写的效应量公式就能替代，MATLAB 的灰色关联度用二十行 numpy 也能写完。真正需要装的环境是 Python 3.8 以上：

python -m venv .venv source .venv/bin/activate # Windows 用 .venv\Scripts\activate pip install pandas numpy scipy scikit-learn matplotlib openpyxl

pandas和openpyxl负责读 xlsx，scikit-learn提供PCA和FactorAnalysis，scipy里有卡方检验和层次聚类，matplotlib用来复刻论文里的饼图和分布图。版本上不用太纠结，这几个库的 API 近几年没大变，反倒是 Excel 文件里的列名带全角括号和单位这种事更容易让脚本在第一行就崩掉。

3. 数据层与问题一：从三张表单到因子分析反推未风化成分

3.1 三张表单怎么拼成一张分析表

表单 2 是主表，一行一个样品，列是十四种化学成分的百分比；表单 1 是属性表，一行一个文物，列是类型、纹饰、颜色、表面风化。两边靠“文物采样点”关联，但这个字段并不干净，里面混着“30 部位 1”“42 未风化点 1”这类带后缀的写法，直接 merge 会大量匹配不上。

处理办法是用正则把编号抽出来当连接键，同时加一个validate参数兜底，防止一对多把行数悄悄撑大：

import pandas as pd import numpy as np form1 = pd.read_excel("data/form1.xlsx", sheet_name=0) form2 = pd.read_excel("data/form2.xlsx", sheet_name=0) form3 = pd.read_excel("data/form3.xlsx", sheet_name=0) def clean_columns(df): # 列头一般长这样：二氧化硅(SiO2)，剥掉括号和单位再对齐 return df.rename(columns=lambda c: str(c).split("(")[0].strip()) form2, form3 = clean_columns(form2), clean_columns(form3) # 从"30 部位 1"这类字符串里抽出第一段数字作为编号 form2["编号"] = form2["文物采样点"].astype(str).str.extract(r"(\d+)").astype(int) form1["编号"] = form1["文物采样点"].astype(str).str.extract(r"(\d+)").astype(int) merged = form2.merge( form1[["编号", "类型", "纹饰", "颜色", "表面风化"]], on="编号", how="left", validate="many_to_one" # 右侧必须唯一，否则直接报错 ) chem_cols = [c for c in form2.columns if c not in ("文物采样点", "编号")] merged["成分合计"] = merged[chem_cols].sum(axis=1) print(merged["成分合计"].describe()) # 明显偏离 100 就说明有缺失或单位问题

这段代码的关键在三个地方。validate="many_to_one"是一道保险，如果表单 1 里同一个编号出现了两行，pandas 会立刻抛异常，而不是给你一份行数翻倍的宽表让你后面慢慢查。正则(\d+)只取第一段数字，正好把“部位”“未风化点”这些后缀丢掉。最后那列“成分合计”不是建模用的特征，是体检指标，正常样品的十四项加起来应该接近 100，如果某行只有 60 多，八成是某个成分缺失被读成了 NaN。

3.2 风化系数：把“风化”从定性描述变成可计算阈值

论文里最有价值的一个构造是风化系数。原始数据里“是否风化”是人工标注的标签，但标注有主观性，论文里也提到样本 20、48、30 部位 1 三个样本的标签和系数对不上，最后以系数为准做了手工调整。所以这个系数本质上是把“风化程度”变成一个连续可比的数值，阈值从数据分布的断档里读出来，而不是拍脑袋定。

两类玻璃的风化机理是相反的：铅钡玻璃风化后铅、钡富集而硅流失，高钾玻璃风化后钾、钙、镁这些助熔剂成分流失而硅相对富集。所以系数公式必须分开写。我按论文表 1 和表 2 的数值反推验证过，铅钡用(PbO + BaO) / SiO2，高钾用SiO2 / (K2O + CaO + MgO + Al2O3 + Fe2O3)，算出来和论文里印的风化系数能对上小数点后三位。

def add_weathering_index(df: pd.DataFrame) -> pd.DataFrame: df = df.copy() eps = 1e-9 # 防止除零，成分含量可能真的是 0 pb = df["类型"] == "铅钡" k = df["类型"] == "高钾" # 铅钡：铅钡富集、硅流失，比值随风化程度单调上升 df.loc[pb, "风化系数"] = ( (df.loc[pb, "PbO"] + df.loc[pb, "BaO"]) / (df.loc[pb, "SiO2"] + eps) ) # 高钾：助熔剂流失、硅富集，硅除以助熔剂之和随风化程度上升 flux = df.loc[k, ["K2O", "CaO", "MgO", "Al2O3", "Fe2O3"]].sum(axis=1) df.loc[k, "风化系数"] = df.loc[k, "SiO2"] / (flux + eps) # 阈值来自排序后系数序列的断档：铅钡在 0.88~1.30 之间跳，取 1.0 # 高钾在 7.23~12.41 之间跳，取 10.0 df["判定风化"] = np.where( (pb & (df["风化系数"] > 1.0)) | (k & (df["风化系数"] > 10.0)), 1, 0 ) return df merged = add_weathering_index(merged) print(pd.crosstab(merged["表面风化"], merged["判定风化"]))

参数上要留意两点。eps取 1e-9 是纯粹的技术兜底，不参与业务含义；两个阈值 1.0 和 10.0 不是超参数，是数据本身给出的分界，换了数据集必须重新画散点图找断档，不能照抄。跑完最后那行交叉表，如果人工标签和判定结果的交叉表对角线占绝大多数，说明系数构造是对的；如果错得离谱，先查类型字段有没有空值，pb和k两个布尔索引都为 False 时那一行算不出系数。

3.3 相关性检验：P 值和效应量要一起读

论文表 3 给的是相关系数矩阵，表 4 是卡方检验，表 5 是效应量。这三张表经常被误读成“P 小于 0.05 就说明强相关”，其实 P 值只回答“有没有关联”，效应量才回答“关联有多强”。纹饰对表面风化的 Phi 是 0.561，颜色是 0.707，类型是 0.272，这三个数放在一起看，才能说颜色和风化的关联最紧。

用 scipy 复刻这套检验很直接，卡方统计量和 P 值现成，效应量要自己算：

from scipy.stats import chi2_contingency def cramers_v(table: pd.DataFrame) -> float: chi2 = chi2_contingency(table)[0] n = table.values.sum() r, k = table.shape return float(np.sqrt(chi2 / (n * (min(r, k) - 1)))) # 取值 0~1 for col in ["纹饰", "颜色", "类型"]: table = pd.crosstab(merged["表面风化"], merged[col]) chi2, p, dof, _ = chi2_contingency(table) print(f"{col}: chi2={chi2:.3f}, p={p:.4f}, V={cramers_v(table):.3f}")

cramers_v的分母里min(r, k) - 1是校正项，列联表不是方阵时必须用它，否则算出来的值会超过 1。chi2_contingency默认带 Yates 校正，2x2 表时和论文里 SPSS 输出的“校正 X²”对得上。这段代码跑完如果 P 值都在 0.05 以下，说明风化状态确实和这三个属性有关联，后面的建模才有依据；如果某个属性的 P 值很大，那在问题一里就不该把它当主要解释变量。

3.4 PCA 降维与因子分析反推未风化成分

问题一的核心诉求是：已知一个风化样品的成分，反推它没风化时长什么样。论文的思路是把成分向量拆成“公共因子 + 特殊因子”，公共因子就是未风化成分，载荷矩阵 A 描述从公共因子到观测成分的投影，那么只要 A 可逆，就能从观测值反解公共因子。

实现上有两个坎。第一，A 必须是方阵才能直接求逆，所以主成分个数要选成和参与建模的变量数一致，论文里选 8 个主成分、保留 97% 以上方差，就是为了凑这个方阵。第二，求逆之前必须检查条件数，成分之间高度共线时 A 接近奇异，反推出来的值会出现负数和离谱的数量级。

from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA, FactorAnalysis cols = ["SiO2","Na2O","K2O","CaO","MgO","Al2O3","Fe2O3","CuO","PbO","BaO","P2O5","SrO"] # 只用未风化样本估计"正常成分"的协方差结构，这是整个反推的基准 sub = merged.query("类型 == '铅钡' and 判定风化 == 0")[cols].fillna(0) scaler = StandardScaler() Xs = scaler.fit_transform(sub) pca = PCA(n_components=0.97, svd_solver="full") # 累计方差 97% 决定保留几个 Xs_pca = pca.fit_transform(Xs) k = pca.n_components_ fa = FactorAnalysis(n_components=k, random_state=0) fa.fit(Xs) A = fa.components_.T # p x k 载荷矩阵 print("主成分个数:", k, "条件数:", np.linalg.cond(A))

n_components=0.97是让 sklearn 自动按累计方差贡献率选个数，比手写循环判断省事；FactorAnalysis的components_形状是k x p，转置之后才是论文里那个p x k的载荷矩阵。条件数这行打印出来如果在几百以内，求逆是安全的；上万就说明变量之间有强共线，得先砍掉几个成分再算。

反推那一步，论文是直接把 A 求逆再乘中心化后的观测向量。实际写的时候我建议用pinv，它在方阵且满秩时和inv结果一致，矩阵病态时退化成最小二乘解，不会直接抛异常：

def reconstruct_unweathered(x_obs, mu, A): """x_obs: 单个风化样品的成分向量; mu: 未风化样本均值; A: p x k 载荷矩阵""" f = np.linalg.pinv(A) @ (x_obs - mu) # 反推公共因子 return np.clip(f, 0, None) # 成分含量不可能为负，裁掉 mu = sub.mean().values weathered = merged.query("类型 == '铅钡' and 判定风化 == 1")[cols].fillna(0).values pred = np.array([reconstruct_unweathered(x, mu, A) for x in weathered]) print(pd.DataFrame(pred, columns=cols).round(2).head())

mu必须用未风化样本的均值，不能用全体样本，否则中心化的基准就偏了。最后的clip是血泪经验：线性反推不保证非负，不加这一步，表 7 里会出现-3.21这种含量，答辩时很难解释。裁完之后再和论文表 7 对照，允许小数点后一位的差异，如果整体偏移很大，回头检查标准化是不是在未风化子集上 fit 的。

4. 问题二到问题四：区间分类、3% 扰动与灰色关联度排序

4.1 主成分筛选与成分区间生成

问题二的第一步是把表单 2 按“类型 × 风化状态”切成四组，每组单独做一次主成分分析，挑出贡献率高的几个成分。论文给出的结果是铅钡未风化取六个成分、累计贡献率 88.89%，高钾未风化取五个、铅钡风化取五个、高钾风化取三个，累计贡献率都在 87% 以上。

第二步是生成区间。做法不是取均值的置信区间，而是把数据排序后去掉最大最小各一个，取剩余部分的极值当上下界。这个做法比较糙，但对小样本很实用，能天然躲开单个异常点的干扰：

def robust_interval(s: pd.Series, trim: int = 1) -> tuple: v = np.sort(pd.Series(s).dropna().values) if len(v) <= 2 * trim: # 样本太少就不裁，否则区间为空 return float(v.min()), float(v.max()) v = v[trim:-trim] return float(v.min()), float(v.max()) groups = { "铅钡-未风化": (merged.query("类型=='铅钡' and 判定风化==0"), ["SiO2","Na2O","K2O","CaO","MgO","Al2O3"]), "高钾-未风化": (merged.query("类型=='高钾' and 判定风化==0"), ["SiO2","Na2O","K2O","CaO","MgO"]), "铅钡-风化": (merged.query("类型=='铅钡' and 判定风化==1"), ["SiO2","Na2O","K2O","CaO","MgO"]), "高钾-风化": (merged.query("类型=='高钾' and 判定风化==1"), ["SiO2","Na2O","K2O"]), } intervals = {name: {c: robust_interval(d[c]) for c in cols_g} for name, (d, cols_g) in groups.items()}

trim=1是去掉一头一尾各一个样本，样本量二三十个时这个力度刚好；样本上百时可以考虑改成按分位数裁。跑完拿结果和论文里的区间对一下，铅钡未风化的 SiO2 应该是[50.61, 69.71]这个量级，高钾风化的是[92.35, 95.02]——风化高钾玻璃的硅含量极高，这是它最显著的特征，也是后面分类能分得开的原因。

4.2 最短距离法聚类出亚类分界点

亚类划分只用了二氧化硅一个维度，理由是四组数据里 SiO2 的数值差异最大、区分度最好。方法是最短距离法，也就是层次聚类里的single链接，把样本并成两类之后取两个聚点的中点当分界。

from scipy.cluster.hierarchy import linkage, fcluster def split_by_single_linkage(values, k: int = 2): v = np.asarray(values, dtype=float).reshape(-1, 1) z = linkage(v, method="single") # 最短距离法 labels = fcluster(z, t=k, criterion="maxclust") centers = np.array([v[labels == i].mean() for i in range(1, k + 1)]) centers.sort() return centers, float(centers.mean()) # 聚点 + 分界点 centers, boundary = split_by_single_linkage( groups["铅钡-未风化"][0]["SiO2"].dropna() ) print("聚点:", centers, "分界:", round(boundary, 2))

method="single"就是最短距离法，两类之间的距离定义为最近两个样本的距离；criterion="maxclust"表示按“最多分几类”来切，不是按距离阈值切。这个组合在一维数据上等价于“找最大的那个间隙”，所以对离群点极其敏感——如果不去掉极值，一个 95 的异常点就能把分界点拽偏十几。论文里算出来的分界点分别是铅钡未风化 58.70、高钾未风化 67.00、铅钡风化 24.80、高钾风化 93.62，正好是两个聚点的算术平均。

分完之后论文用饼图看两类样本的分布比例，判断有没有偏得离谱。这个检查比看轮廓系数直观，如果一类只占 5%，说明聚类结果没意义，得换个成分或者换个 k 值再试。

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

表单 3 是未知类别的样品，判别逻辑就是把每行代入四组区间，全成分落在区间内就算命中。论文里说“数值在这个区间里面或者靠近这个区间就属于这个种类”，这个“靠近”没有给严格定义，实现时一般给个外扩余量margin：

def classify_by_intervals(row, intervals: dict, margin: float = 0.0): hits = [] for name, rule in intervals.items(): if all(lo - margin <= row[c] <= hi + margin for c, (lo, hi) in rule.items()): hits.append(name) return hits[0] if hits else "unknown"

margin取 0 时是严格落在区间内，命中率会偏低；放大到区间宽度的 5% 到 10% 能提高召回，但会带来误判。实际做的时候我一般先设 0 跑一遍，看有多少条判成 unknown，再决定要不要放宽。四组区间里高钾风化的成分最少（只有三项），所以它的判别条件最宽松，容易被误命中，排序上要把它放在最后判断。

敏感性分析是给数据乘一个 3% 的随机扰动，再看判别结果翻不翻。用乘性扰动而不是加性，是因为成分是比例数据，加一个绝对值对不同量级的成分影响完全不成比例：

def perturbation_test(X, classify_fn, ratio: float = 0.03, repeat: int = 30, seed: int = 2022): rng = np.random.default_rng(seed) base = [classify_fn(r) for r in X] flips = [] for _ in range(repeat): noisy = X * rng.normal(1.0, ratio, size=X.shape) now = [classify_fn(r) for r in noisy] flips.append(np.mean([a != b for a, b in zip(base, now)])) return float(np.mean(flips)), float(np.std(flips))

seed必须固定，否则每次跑出来波动率都不一样，没法写进报告；repeat=30是个折中，跑一百次结果更稳但没必要。返回值报“均值 ± 标准差”而不是只报一个数，因为扰动本身是随机的，单次结果说明不了问题。论文里表单 3 的判别在扰动后基本不翻，说明分类边界离样本点比较远；而问题二里对表单 2 重新分类时误差超过 20%，这两处结论方向不一致，第 5 章会专门说这个。

4.4 灰色关联度：排序与类别差异

问题四要回答的是“同类玻璃内部，哪些成分的变化趋势更同步”。灰色关联度的做法是选一个成分当参考序列，其余成分当比较序列，先无量纲化，再逐点算关联系数，最后取均值当关联度。

def grey_degree(ref, comp, rho: float = 0.5): """ref: (n,) 参考序列; comp: (n, m) 比较序列; 返回每个比较序列的关联度""" ref = np.asarray(ref, float) comp = np.asarray(comp, float) ref_n = ref / (ref.mean() + 1e-12) # 均值化无量纲 comp_n = comp / (comp.mean(axis=0) + 1e-12) delta = np.abs(ref_n[:, None] - comp_n) dmin, dmax = delta.min(), delta.max() xi = (dmin + rho * dmax) / (delta + rho * dmax) # 关联系数 return xi.mean(axis=0) # 关联度

rho是分辨系数，取值在 0 到 1 之间，越小分辨力越强，论文取 0.5 是通用做法。均值化那一步不能省，成分含量的量级差着两个数量级，不归一化的话绝对值大的成分会主导整个计算。返回的关联度排序，论文给的结论是铅钡未风化为 SiO2 > CaO > K2O > Al2O3 > MgO > Na2O，高钾那边 SiO2 同样排第一但后面几位顺序不同——这个差异就是“两类玻璃成分关联关系不同”的证据。

有个细节值得注意：灰色关联度算的是序列之间的形状相似度，不是统计相关，所以它对样本量要求低，七八条样本也能算。但也正因为如此，样本顺序会影响结果，把数据表按不同列排序后

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