简介2022年全国大学生数学建模竞赛C题的完整解题资料包面向备赛学生与数学建模爱好者围绕古代玻璃文物成分分析与鉴别问题系统梳理了数据预处理、数据探索、特征工程、主成分分析与因子分析等关键建模方法。资源为单个PDF文档大小3.55MB虽以文档形式交付但内容涵盖完整建模思路、核心代码实现、结果图表与结论便于对照学习。包内重点展示高钾玻璃与铅钡玻璃的风化系数定义、亚类划分区间、敏感性分析与灰色关联度分析过程并给出四问的详细解答与模型解释适合需要参考优秀论文结构、理解降维与关联分析应用的读者目前已有927人学习下载。整体内容紧凑精炼可直接作为赛题复盘、论文写作或算法实战的参考资料。1. 项目速览数学建模C题如何给古代玻璃“去风化”2022年全国大学生数学建模C题的核心是给一批风化了几百上千年的古代玻璃文物做“成分还原”和“自动鉴别”。这道题乍一看是考古学问题实际拆开全是数据活儿三张表单、十几个化学成分指标、一堆风化与未风化的标记要用主成分分析、因子分析、最短距离聚类、灰色关联度把这堆数字变成可以解释的结论。这个资源就是整套论文与可运行代码的打包论文负责讲清建模思路代码负责把四个问题的结果完整复现。适合两类人一是备战数学建模竞赛、想找一份完整C题范式的同学二是对成分数据做降维、聚类、敏感性分析套路感兴趣的开发者。2. 数据底子与预处理三张表单与两类风化系数2.1 三张表单的字段关系与合并思路这道题的数据核心是三张表先把结构理清后面所有模型都建立在这三张表之上表单1文物采样点的档案包含样品编号、玻璃类型高钾玻璃/铅钡玻璃、是否风化、颜色、纹饰表单2每个采样点的化学成分检测结果从B列到O列共十几个成分包括二氧化硅SiO2、氧化钠Na2O、氧化钾K2O、氧化钙CaO、氧化镁MgO、氧化铝Al2O3、氧化铁Fe2O3、氧化铜CuO、氧化铅PbO、氧化钡BaO、五氧化二磷P2O5、氧化锶SrO等表单3一批未知类别的玻璃文物成分数据待判别所属类型。原始文档里用Excel做了大量预处理在B列到O列后面增加一个合计列P验证每行成分百分比加起来是否接近100%从A列抽取样品编号放入Q列再把表单1中的“玻璃类型、是否风化、颜色、纹饰”四项信息按样品编号关联到表单2形成R、S、T、U四列。这一步本质上是两张表按主键做left join用pandas实现是更常见的做法import pandas as pd info pd.read_excel(表单1.xlsx) chem pd.read_excel(表单2.xlsx) # 从A列抽出纯数字编号兼容“XX部位1”这类复合编号 chem[sample_id] chem[A列].astype(str).str.extract(r(\d)) # 计算各成分合计做数据质量校验 chem[total] chem.iloc[:, 1:15].sum(axis1) merged chem.merge(info, left_onsample_id, right_on样品编号, howleft)这里的sample_id提取是关键表单2里同一件文物会有“部位1”“部位2”“严重风化点”等多个采样点直接用原始文本做关联会失败先抽编号再关联是这步的通用做法。total列用来检查成分合计是否在100%附近正常玻璃样本的检测值合计应在100%上下偏差超过5%的样本多半存在复制时的乱码或串列问题要单独标记出来。howleft保证化学表里每一行的采样数据都保留下来避免有效样本被内连接悄悄过滤掉。2.2 两类风化系数的定义与阈值确定预处理的关键一步是定义风化系数。论文没有用统一公式而是给两类玻璃分别定义因为这两类玻璃的风化规律恰好相反铅钡玻璃风化系数 (PbO BaO) / SiO2 高钾玻璃风化系数 SiO2 / (K2O CaO MgO Al2O3 Fe2O3 CuO)从物理意义上看铅钡玻璃风化后SiO2占比显著流失氧化铅、氧化钡相对富集比值变大高钾玻璃风化后SiO2占比显著提升钾、钙、镁、铝等助熔剂氧化物流失比值同样变大。两边风化系数都是“越大越风化”但分子分母的含义相反这是这个设计最精巧的地方。再看阈值是怎么定的。论文对“铅钡玻璃-A纹”这组样本按风化系数排序发现0.879到1.3之间存在明显断层于是把1.0作为判别阈值对高钾玻璃7.23到12.41之间存在断层取10.0作为阈值。这种取法非常直观本质就是找排序后相邻值的最大空档用代码可以写成import numpy as np def split_by_gap(values): s np.sort(values) gaps np.diff(s) idx np.argmax(gaps) threshold (s[idx] s[idx1]) / 2 return threshold, s, gaps参数values传入同一类型样本的风化系数数组返回的是断层中点作为判别阈值以及排序后的序列和相邻差值。用这个函数对三类样本分别计算得到的结果与论文里手工观察的阈值一致。实际使用中我一般会先画一张排序后的折线图肉眼看断层在哪儿再用这个函数辅助定位避免把正常波动误判成阈值点。2.3 数据探索挖出的三条规律对预处理后的数据做分组透视原文得出三个值得记住的观察纹饰B的玻璃文物全部是高钾玻璃、全部风化、颜色全是蓝绿色纹饰A的高钾玻璃全部无风化且为蓝绿色铅钡玻璃的风化与无风化比例是11:5纹饰C的高钾玻璃全部无风化铅钡玻璃的风化与无风化比例是17:73。由此推出两个结论同等条件同纹饰意味着年代和地域可能相近下铅钡玻璃比高钾玻璃更容易风化纹饰B代表的埋藏环境风化效应更强连相对耐风化的高钾玻璃都被全部风化。这里有个容易忽略的细节类型、纹饰、颜色、是否风化四个字段的关联不是只看交叉表就够的。论文把表单1的字段量化之后丢进SPSS做了相关分析和卡方检验结果表明“纹饰与类型”“颜色与类型”存在显著关联而“纹饰与是否风化”的显著性较弱。这种相关性检验加显著性判断的做法比肉眼观察交叉表靠谱得多值得在类似竞赛场景里沿用。3. 问题一拆解PCA与因子分析反推风化前成分3.1 为什么不能直接线性外推问题一的原始表述是根据风化点检测数据预测其风化前的化学成分含量。最朴素的想法是把风化样本与未风化样本放一起做线性回归但论文选择了因子分析原因在于风化是一个持续叠加的过程不同埋藏环境下各种化学成分的流失与富集速度不同观测到的成分是“原始成分经过多个因子线性变换”之后的结果。直接外推忽略各成分之间的协方差结构预测结果很难落在合理范围内。因子分析模型写成这样X μ AF ε其中X是观测到的风化后成分向量μ是期望向量F是公共因子向量对应未风化时的成分构成A是负载矩阵ε是特殊因子向量。这个式子的含义是风化后的每个成分都可以看成若干公共因子按不同权重线性叠加的结果。反过来只要估计出A并求逆就能用风化后的X反推F得到未风化时的成分。3.2 从相关系数矩阵特征分解得到负载矩阵求负载矩阵A的常见做法是先对风化样本的相关系数矩阵做特征分解得到特征值λ1≥λ2≥…≥λm和标准正交特征向量u1、u2、…、um然后取前k个主成分构造A [√λ1·u1, √λ2·u2, …, √λk·uk]特征值代表对应方向上的方差贡献√λ乘上特征向量等于把主成分方向还原成原始量纲的因子载荷。论文对“铅钡玻璃-C纹”样本取前8个主成分累计贡献率超过97%说明8个因子基本能还原原始成分的几乎全部信息。用Python实现如下import numpy as np import pandas as pd df pd.read_excel(data/lead_barium_C_weathered.xlsx) cols df.columns[1:13] # 前12个氧化物成分列 X df[cols].values # 标准化 X_std (X - X.mean(axis0)) / X.std(axis0, ddof0) # 相关系数矩阵的特征分解 R np.corrcoef(X_std.T) eigvals, eigvecs np.linalg.eigh(R) idx np.argsort(eigvals)[::-1] eigvals, eigvecs eigvals[idx], eigvecs[:, idx] # 取累计贡献率超过97%的主成分数量 cum_ratio np.cumsum(eigvals) / eigvals.sum() k int(np.argmax(cum_ratio 0.97)) 1 # 负载矩阵A A eigvecs[:, :k] * np.sqrt(eigvals[:k])注意np.linalg.eigh返回的特征值默认升序必须按降序重排否则取前k个主成分会拿错方向。A的shape是12×k每一列是一个因子的负载向量绝对值越大说明该成分对这个因子越敏感。这里手动做了标准化也可以用sklearn.preprocessing.StandardScaler替代效果一样但要注意ddof参数对标准差的影响。3.3 转换矩阵求逆与未风化成分预测有了负载矩阵A用最小二乘意义下的伪逆求转换矩阵把公共因子F表示为风化后成分X的线性组合F A⁺(X − μ)代入风化样本得到F就是剥离风化效应之后的预测值。论文给铅钡玻璃-C纹的“严重风化点”预测出SiO2为51.62%明显高于风化样本实测的36%上下与未风化样本的分布区间吻合说明因子分析的反推方向是正确的。# 求伪逆并反推公共因子 A_pinv np.linalg.pinv(A) F A_pinv (X_std - X_std.mean(axis0)).T # 还原到原始成分量纲 X_recover F.T A_pinv.T * X.std(axis0, ddof0) X.mean(axis0)这里有一个容易翻车的细节标准化时记录的均值和标准差反推时先减均值再乘伪逆最后要乘回标准差并加回均值顺序反了预测值会整体偏移。我一般会把还原结果和未风化样本的均值做一个对比散点图一眼就能看出预测是否落在合理区间比看单个数字靠谱。3.4 先PCA再因子分析的设计意图论文在问题一里的顺序是先做主成分分析降维再在低维空间做因子分析。从工程角度看有两个好处一是14个化学成分里不少变量高度相关比如PbO和BaO经常同向变化直接求逆矩阵会病态二是降维之后只保留累计贡献率高的主成分噪声和与风化无关的微量成分被过滤掉预测结果更稳定。这个设计也提醒我们一件事因子分析不是上来就对原始数据跑先做一次PCA看方差结构再决定保留几个因子是更稳妥的建模顺序。很多同学直接把原始成分丢进因子分析得到的结果往往被少数高方差成分主导反而不如先降维再分解的做法干净。4. 分类建模链路区间划分、亚类聚类与灰色关联度4.1 四类样本的主成分区间划分问题二的目标是找到高钾玻璃、铅钡玻璃的分类规律并做亚类划分。论文的思路很工程化先把表单2拆成四类——铅钡风化、铅钡未风化、高钾风化、高钾未风化对每一类单独做主成分分析选出累计贡献率超过85%的主要化学成分然后“去掉最大值和最小值之后将剩下的最大值和最小值形成一个区间”。以铅钡未风化为例选出的6个成分是SiO2、Na2O、K2O、CaO、MgO、Al2O3贡献率88.89%对应区间分别是[50.61, 69.71]、[0.8, 5.74]、[0.15, 0.32]、[0.38, 2.82]、[0.71, 1.49]、[1.44, 13.65]。高钾未风化选5个成分贡献率88.89%铅钡风化选5个成分贡献率87.36%高钾风化只选3个成分贡献率95.41%。这个区间法的本质是一个非参数分类器不假设样本服从正态分布只用去除极值后的覆盖区间描述一类样本的分布范围。好处是简单直观、可解释性强坏处是区间之间可能重叠边界样本判断不稳所以论文在问题三里专门做了敏感性分析来验证。实现上用排序取中间值即可def build_interval(data): s np.sort(data) return s[1], s[-2]这里data是一类样本在某个化学成分上的取值数组去掉一个最大值和一个最小值后取新极值作为区间下界和上界。样本量较小的时候这种去极值操作对区间宽度影响很大换数据集时要注意样本量是否足够支撑这个操作。4.2 最短距离法做亚类划分与敏感性分析论文选出SiO2作为差异最大的成分做亚类划分用最短距离法单链接聚类把每个大类的SiO2聚成两类得到两个聚点再取两个聚点的平均值作为分割点划分成两个区间。例如铅钡未风化的聚点是64.82和52.49分割点取中间值58.70左右得到[50.61, 58.70)和[58.70, 69.71]两个亚类区间。最短距离法的特点是类与类的距离取“最近两点的距离”计算简单聚类速度快但容易受离群点影响形成链式结构。在这道题里样本量不大SiO2区分度明显用最短距离法结果合理。用scipy实现from scipy.cluster.hierarchy import fcluster, linkage from scipy.spatial.distance import pdist def cluster_split(values, n_clusters2): X np.array(values, dtypefloat).reshape(-1, 1) Z linkage(pdist(X), methodsingle) labels fcluster(Z, n_clusters, criterionmaxclust) centers [X[labels i].mean() for i in range(1, n_clusters 1)] return centers参数values传入某个大类所有样本的SiO2含量n_clusters指定分成几类。methodsingle就是最短距离法pdist先计算样本两两距离linkage做层次聚类fcluster按最大聚类数切分返回每类的聚点。敏感性分析的做法是给SiO2数据加3%的随机扰动再走一遍亚类划分统计样本落入第一类或第二类的比例变化。论文发现三类数据的误差超过20%说明这种亚类划分对噪声比较敏感实际使用时不能把亚类标签当稳定结论来用。4.3 问题三的归属判别问题三相对简单把表单3里的未知样本成分按问题二求得的四类区间逐一比对数值落在哪个区间或最靠近哪个区间就判定为哪一类。论文将表单3样本判别为铅钡玻璃类并且做了3%随机扰动后判别结果不变说明未知样本落在区间内部较深的位置鲁棒性好。从工程角度区间判别可以写成一个查表函数def classify_by_intervals(row, intervals): best_label, best_score None, float(-inf) for label, comp_intervals in intervals.items(): score 0 for comp, (lo, hi) in comp_intervals.items(): if lo row[comp] hi: score 1 else: score - min(abs(row[comp] - lo), abs(row[comp] - hi)) if score best_score: best_score, best_label score, label return best_label这里的row是未知样本的成分字典intervals是四类样本的成分区间字典。评分逻辑不是简单的命中计数而是“命中一个成分加1未命中按到边界的距离扣分”能有效区分“所有成分都在区间内”和“恰好压线”两种情形对边界样本的判断更合理。4.4 灰色关联度分析成分关联关系问题四要求分析不同类别玻璃化学成分之间的关联关系。论文用的是灰色关联度分析五步走确定参考数列与比较数列、无量纲化、计算关联系数、计算关联度、按关联度排序。以铅钡未风化数据为例15个样品作为参考数列6个主要化学成分作为比较数列。先做均值化无量纲化再按公式计算每个成分与参考数列的关联系数取平均得到关联度。分辨系数ρ取0.5这是灰色关联度的工程默认值ρ越小分辨力越大但抗干扰能力下降。Python实现def grey_relational_grade(X, rho0.5): # 均值化无量纲化每个特征除以自身均值 X_norm X / X.mean(axis0, keepdimsTrue) # 参考数列每个样本在所有特征上的均值 ref X_norm.mean(axis1) delta np.abs(X_norm - ref[:, None]) delta_min, delta_max delta.min(), delta.max() gamma (delta_min rho * delta_max) / (delta rho * delta_max) return gamma.mean(axis0)参数X的shape是(n_samples, n_features)每一列是一个化学成分。rho是分辨系数默认0.5可以保证分辨力与抗干扰能力的平衡。最终结论是铅钡玻璃的关联度排序为SiO2CaOK2OAl2O3MgONa2O高钾玻璃为SiO2Na2OK2OCaOMgO。两类玻璃中SiO2的关联性都是最强但第二、第三名的成分明显不同铅钡玻璃里CaO、K2O靠前高钾玻璃里Na2O、K2O靠前反映了两类玻璃助熔剂体系的差异。论文同时指出四类样本里SiO2的关联性始终最强其他成分的关联度差异不算大这个结论也可以作为后续分析的重点方向。提示灰色关联度的ρ一般取0.5当数据噪声较大时可以适当提高ρ到0.6左右来平滑关联系数但分辨力会随之下降需要根据数据情况权衡。5. 避坑指南复现这套过程最容易踩的五个坑5.1 OCR乱码导致的成分字段错位现象拿到的文档里很多化学式显示为乱码占位符复制到Excel后跑PCA主成分贡献率全乱特征值分布异常。原因原始文档在格式转换或复制粘贴时带上下标的化学式丢失了Unicode映射导致列错位、数值串行。解决所有成分数据以原始Excel表单为准不要从论文正文里复制数字。对每个样本做合计校验偏离100%超过5%的要逐行定位是哪一列串位。宁可丢掉个别异常样本也不要带病数据进模型。5.2 风化系数公式用反现象复现者想“统一公式”把两类玻璃的风化系数都定义成同一个表达式结果高钾玻璃的判别阈值怎么调都分不开样本。原因这两类玻璃的风化机理正好相反一个SiO2流失、一个SiO2富集共用公式会让风化方向互相抵消。解决严格按论文分开定义铅钡玻璃用(PbOBaO)/SiO2高钾玻璃用SiO2/(K2OCaOMgOAl2O3Fe2O3CuO)。换数据集时要重新画排序折线图找断层不要沿用论文里的阈值。5.3 直接把风化样本和未风化样本混在一起做PCA现象主成分载荷全部集中在前两个成分上各类样本的区间严重重叠分类几乎失效。原因风化样本与未风化样本的均值差异太大PCA会把“风化vs未风化”当成第一主成分淹没类型之间的真实差异。解决先按“类型是否风化”拆成四组每组单独做PCA再比较各组的主成分区间。这才是论文采用的正确做法也是复现时最容易忽略的一步。5.4 灰色关联度无量纲化方式不统一现象同一份数据用初值化和均值化两种方式处理算出的关联度排序不一致甚至结论相反。原因初值化适合有趋势的时间序列数据均值化适合水平型数据。化学成分含量相差几十倍初值化会让含量小的成分被忽略。解决化学成分这类“考察绝对偏差、无时间趋势”的数据用均值化处理并且在文档里写清楚用的是什么方法保证结论可复现。5.5 随机扰动不固定随机种子现象别人复现敏感性分析时3%扰动后的类别变化结论总是对不上同一份数据两次跑结果不同。原因随机数生成器没有固定种子每次生成的扰动序列不同边界样本在区间两侧跳来跳去。解决敏感性分析前先固定np.random.seed(42)并把扰动后的数据存一份副本结论才能稳定复现。这是很多人忽略但影响极大的细节。6. 二次开发把区间分类法提成可用的玻璃鉴别函数库这套资源里的建模过程虽然是围绕竞赛题展开的但核心算法——风化系数、PCA区间分类、最短距离亚类、灰色关联度——都是独立可复用的模块。如果要把它们用到真实文物数据或其他类似的风化成分研究我会改造成两层结构的函数库。第一层是通用算法函数层比如pca_interval_selection()、single_linkage_subclass()、grey_relational_grade()这些函数只接受numpy数组和参数不关心数据来自哪张表单第二层是领域适配层把表单列名映射成语义字段例如把“二氧化硅(SiO2)”统一为componentSiO2把“是否风化”统一为labelweathered。这样换一批新样本时只需要改适配层算法层完全不用动。验证方法上我一般会在区间分类器之外加一个留一法交叉验证每次留出一个样本用其余样本重新构造区间看被留出的样本能否被正确分类。论文里的敏感性分析回答的是“数据扰动后稳不稳”留一法回答的是“样本依赖性强不强”两者互补。如果表单3里的样本存在已知真值标签这是最优先的验证手段。具体可改的点有三处一是把硬编码的区间改成数据库或配置文件新增样本时自动重算区间而不是手动改代码二是把“命中区间加分、未命中扣距离分”的判别函数升级成加权投票让多个成分的贡献更均衡三是灰色关联度的分辨系数ρ不固定为0.5改成按数据噪声水平做网格搜索选择排序结果最稳定的值。我在改造这些代码时的教训是区间法这种朴素分类器最大的风险在边界处样本训练数据一变整个区间就要重新算。从那以后我每次拆这类“论文代码”资源时都强制走一遍“原始表单→合计校验→分组PCA→区间重算”的完整流程确保拿到的是可复现的结果而不是一个黑盒输出。这套代码本身门槛不高但它把PCA、聚类、灰色关联度这些散装算法串成了一条完整的问题解决链路值得对着数据跑一遍再动手改。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/17 19:14:44 · 阅读 1.2k 数学建模C题：主成分分析与聚类破解古代玻璃风化之谜
简介2022年全国大学生数学建模竞赛C题的完整解题资料包面向备赛学生与数学建模爱好者围绕古代玻璃文物成分分析与鉴别问题系统梳理了数据预处理、数据探索、特征工程、主成分分析与因子分析等关键建模方法。资源为单个PDF文档大小3.55MB虽以文档形式交付但内容涵盖完整建模思路、核心代码实现、结果图表与结论便于对照学习。包内重点展示高钾玻璃与铅钡玻璃的风化系数定义、亚类划分区间、敏感性分析与灰色关联度分析过程并给出四问的详细解答与模型解释适合需要参考优秀论文结构、理解降维与关联分析应用的读者目前已有927人学习下载。整体内容紧凑精炼可直接作为赛题复盘、论文写作或算法实战的参考资料。1. 项目速览数学建模C题如何给古代玻璃“去风化”2022年全国大学生数学建模C题的核心是给一批风化了几百上千年的古代玻璃文物做“成分还原”和“自动鉴别”。这道题乍一看是考古学问题实际拆开全是数据活儿三张表单、十几个化学成分指标、一堆风化与未风化的标记要用主成分分析、因子分析、最短距离聚类、灰色关联度把这堆数字变成可以解释的结论。这个资源就是整套论文与可运行代码的打包论文负责讲清建模思路代码负责把四个问题的结果完整复现。适合两类人一是备战数学建模竞赛、想找一份完整C题范式的同学二是对成分数据做降维、聚类、敏感性分析套路感兴趣的开发者。2. 数据底子与预处理三张表单与两类风化系数2.1 三张表单的字段关系与合并思路这道题的数据核心是三张表先把结构理清后面所有模型都建立在这三张表之上表单1文物采样点的档案包含样品编号、玻璃类型高钾玻璃/铅钡玻璃、是否风化、颜色、纹饰表单2每个采样点的化学成分检测结果从B列到O列共十几个成分包括二氧化硅SiO2、氧化钠Na2O、氧化钾K2O、氧化钙CaO、氧化镁MgO、氧化铝Al2O3、氧化铁Fe2O3、氧化铜CuO、氧化铅PbO、氧化钡BaO、五氧化二磷P2O5、氧化锶SrO等表单3一批未知类别的玻璃文物成分数据待判别所属类型。原始文档里用Excel做了大量预处理在B列到O列后面增加一个合计列P验证每行成分百分比加起来是否接近100%从A列抽取样品编号放入Q列再把表单1中的“玻璃类型、是否风化、颜色、纹饰”四项信息按样品编号关联到表单2形成R、S、T、U四列。这一步本质上是两张表按主键做left join用pandas实现是更常见的做法import pandas as pd info pd.read_excel(表单1.xlsx) chem pd.read_excel(表单2.xlsx) # 从A列抽出纯数字编号兼容“XX部位1”这类复合编号 chem[sample_id] chem[A列].astype(str).str.extract(r(\d)) # 计算各成分合计做数据质量校验 chem[total] chem.iloc[:, 1:15].sum(axis1) merged chem.merge(info, left_onsample_id, right_on样品编号, howleft)这里的sample_id提取是关键表单2里同一件文物会有“部位1”“部位2”“严重风化点”等多个采样点直接用原始文本做关联会失败先抽编号再关联是这步的通用做法。total列用来检查成分合计是否在100%附近正常玻璃样本的检测值合计应在100%上下偏差超过5%的样本多半存在复制时的乱码或串列问题要单独标记出来。howleft保证化学表里每一行的采样数据都保留下来避免有效样本被内连接悄悄过滤掉。2.2 两类风化系数的定义与阈值确定预处理的关键一步是定义风化系数。论文没有用统一公式而是给两类玻璃分别定义因为这两类玻璃的风化规律恰好相反铅钡玻璃风化系数 (PbO BaO) / SiO2 高钾玻璃风化系数 SiO2 / (K2O CaO MgO Al2O3 Fe2O3 CuO)从物理意义上看铅钡玻璃风化后SiO2占比显著流失氧化铅、氧化钡相对富集比值变大高钾玻璃风化后SiO2占比显著提升钾、钙、镁、铝等助熔剂氧化物流失比值同样变大。两边风化系数都是“越大越风化”但分子分母的含义相反这是这个设计最精巧的地方。再看阈值是怎么定的。论文对“铅钡玻璃-A纹”这组样本按风化系数排序发现0.879到1.3之间存在明显断层于是把1.0作为判别阈值对高钾玻璃7.23到12.41之间存在断层取10.0作为阈值。这种取法非常直观本质就是找排序后相邻值的最大空档用代码可以写成import numpy as np def split_by_gap(values): s np.sort(values) gaps np.diff(s) idx np.argmax(gaps) threshold (s[idx] s[idx1]) / 2 return threshold, s, gaps参数values传入同一类型样本的风化系数数组返回的是断层中点作为判别阈值以及排序后的序列和相邻差值。用这个函数对三类样本分别计算得到的结果与论文里手工观察的阈值一致。实际使用中我一般会先画一张排序后的折线图肉眼看断层在哪儿再用这个函数辅助定位避免把正常波动误判成阈值点。2.3 数据探索挖出的三条规律对预处理后的数据做分组透视原文得出三个值得记住的观察纹饰B的玻璃文物全部是高钾玻璃、全部风化、颜色全是蓝绿色纹饰A的高钾玻璃全部无风化且为蓝绿色铅钡玻璃的风化与无风化比例是11:5纹饰C的高钾玻璃全部无风化铅钡玻璃的风化与无风化比例是17:73。由此推出两个结论同等条件同纹饰意味着年代和地域可能相近下铅钡玻璃比高钾玻璃更容易风化纹饰B代表的埋藏环境风化效应更强连相对耐风化的高钾玻璃都被全部风化。这里有个容易忽略的细节类型、纹饰、颜色、是否风化四个字段的关联不是只看交叉表就够的。论文把表单1的字段量化之后丢进SPSS做了相关分析和卡方检验结果表明“纹饰与类型”“颜色与类型”存在显著关联而“纹饰与是否风化”的显著性较弱。这种相关性检验加显著性判断的做法比肉眼观察交叉表靠谱得多值得在类似竞赛场景里沿用。3. 问题一拆解PCA与因子分析反推风化前成分3.1 为什么不能直接线性外推问题一的原始表述是根据风化点检测数据预测其风化前的化学成分含量。最朴素的想法是把风化样本与未风化样本放一起做线性回归但论文选择了因子分析原因在于风化是一个持续叠加的过程不同埋藏环境下各种化学成分的流失与富集速度不同观测到的成分是“原始成分经过多个因子线性变换”之后的结果。直接外推忽略各成分之间的协方差结构预测结果很难落在合理范围内。因子分析模型写成这样X μ AF ε其中X是观测到的风化后成分向量μ是期望向量F是公共因子向量对应未风化时的成分构成A是负载矩阵ε是特殊因子向量。这个式子的含义是风化后的每个成分都可以看成若干公共因子按不同权重线性叠加的结果。反过来只要估计出A并求逆就能用风化后的X反推F得到未风化时的成分。3.2 从相关系数矩阵特征分解得到负载矩阵求负载矩阵A的常见做法是先对风化样本的相关系数矩阵做特征分解得到特征值λ1≥λ2≥…≥λm和标准正交特征向量u1、u2、…、um然后取前k个主成分构造A [√λ1·u1, √λ2·u2, …, √λk·uk]特征值代表对应方向上的方差贡献√λ乘上特征向量等于把主成分方向还原成原始量纲的因子载荷。论文对“铅钡玻璃-C纹”样本取前8个主成分累计贡献率超过97%说明8个因子基本能还原原始成分的几乎全部信息。用Python实现如下import numpy as np import pandas as pd df pd.read_excel(data/lead_barium_C_weathered.xlsx) cols df.columns[1:13] # 前12个氧化物成分列 X df[cols].values # 标准化 X_std (X - X.mean(axis0)) / X.std(axis0, ddof0) # 相关系数矩阵的特征分解 R np.corrcoef(X_std.T) eigvals, eigvecs np.linalg.eigh(R) idx np.argsort(eigvals)[::-1] eigvals, eigvecs eigvals[idx], eigvecs[:, idx] # 取累计贡献率超过97%的主成分数量 cum_ratio np.cumsum(eigvals) / eigvals.sum() k int(np.argmax(cum_ratio 0.97)) 1 # 负载矩阵A A eigvecs[:, :k] * np.sqrt(eigvals[:k])注意np.linalg.eigh返回的特征值默认升序必须按降序重排否则取前k个主成分会拿错方向。A的shape是12×k每一列是一个因子的负载向量绝对值越大说明该成分对这个因子越敏感。这里手动做了标准化也可以用sklearn.preprocessing.StandardScaler替代效果一样但要注意ddof参数对标准差的影响。3.3 转换矩阵求逆与未风化成分预测有了负载矩阵A用最小二乘意义下的伪逆求转换矩阵把公共因子F表示为风化后成分X的线性组合F A⁺(X − μ)代入风化样本得到F就是剥离风化效应之后的预测值。论文给铅钡玻璃-C纹的“严重风化点”预测出SiO2为51.62%明显高于风化样本实测的36%上下与未风化样本的分布区间吻合说明因子分析的反推方向是正确的。# 求伪逆并反推公共因子 A_pinv np.linalg.pinv(A) F A_pinv (X_std - X_std.mean(axis0)).T # 还原到原始成分量纲 X_recover F.T A_pinv.T * X.std(axis0, ddof0) X.mean(axis0)这里有一个容易翻车的细节标准化时记录的均值和标准差反推时先减均值再乘伪逆最后要乘回标准差并加回均值顺序反了预测值会整体偏移。我一般会把还原结果和未风化样本的均值做一个对比散点图一眼就能看出预测是否落在合理区间比看单个数字靠谱。3.4 先PCA再因子分析的设计意图论文在问题一里的顺序是先做主成分分析降维再在低维空间做因子分析。从工程角度看有两个好处一是14个化学成分里不少变量高度相关比如PbO和BaO经常同向变化直接求逆矩阵会病态二是降维之后只保留累计贡献率高的主成分噪声和与风化无关的微量成分被过滤掉预测结果更稳定。这个设计也提醒我们一件事因子分析不是上来就对原始数据跑先做一次PCA看方差结构再决定保留几个因子是更稳妥的建模顺序。很多同学直接把原始成分丢进因子分析得到的结果往往被少数高方差成分主导反而不如先降维再分解的做法干净。4. 分类建模链路区间划分、亚类聚类与灰色关联度4.1 四类样本的主成分区间划分问题二的目标是找到高钾玻璃、铅钡玻璃的分类规律并做亚类划分。论文的思路很工程化先把表单2拆成四类——铅钡风化、铅钡未风化、高钾风化、高钾未风化对每一类单独做主成分分析选出累计贡献率超过85%的主要化学成分然后“去掉最大值和最小值之后将剩下的最大值和最小值形成一个区间”。以铅钡未风化为例选出的6个成分是SiO2、Na2O、K2O、CaO、MgO、Al2O3贡献率88.89%对应区间分别是[50.61, 69.71]、[0.8, 5.74]、[0.15, 0.32]、[0.38, 2.82]、[0.71, 1.49]、[1.44, 13.65]。高钾未风化选5个成分贡献率88.89%铅钡风化选5个成分贡献率87.36%高钾风化只选3个成分贡献率95.41%。这个区间法的本质是一个非参数分类器不假设样本服从正态分布只用去除极值后的覆盖区间描述一类样本的分布范围。好处是简单直观、可解释性强坏处是区间之间可能重叠边界样本判断不稳所以论文在问题三里专门做了敏感性分析来验证。实现上用排序取中间值即可def build_interval(data): s np.sort(data) return s[1], s[-2]这里data是一类样本在某个化学成分上的取值数组去掉一个最大值和一个最小值后取新极值作为区间下界和上界。样本量较小的时候这种去极值操作对区间宽度影响很大换数据集时要注意样本量是否足够支撑这个操作。4.2 最短距离法做亚类划分与敏感性分析论文选出SiO2作为差异最大的成分做亚类划分用最短距离法单链接聚类把每个大类的SiO2聚成两类得到两个聚点再取两个聚点的平均值作为分割点划分成两个区间。例如铅钡未风化的聚点是64.82和52.49分割点取中间值58.70左右得到[50.61, 58.70)和[58.70, 69.71]两个亚类区间。最短距离法的特点是类与类的距离取“最近两点的距离”计算简单聚类速度快但容易受离群点影响形成链式结构。在这道题里样本量不大SiO2区分度明显用最短距离法结果合理。用scipy实现from scipy.cluster.hierarchy import fcluster, linkage from scipy.spatial.distance import pdist def cluster_split(values, n_clusters2): X np.array(values, dtypefloat).reshape(-1, 1) Z linkage(pdist(X), methodsingle) labels fcluster(Z, n_clusters, criterionmaxclust) centers [X[labels i].mean() for i in range(1, n_clusters 1)] return centers参数values传入某个大类所有样本的SiO2含量n_clusters指定分成几类。methodsingle就是最短距离法pdist先计算样本两两距离linkage做层次聚类fcluster按最大聚类数切分返回每类的聚点。敏感性分析的做法是给SiO2数据加3%的随机扰动再走一遍亚类划分统计样本落入第一类或第二类的比例变化。论文发现三类数据的误差超过20%说明这种亚类划分对噪声比较敏感实际使用时不能把亚类标签当稳定结论来用。4.3 问题三的归属判别问题三相对简单把表单3里的未知样本成分按问题二求得的四类区间逐一比对数值落在哪个区间或最靠近哪个区间就判定为哪一类。论文将表单3样本判别为铅钡玻璃类并且做了3%随机扰动后判别结果不变说明未知样本落在区间内部较深的位置鲁棒性好。从工程角度区间判别可以写成一个查表函数def classify_by_intervals(row, intervals): best_label, best_score None, float(-inf) for label, comp_intervals in intervals.items(): score 0 for comp, (lo, hi) in comp_intervals.items(): if lo row[comp] hi: score 1 else: score - min(abs(row[comp] - lo), abs(row[comp] - hi)) if score best_score: best_score, best_label score, label return best_label这里的row是未知样本的成分字典intervals是四类样本的成分区间字典。评分逻辑不是简单的命中计数而是“命中一个成分加1未命中按到边界的距离扣分”能有效区分“所有成分都在区间内”和“恰好压线”两种情形对边界样本的判断更合理。4.4 灰色关联度分析成分关联关系问题四要求分析不同类别玻璃化学成分之间的关联关系。论文用的是灰色关联度分析五步走确定参考数列与比较数列、无量纲化、计算关联系数、计算关联度、按关联度排序。以铅钡未风化数据为例15个样品作为参考数列6个主要化学成分作为比较数列。先做均值化无量纲化再按公式计算每个成分与参考数列的关联系数取平均得到关联度。分辨系数ρ取0.5这是灰色关联度的工程默认值ρ越小分辨力越大但抗干扰能力下降。Python实现def grey_relational_grade(X, rho0.5): # 均值化无量纲化每个特征除以自身均值 X_norm X / X.mean(axis0, keepdimsTrue) # 参考数列每个样本在所有特征上的均值 ref X_norm.mean(axis1) delta np.abs(X_norm - ref[:, None]) delta_min, delta_max delta.min(), delta.max() gamma (delta_min rho * delta_max) / (delta rho * delta_max) return gamma.mean(axis0)参数X的shape是(n_samples, n_features)每一列是一个化学成分。rho是分辨系数默认0.5可以保证分辨力与抗干扰能力的平衡。最终结论是铅钡玻璃的关联度排序为SiO2CaOK2OAl2O3MgONa2O高钾玻璃为SiO2Na2OK2OCaOMgO。两类玻璃中SiO2的关联性都是最强但第二、第三名的成分明显不同铅钡玻璃里CaO、K2O靠前高钾玻璃里Na2O、K2O靠前反映了两类玻璃助熔剂体系的差异。论文同时指出四类样本里SiO2的关联性始终最强其他成分的关联度差异不算大这个结论也可以作为后续分析的重点方向。提示灰色关联度的ρ一般取0.5当数据噪声较大时可以适当提高ρ到0.6左右来平滑关联系数但分辨力会随之下降需要根据数据情况权衡。5. 避坑指南复现这套过程最容易踩的五个坑5.1 OCR乱码导致的成分字段错位现象拿到的文档里很多化学式显示为乱码占位符复制到Excel后跑PCA主成分贡献率全乱特征值分布异常。原因原始文档在格式转换或复制粘贴时带上下标的化学式丢失了Unicode映射导致列错位、数值串行。解决所有成分数据以原始Excel表单为准不要从论文正文里复制数字。对每个样本做合计校验偏离100%超过5%的要逐行定位是哪一列串位。宁可丢掉个别异常样本也不要带病数据进模型。5.2 风化系数公式用反现象复现者想“统一公式”把两类玻璃的风化系数都定义成同一个表达式结果高钾玻璃的判别阈值怎么调都分不开样本。原因这两类玻璃的风化机理正好相反一个SiO2流失、一个SiO2富集共用公式会让风化方向互相抵消。解决严格按论文分开定义铅钡玻璃用(PbOBaO)/SiO2高钾玻璃用SiO2/(K2OCaOMgOAl2O3Fe2O3CuO)。换数据集时要重新画排序折线图找断层不要沿用论文里的阈值。5.3 直接把风化样本和未风化样本混在一起做PCA现象主成分载荷全部集中在前两个成分上各类样本的区间严重重叠分类几乎失效。原因风化样本与未风化样本的均值差异太大PCA会把“风化vs未风化”当成第一主成分淹没类型之间的真实差异。解决先按“类型是否风化”拆成四组每组单独做PCA再比较各组的主成分区间。这才是论文采用的正确做法也是复现时最容易忽略的一步。5.4 灰色关联度无量纲化方式不统一现象同一份数据用初值化和均值化两种方式处理算出的关联度排序不一致甚至结论相反。原因初值化适合有趋势的时间序列数据均值化适合水平型数据。化学成分含量相差几十倍初值化会让含量小的成分被忽略。解决化学成分这类“考察绝对偏差、无时间趋势”的数据用均值化处理并且在文档里写清楚用的是什么方法保证结论可复现。5.5 随机扰动不固定随机种子现象别人复现敏感性分析时3%扰动后的类别变化结论总是对不上同一份数据两次跑结果不同。原因随机数生成器没有固定种子每次生成的扰动序列不同边界样本在区间两侧跳来跳去。解决敏感性分析前先固定np.random.seed(42)并把扰动后的数据存一份副本结论才能稳定复现。这是很多人忽略但影响极大的细节。6. 二次开发把区间分类法提成可用的玻璃鉴别函数库这套资源里的建模过程虽然是围绕竞赛题展开的但核心算法——风化系数、PCA区间分类、最短距离亚类、灰色关联度——都是独立可复用的模块。如果要把它们用到真实文物数据或其他类似的风化成分研究我会改造成两层结构的函数库。第一层是通用算法函数层比如pca_interval_selection()、single_linkage_subclass()、grey_relational_grade()这些函数只接受numpy数组和参数不关心数据来自哪张表单第二层是领域适配层把表单列名映射成语义字段例如把“二氧化硅(SiO2)”统一为componentSiO2把“是否风化”统一为labelweathered。这样换一批新样本时只需要改适配层算法层完全不用动。验证方法上我一般会在区间分类器之外加一个留一法交叉验证每次留出一个样本用其余样本重新构造区间看被留出的样本能否被正确分类。论文里的敏感性分析回答的是“数据扰动后稳不稳”留一法回答的是“样本依赖性强不强”两者互补。如果表单3里的样本存在已知真值标签这是最优先的验证手段。具体可改的点有三处一是把硬编码的区间改成数据库或配置文件新增样本时自动重算区间而不是手动改代码二是把“命中区间加分、未命中扣距离分”的判别函数升级成加权投票让多个成分的贡献更均衡三是灰色关联度的分辨系数ρ不固定为0.5改成按数据噪声水平做网格搜索选择排序结果最稳定的值。我在改造这些代码时的教训是区间法这种朴素分类器最大的风险在边界处样本训练数据一变整个区间就要重新算。从那以后我每次拆这类“论文代码”资源时都强制走一遍“原始表单→合计校验→分组PCA→区间重算”的完整流程确保拿到的是可复现的结果而不是一个黑盒输出。这套代码本身门槛不高但它把PCA、聚类、灰色关联度这些散装算法串成了一条完整的问题解决链路值得对着数据跑一遍再动手改。希望帮到你。本文还有配套的精品资源点击获取